proteus.interior_struct.zalmoxis¶

def get_zalmoxis_output_filepath(outdir: str):
    """Returns the output file path for Zalmoxis data.
    Args:
        outdir (str): Output directory.
    Returns:
        str: Path to the output file.
    """
    return os.path.join(outdir, 'data', 'zalmoxis_output.dat')

def validate_zalmoxis_output_schema(
    output_path: str,
    hf_row: dict,
    rtol_radius: float = 1e-6,
    rtol_mass: float = 5e-2,
    mantle_mass_ref: float | None = None,
) -> None:
    """Verify zalmoxis_output.dat is consistent with hf_row scalars.

    The file is the contract Aragog reads inside ``solver.reset()``
    (eos_method=2). This check confirms the file's last r matches
    ``hf_row['R_int']``, and that the mantle mass matches
    ``hf_row['M_int'] - hf_row['M_core']``.
    Catches file I/O corruption, column-order mistakes, truncation,
    encoding drift, and Aragog/Zalmoxis schema desync at the file
    handover boundary.

    The mantle mass is taken from ``mantle_mass_ref`` when the caller supplies
    it: the structure's accumulator total ``mass_enclosed[-1]`` minus the exact
    core-mass target ``cmb_mass``. Re-integrating the coarse output nodes with a
    grid trapezoid instead diverges from the sub-grid-substepped RK45 integral
    across the steep interior density profile, reaching ~10% at high planet mass,
    which would false-reject a structure that actually conserves mass. With the
    reference supplied, part (b) reduces to a check that the core-mantle split
    point is well resolved: the residual it measures is the CMB-node snap
    overshoot ``mass_enclosed[cmb_index] - cmb_mass`` relative to the mantle
    mass, a fraction of one boundary shell. File-density-column corruption is
    caught downstream by Aragog's EOS-vs-mesh consistency at ``solver.reset()``;
    the radius check (part (a)) catches truncation and column swaps directly.
    When no reference is supplied the check falls back to the grid-trapezoidal
    shell-sum.

    Parameters
    ----------
    output_path : str
        Path to the just-written ``zalmoxis_output.dat``.
    hf_row : dict
        PROTEUS hf_row holding the scalar truth (R_int, M_int, M_core).
    rtol_radius : float
        Relative tolerance for the top-of-mantle vs. R_int check.
        Default 1e-6: the file's last r and ``hf_row['R_int']`` come
        from the same variable in zalmoxis_solver, so equality is
        exact modulo float-string round-trip noise. Tight tolerance
        catches truncation, last-line corruption, and column-swap
        bugs at the bit level.
    rtol_mass : float
        Relative tolerance for the integrated mantle mass vs.
        ``M_int - M_core``. Default 5e-2 (5%) reflects two stacked
        sources of legitimate mismatch:
        (a) integrator-method difference: Zalmoxis' ``mass_enclosed``
        is the ODE state from RK45 with sub-grid substepping, while
        the schema check re-integrates via a grid-trapezoidal
        shell-sum. On stiff CHILI density profiles this shows
        about 0.8 to 2.0 %.
        (b) ``blend_mesh_files`` post-write modifies the file in
        place but does not update hf_row scalars. When blending fires
        with alpha < 1 (capping large R-shifts), the file's integrated
        mass can drift up to about 5% from the unblended
        hf_row['M_int'], so this check has to tolerate it.
        5e-2 keeps a >2x margin over the worst legitimate noise
        while still catching gross corruption (column swap,
        truncation, byte-flip) at >>5 %. **The tight mass-conservation
        contract (<0.1 %) lives in the wrapper-level mass-anchor
        check on hf_row['M_int'] / hf_row['M_int_target']**, not here.

    Raises
    ------
    RuntimeError
        On any violation (file unreadable, wrong shape, radius
        mismatch, mass mismatch). Caller (the wrapper's except-block)
        catches this and routes through the fall-back path.
    """
    try:
        data = np.loadtxt(output_path)
    except Exception as exc:
        raise RuntimeError(
            'zalmoxis_output.dat schema violation: could not reload '
            f'the just-written file ({output_path}): {exc}'
        )
    if data.size == 0 or data.ndim != 2 or data.shape[1] != 5:
        raise RuntimeError(
            'zalmoxis_output.dat schema violation: unexpected shape '
            f'{data.shape if data.size else "empty"} '
            f'(expected (N, 5)) for {output_path}'
        )

    r_file = data[:, 0]
    rho_file = data[:, 2]

    # (a) top-of-mantle radius == hf_row['R_int']
    r_top = float(r_file[-1])
    r_int_hf = float(hf_row.get('R_int', 0.0))
    if r_int_hf > 0:
        r_rel = abs(r_top / r_int_hf - 1.0)
        if r_rel > rtol_radius:
            raise RuntimeError(
                'zalmoxis_output.dat schema violation: top-of-mantle '
                f'r={r_top:.6e} from file differs from '
                f'hf_row[R_int]={r_int_hf:.6e} '
                f'(rel={r_rel:.3e} > {rtol_radius:.1e})'
            )

    # (b) mantle integrated mass == hf_row['M_int'] - hf_row['M_core'].
    # Prefer the structure's own RK45 accumulator mantle mass when the
    # caller supplies it (mantle_mass_ref). That value is the sub-grid
    # substepped ODE integral, exact to the solver tolerance. Re-integrating
    # the coarse output nodes with a grid trapezoid instead diverges from it
    # across the steep core-mantle density jump (the CMB-node snap attributes
    # a whole boundary shell to one side), reaching ~10% at high planet mass,
    # which would false-reject a structure that actually conserves mass. Only
    # when no reference is available do we fall back to the grid-trapezoidal
    # shell-sum (4/3 pi (r2^3 - r1^3) * rho_avg).
    if mantle_mass_ref is not None:
        mantle_mass = float(mantle_mass_ref)
        mass_source = 'accumulator'
    else:
        shells = (
            (4.0 / 3.0)
            * np.pi
            * (r_file[1:] ** 3 - r_file[:-1] ** 3)
            * 0.5
            * (rho_file[1:] + rho_file[:-1])
        )
        mantle_mass = float(np.sum(shells))
        mass_source = 'trapezoid'
    M_int_hf = float(hf_row.get('M_int', 0.0))
    M_core_hf = float(hf_row.get('M_core', 0.0))
    expected_mantle = M_int_hf - M_core_hf
    if expected_mantle > 0:
        m_rel = abs(mantle_mass / expected_mantle - 1.0)
        log.debug(
            'zalmoxis mass-closure: mantle split m_rel=%.3e (%s, rtol_mass=%.1e)',
            m_rel,
            mass_source,
            rtol_mass,
        )
        if m_rel > rtol_mass:
            raise RuntimeError(
                'zalmoxis_output.dat schema violation: '
                f'mantle mass ({mass_source})={mantle_mass:.6e} kg '
                'differs from hf_row[M_int - M_core]='
                f'{expected_mantle:.6e} kg '
                f'(rel={m_rel:.3e} > {rtol_mass:.1e})'
            )

def zalmoxis_mesh_gaps(output_path: str, hf_row: dict) -> tuple[float, float, float] | None:
    """Compare a mesh file's radial bounds with a helpfile row.

    The tolerance is that of Aragog's ``EntropySolver.reset()``,
    ``max(1 m, 1e-9 * (R_int - R_core))``, applied to both bounds in both
    directions; an 11-digit helpfile rounds the radii by at most 5e-5 m below
    1e7 m and 5e-4 m up to 1e8 m.

    Parameters
    ----------
    output_path : str
        Path to ``zalmoxis_output.dat``, its ``.prev`` backup or a saved copy.
    hf_row : dict
        Helpfile row holding ``R_core`` and ``R_int`` [m].

    Returns
    -------
    tuple of float or None
        ``(file r[0] - R_core, file r[-1] - R_int, tolerance)`` in metres, or
        None unless the file ends with a newline and holds at least two finite
        5-column rows with strictly increasing radii.
    """
    try:
        text = Path(output_path).read_text()
        data = np.loadtxt(text.splitlines(), ndmin=2) if text.strip() else None
    except (OSError, ValueError):
        return None
    if data is None:
        return None
    r = data[:, 0]
    if (
        not text.endswith('\n')
        or data.shape[1] != 5
        or r.size < 2
        or not np.isfinite(data).all()
        or (np.diff(r) <= 0).any()
    ):
        return None
    atol = max(1.0, 1.0e-9 * (hf_row['R_int'] - hf_row['R_core']))
    return r[0] - hf_row['R_core'], r[-1] - hf_row['R_int'], atol

def copy_zalmoxis_output(src: str, dst: str) -> None:
    """Copy through a temporary file, so ``dst`` is never partial if the process stops."""
    tmp = dst + '.tmp'
    try:
        shutil.copy2(src, tmp)
        os.replace(tmp, dst)
    except BaseException:
        Path(tmp).unlink(missing_ok=True)
        raise

def save_zalmoxis_output_snapshot(outdir: str, time: float) -> None:
    """Copy ``zalmoxis_output.dat`` to ``data/<time>_zalmoxis.dat``.

    The copy is written only next to the row's ``<time>_int.nc``, so archiving
    and pruning treat it as part of that snapshot.

    Parameters
    ----------
    outdir : str
        Run output directory.
    time : float
        Simulated time of the helpfile row being written [yr].
    """
    src = get_zalmoxis_output_filepath(outdir)
    data = os.path.join(outdir, 'data')
    if os.path.isfile(src) and os.path.isfile(snapshot_path_for_time(data, time, '_int.nc')):
        copy_zalmoxis_output(
            src, os.path.join(data, format_subyear_time(time) + '_zalmoxis.dat')
        )

def build_volatile_profile(hf_row: dict, mantle_eos: str):
    """Build a VolatileProfile from helpfile volatile masses.

    Computes per-phase (liquid/solid) mass fractions for dissolved volatiles
    that have Zalmoxis EOS tables. Returns None if no volatiles are dissolved
    or if the mantle liquid/solid masses are unavailable.

    Ownership of the liquid/solid split: the coupled PROTEUS path takes
    the per-phase masses from the outgassing chemistry (CALLIOPE/Aragog
    equilibrium via the helpfile), which already solved solubility.
    Zalmoxis's own ``partition_rule`` hook is a structure-side
    idealization for standalone Zalmoxis runs with no chemistry coupled;
    it is deliberately not used here.

    Parameters
    ----------
    hf_row : dict
        Current helpfile row with volatile mass keys (e.g. ``H2O_kg_liquid``,
        ``H2O_kg_solid``) and mantle mass keys: ``M_mantle`` with
        ``Phi_global`` (preferred, structural), or ``M_mantle_liquid`` and
        ``M_mantle_solid`` as a fallback.
    mantle_eos : str
        Primary mantle EOS identifier (e.g. ``'PALEOS:MgSiO3'``).

    Returns
    -------
    VolatileProfile or None
        Profile with per-phase fractions, or None if not applicable.
    """
    from zalmoxis.mixing import VolatileProfile

    # Denominators: the melt concentration must be taken against the same
    # reservoir the chemistry used and the structure carries. CALLIOPE
    # dissolves into M_mantle * Phi_global, and the structure integrates
    # the (structural) mantle mass, so split the structural M_mantle by
    # Phi_global. Aragog's M_mantle_liquid/M_mantle_solid use its PALEOS
    # EOS cell masses instead, ~7% denser than the wet structural mantle,
    # which diluted w_liquid and lost the same fraction of dissolved
    # water from the structure. Fall back to the per-phase masses when
    # the structural fields are absent (e.g. minimal test rows).
    M_mantle = float(hf_row.get('M_mantle', 0.0))
    Phi_global = float(hf_row.get('Phi_global', float('nan')))
    if M_mantle > 0 and np.isfinite(Phi_global):
        Phi_global = min(max(Phi_global, 0.0), 1.0)
        M_liq = Phi_global * M_mantle
        M_sol = (1.0 - Phi_global) * M_mantle
    else:
        M_liq = float(hf_row.get('M_mantle_liquid', 0.0))
        M_sol = float(hf_row.get('M_mantle_solid', 0.0))

    # Need mantle mass data to compute fractions
    if M_liq + M_sol <= 0:
        return None

    w_liquid = {}
    w_solid = {}
    has_nonzero = False

    for species, eos_name in VOLATILE_EOS_MAP.items():
        kg_liq = float(hf_row.get(f'{species}_kg_liquid', 0.0))
        kg_sol = float(hf_row.get(f'{species}_kg_solid', 0.0))

        # Only include species with meaningful dissolved mass
        if kg_liq + kg_sol <= 0:
            continue

        # Mass fraction in liquid phase
        w_l = kg_liq / M_liq if M_liq > 0 else 0.0
        # Mass fraction in solid phase
        w_s = kg_sol / M_sol if M_sol > 0 else 0.0

        w_liquid[eos_name] = w_l
        w_solid[eos_name] = w_s
        has_nonzero = True

    if not has_nonzero:
        return None

    # Normalize: total volatile fraction in each phase must not exceed 0.95
    # (at least 5% silicate). Clamp proportionally if sum exceeds the limit.
    for w_dict in (w_liquid, w_solid):
        total = sum(w_dict.values())
        max_volatile_frac = 0.95
        if total > max_volatile_frac:
            scale = max_volatile_frac / total
            for k in w_dict:
                w_dict[k] *= scale

    log.info(
        'Built VolatileProfile: liquid=%s, solid=%s',
        {k: f'{v:.4f}' for k, v in w_liquid.items()},
        {k: f'{v:.4f}' for k, v in w_solid.items()},
    )

    return VolatileProfile(
        w_liquid=w_liquid,
        w_solid=w_solid,
        primary_component=mantle_eos,
    )

def extend_mantle_eos_with_volatiles(mantle_eos: str, volatile_profile) -> str:
    """Extend a single-component mantle EOS string with volatile components.

    If a VolatileProfile is provided and the mantle EOS is a single component,
    this adds the volatile EOS components with small placeholder fractions.
    The actual fractions are overridden at each ODE step by the VolatileProfile.

    Parameters
    ----------
    mantle_eos : str
        Base mantle EOS string (e.g. ``'PALEOS:MgSiO3'``).
    volatile_profile : VolatileProfile or None
        Profile containing volatile EOS component names.

    Returns
    -------
    str
        Extended EOS string (e.g.
        ``'PALEOS:MgSiO3:0.98+PALEOS:H2O:0.01+Chabrier:H:0.01'``),
        or the original string if no extension needed.
    """
    if volatile_profile is None:
        return mantle_eos

    # Don't modify if already multi-component
    if '+' in mantle_eos:
        return mantle_eos

    # Collect all volatile EOS components from the profile
    all_vol_components = set()
    for d in (volatile_profile.w_liquid, volatile_profile.w_solid):
        all_vol_components.update(d.keys())

    if not all_vol_components:
        return mantle_eos

    # Build extended string with small placeholder fractions
    # (actual fractions set by VolatileProfile at each radius)
    n_vol = len(all_vol_components)
    placeholder = 0.01  # 1% each
    primary_frac = max(0.5, 1.0 - n_vol * placeholder)
    parts = [f'{mantle_eos}:{primary_frac:.4f}']
    for comp in sorted(all_vol_components):
        parts.append(f'{comp}:{placeholder:.4f}')

    extended = '+'.join(parts)
    log.info('Extended mantle EOS: %s -> %s', mantle_eos, extended)
    return extended

def solve_superliquidus_adiabat(config: Config, hf_row: dict | None) -> dict:
    """Memoised super-liquidus anchor solve, see ``_solve_superliquidus_adiabat``.

    A failure of the integration with ``ValueError``, ``KeyError``,
    ``IndexError``, ``ArithmeticError`` or ``RuntimeError``
    (``common.ANCHOR_NUMERICAL_ERRORS``) is raised as a chained
    ``InitialConditionError`` that names the type and P_cmb.
    ``NotImplementedError``, ``RecursionError`` and ``UnicodeError``
    (``common.ANCHOR_PASSTHROUGH_ERRORS``), and every other type such as
    ``TypeError``, ``AttributeError``, ``OSError`` or ``MemoryError``,
    propagate unchanged and are not memoised. An ``InitialConditionError``
    from the solve is raised as is. Both are memoised under the same key as
    a success, as a copy without traceback, so repeated structure solves at
    one P_cmb do not repeat a failing scan.

    Parameters
    ----------
    config : Config
        PROTEUS configuration.
    hf_row : dict or None
        Helpfile row; ``hf_row['P_cmb']`` is used when populated.

    Returns
    -------
    dict
        See ``_solve_superliquidus_adiabat``.

    Raises
    ------
    InitialConditionError
        See ``_solve_superliquidus_adiabat``; also for a numerical failure of
        the integration (chained). A memoised failure is raised again with
        the same message, chained to the stored copy.
    FileNotFoundError
        If an EOS or melting-curve file is missing.
    """
    from proteus.interior_energetics.common import (
        ANCHOR_NUMERICAL_ERRORS,
        ANCHOR_PASSTHROUGH_ERRORS,
        InitialConditionError,
    )
    from proteus.utils.structure_estimate import resolve_P_cmb

    P_cmb = resolve_P_cmb(hf_row, config)[0]
    key = _superliq_cache_key(config, P_cmb)
    if key in _SUPERLIQ_FAILED:
        first = _SUPERLIQ_FAILED[key]
        raise InitialConditionError(str(first)) from first
    try:
        return _solve_superliquidus_adiabat(config, hf_row)
    except InitialConditionError as exc:
        _SUPERLIQ_FAILED[key] = InitialConditionError(str(exc))
        raise
    except (ZalmoxisMissingEOSFilesError, *ANCHOR_PASSTHROUGH_ERRORS):
        raise
    except ANCHOR_NUMERICAL_ERRORS as exc:
        msg = (
            f'liquidus_super: the PALEOS P-T anchor failed at P_cmb={P_cmb / 1e9:.0f} GPa '
            f'({type(exc).__name__}: {exc}).'
        )
        _SUPERLIQ_FAILED[key] = InitialConditionError(msg)
        raise InitialConditionError(msg) from exc

def load_zalmoxis_configuration(
    config: Config,
    hf_row: dict,
    temperature_mode_override: str | None = None,
    external_temperature_source: bool = False,
):
    """Loads the model configuration for Zalmoxis and calculates the dry mass of the planet based on the total mass and the mass of volatiles.
    Args:
        config (Config): The configuration object containing the Zalmoxis parameters.
        hf_row (dict): A dictionary containing the mass of volatiles and other parameters.
        temperature_mode_override: Optional local override for
            ``config.planet.temperature_mode``. Lets callers force a different
            structure-solve mode (e.g. 'adiabatic' for SPIDER coupling with
            T-dependent mantle EOS) without mutating the shared Config object.
            When None, falls back to ``config.planet.temperature_mode``.
        external_temperature_source: True when the caller passes an external
            temperature function or array so the structure solve follows an
            evolved T(r) instead of the internal temperature-mode dispatch. In
            'liquidus_super' mode this lets ``cmb_temperature`` reuse the last
            solved super-liquidus anchor rather than re-solving it, since the
            anchor is not consumed when an external source drives the solve.
    Returns:
        dict: A dictionary containing the Zalmoxis configuration parameters.
    """

    # The Zalmoxis solver consumes core_frac as a mass fraction and does
    # not read core_frac_mode. Warn if a radius fraction was requested, so
    # the user is not surprised by the mass-fraction interpretation.
    if config.interior_struct.core_frac_mode == 'radius':
        log.warning(
            'interior_struct.core_frac_mode = "radius" has no effect with the '
            'zalmoxis module: core_frac (%.3f) is interpreted as a mass '
            'fraction. Use core_frac_mode = "mass", or switch to the dummy or '
            'spider module if a radius fraction is intended.',
            config.interior_struct.core_frac,
        )

    # Setup target planet mass (input parameter) as the total mass of the planet (dry mass + volatiles) [kg],
    # plus the volatile mass impacts and escape moved, which mass_tot (rock anchor) leaves out.
    from proteus.interior_struct.common import volatile_mass_change

    total_planet_mass = config.planet.mass_tot * M_earth + volatile_mass_change(hf_row)

    log.debug(
        'Total target planet mass (dry mass + volatiles): %s kg '
        'with EOS: core=%s, mantle=%s, ice=%s',
        total_planet_mass,
        config.interior_struct.zalmoxis.core_eos,
        config.interior_struct.zalmoxis.mantle_eos,
        config.interior_struct.zalmoxis.ice_layer_eos or 'none',
    )

    # Volatile mass (O included) excluded from the structure target: the full inventory
    # with dry_mantle, else only the atmosphere, since the wet-mantle EOS holds the dissolved
    # mass. Element columns absent before the first inventory count as zero (.get below).
    dry_mantle = config.interior_struct.zalmoxis.dry_mantle
    M_volatiles = 0.0
    for e in element_list:
        if dry_mantle:
            M_volatiles += float(hf_row.get(e + '_kg_total', 0.0))
        else:
            M_volatiles += float(hf_row.get(e + '_kg_atm', 0.0))

    log.debug(f'Volatile mass: {M_volatiles} kg')
    log.debug(
        'Mass budget: total=%.6e kg (%.4f M_earth), volatiles=%.6e kg (%.2f%%)',
        total_planet_mass,
        total_planet_mass / M_earth,
        M_volatiles,
        100.0 * M_volatiles / total_planet_mass if total_planet_mass > 0 else 0,
    )

    # Calculate the target planet mass (dry mass) by subtracting the mass of volatiles from the total planet mass
    planet_mass = total_planet_mass - M_volatiles

    log.debug(f'Target planet mass (dry mass): {planet_mass} kg ')

    # Build per-layer EOS config dict from PROTEUS config fields
    layer_eos_config = {
        'core': config.interior_struct.zalmoxis.core_eos,
        'mantle': config.interior_struct.zalmoxis.mantle_eos,
    }
    if config.interior_struct.zalmoxis.ice_layer_eos is not None:
        layer_eos_config['ice_layer'] = config.interior_struct.zalmoxis.ice_layer_eos

    # This dict only sees the dry layer_eos_config here; zalmoxis_solver
    # rebuilds it after extend_mantle_eos_with_volatiles adds a
    # dissolved-volatile component, so that component gets the real mzf too.
    mzf = config.interior_struct.zalmoxis.mushy_zone_factor
    mushy_zone_factors = _build_mushy_zone_factors(layer_eos_config, mzf)

    zc = config.interior_struct.zalmoxis
    log.debug(
        'Zalmoxis config: num_layers=%d, outer_solver=%s, tol_outer=%.1e, '
        'tol_inner=%.1e, use_jax=%s, use_anderson=%s',
        zc.num_levels,
        zc.outer_solver,
        zc.solver_tol_outer,
        zc.solver_tol_inner,
        zc.use_jax,
        zc.use_anderson,
    )

    # Core fraction. The Zalmoxis solver consumes `core_mass_fraction`
    # literally as a mass fraction and does not read `core_frac_mode`, so
    # `core_frac` is always interpreted as a mass fraction here regardless
    # of the mode. `core_frac_mode = "radius"` is only honoured by the
    # dummy and spider structure modules; a warning is emitted above when
    # it is combined with zalmoxis.

    # 'liquidus_super' can read the mantle EOS table here before
    # zalmoxis_solver's own pre-flight check runs, raising a bare
    # FileNotFoundError. That check never inspects the melting-curve
    # directory a WolfBower2018/RTPress100TPa mantle also needs, so a
    # missing one re-raises unconverted here and in zalmoxis_solver itself.
    try:
        cmb_temperature = _resolve_zalmoxis_cmb_temperature(
            config,
            hf_row,
            temperature_mode_override or config.planet.temperature_mode,
            external_temperature_source=external_temperature_source,
        )
    except FileNotFoundError:
        check_zalmoxis_eos_files(layer_eos_config, load_zalmoxis_material_dictionaries())
        raise

    return {
        'planet_mass': planet_mass,
        'core_mass_fraction': config.interior_struct.core_frac,
        'core_frac_mode': config.interior_struct.core_frac_mode,
        'mantle_mass_fraction': config.interior_struct.zalmoxis.mantle_mass_fraction,
        # For the structure solve, 'accretion' and 'isentropic' both reduce
        # to 'adiabatic' inside Zalmoxis. Zalmoxis only solves the structure
        # (M-R via hydrostatic + EOS); the entropy IC for Aragog/SPIDER is
        # set independently from config.planet.ini_entropy. 'accretion'
        # delays White+Li T-profile until after structure converges;
        # 'isentropic' (CHILI protocol) means the energetics solver
        # consumes ini_entropy, not the Zalmoxis T-profile, so the choice
        # of structure-solve T is decoupled from the IC. 'adiabatic_from_cmb'
        # is passed through to Zalmoxis with the CMB-anchor temperature so
        # the structure-side T(r) integrates upward from T_cmb at R_cmb,
        # matching the entropy that the energetics solver receives via
        # compute_initial_entropy. 'liquidus_super' maps to
        # 'adiabatic_from_cmb' here, with cmb_temperature derived from the
        # Fei+2021 liquidus at the converged P_cmb (or a Noack & Lasbleis
        # (2020) mass-aware P_cmb estimate on the very first call before
        # Zalmoxis has populated P_cmb) plus delta_T_super. The anchor is
        # recomputed against the converged P_cmb on the next structure
        # iteration. temperature_mode_override lets SPIDER coupling force
        # adiabatic without mutating the shared Config object (see proteus
        # rules §"Config mutability").
        'temperature_mode': _resolve_zalmoxis_temperature_mode(
            temperature_mode_override or config.planet.temperature_mode
        ),
        'surface_temperature': config.planet.tsurf_init,
        'cmb_temperature': cmb_temperature,
        'center_temperature': config.planet.tcenter_init,
        'temp_profile_file': None,
        'layer_eos_config': layer_eos_config,
        'mushy_zone_factor': mzf,
        'mushy_zone_factors': mushy_zone_factors,
        'num_layers': config.interior_struct.zalmoxis.num_levels,
        'target_surface_pressure': _get_target_surface_pressure(config, hf_row),
        # Solver tolerances and iteration limits
        'tolerance_outer': config.interior_struct.zalmoxis.solver_tol_outer,
        'tolerance_inner': config.interior_struct.zalmoxis.solver_tol_inner,
        'max_iterations_outer': config.interior_struct.zalmoxis.solver_max_iter_outer,
        'max_iterations_inner': config.interior_struct.zalmoxis.solver_max_iter_inner,
        # JAX diffrax structure path and Anderson Picard acceleration both
        # default on. See Zalmoxis.use_jax and use_anderson in config._struct.
        'use_jax': config.interior_struct.zalmoxis.use_jax,
        'use_anderson': config.interior_struct.zalmoxis.use_anderson,
        # outer mass-radius solver dispatch ('newton' default |
        # 'picard'). When 'newton', Zalmoxis uses Newton + brentq
        # bracketing on f(R) = M(R) - M_target instead of the
        # damped-Picard fixed-point loop. Newton requires tight
        # integrator tolerances; we auto-apply newton_relative_tolerance
        # / newton_absolute_tolerance when the Newton path is selected.
        'outer_solver': config.interior_struct.zalmoxis.outer_solver,
        # Newton-specific knobs (newton_max_iter, newton_tol) AND
        # tightened integrator tolerances are passed ONLY when the
        # Newton path is selected. Picard runs see the dict without
        # any Newton keys and without tightened tolerances, so a
        # future Zalmoxis guard against unknown keys would not break
        # Picard callers.
        **(
            {
                'newton_max_iter': (config.interior_struct.zalmoxis.newton_max_iter),
                'newton_tol': config.interior_struct.zalmoxis.newton_tol,
                'relative_tolerance': (
                    config.interior_struct.zalmoxis.newton_relative_tolerance
                ),
                'absolute_tolerance': (
                    config.interior_struct.zalmoxis.newton_absolute_tolerance
                ),
            }
            if config.interior_struct.zalmoxis.outer_solver == 'newton'
            else {}
        ),
    }

def load_zalmoxis_material_dictionaries():
    """Build an EOS registry dict with file paths pointing to FWL_DATA.

    Returns the same dict format as Zalmoxis ``EOS_REGISTRY``, but with
    every ``eos_file`` path resolved to its fwl-io dataset directory in
    ``FWL_DATA`` instead of ``ZALMOXIS_ROOT/data/``. The data root is the one
    ``download_zalmoxis_eos()`` fetches into, so Zalmoxis, when called from
    PROTEUS, reads the files that were fetched.

    Returns
    -------
    dict
        Flat dict keyed by EOS identifier string (e.g.
        ``"Seager2007:iron"``, ``"PALEOS:MgSiO3"``, ``"Chabrier:H"``).
    """
    root = GetFWLData()
    seager_dir = dataset_dir(EOS_SEAGER_2007, data_root=root)

    _seager_iron = {'eos_file': str(seager_dir / 'eos_seager07_iron.txt')}
    _seager_silicate = {'eos_file': str(seager_dir / 'eos_seager07_silicate.txt')}
    _seager_water = {'eos_file': str(seager_dir / 'eos_seager07_water.txt')}

    # Wolf & Bower 2018
    wb_dir = dataset_dir(EOS_WOLF_BOWER_2018, data_root=root)
    _wb_melted = {
        'eos_file': str(wb_dir / 'density_melt.dat'),
        'adiabat_grad_file': str(wb_dir / 'adiabat_temp_grad_melt.dat'),
    }
    _wb_solid = {'eos_file': str(wb_dir / 'density_solid.dat')}

    # RTPress 100 TPa
    rt_dir = dataset_dir(EOS_RTPRESS_100TPA, data_root=root)
    _rt_melted = {
        'eos_file': str(rt_dir / 'density_melt.dat'),
        'adiabat_grad_file': str(rt_dir / 'adiabat_temp_grad_melt.dat'),
    }

    # PALEOS 2-phase MgSiO3 (separate solid/liquid, Zenodo 19680050): the
    # 150 pts/decade (default) and 600 pts/decade (highres) pairs share one dataset.
    paleos2ph_dir = dataset_dir(EOS_PALEOS_MGSIO3, data_root=root)
    paleos2ph_hr_dir = paleos2ph_dir
    _paleos2ph_melted = {
        'eos_file': str(paleos2ph_dir / 'paleos_mgsio3_tables_pt_proteus_liquid.dat'),
        'format': 'paleos',
    }
    _paleos2ph_solid = {
        'eos_file': str(paleos2ph_dir / 'paleos_mgsio3_tables_pt_proteus_solid.dat'),
        'format': 'paleos',
    }
    _paleos2ph_melted_highres = {
        'eos_file': str(
            paleos2ph_hr_dir / 'paleos_mgsio3_tables_pt_proteus_liquid_highres.dat'
        ),
        'format': 'paleos',
    }
    _paleos2ph_solid_highres = {
        'eos_file': str(paleos2ph_hr_dir / 'paleos_mgsio3_tables_pt_proteus_solid_highres.dat'),
        'format': 'paleos',
    }

    # PALEOS unified tables, one dataset per material
    _paleos_iron = {
        'eos_file': str(
            dataset_dir(EOS_PALEOS_IRON, data_root=root) / 'paleos_iron_eos_table_pt.dat'
        ),
        'format': 'paleos_unified',
    }
    _paleos_mgsio3 = {
        'eos_file': str(
            dataset_dir(EOS_PALEOS_MGSIO3_UNIFIED, data_root=root)
            / 'paleos_mgsio3_eos_table_pt.dat'
        ),
        'format': 'paleos_unified',
    }
    _paleos_h2o = {
        'eos_file': str(
            dataset_dir(EOS_PALEOS_H2O, data_root=root) / 'paleos_water_eos_table_pt.dat'
        ),
        'format': 'paleos_unified',
    }

    # Chabrier H/He
    chabrier_dir = dataset_dir(EOS_CHABRIER_2021, data_root=root)
    _chabrier_h = {
        'eos_file': str(chabrier_dir / 'EOS_Chabrier2021_HHe' / 'chabrier2021_H.dat'),
        'format': 'paleos_unified',
    }

    # PALEOS-API live tabulation entries. These carry only a GridSpec at
    # build time; the dispatch layer (``zalmoxis.eos.paleos_api_cache``)
    # mutates them in place on first density query to populate ``eos_file``
    # and rewrite ``format`` to the downstream value (``paleos_unified`` or
    # ``paleos``). Cache keys are SHA+grid-hash under
    # ``$ZALMOXIS_ROOT/data/EOS_PALEOS_API/``; cold-cache cost is a one-time
    # generator run (see ``zalmoxis.eos.paleos_api``).
    from zalmoxis.eos.paleos_api import (
        make_default_grid_h2o,
        make_default_grid_iron,
        make_default_grid_mgsio3,
    )

    _paleos_api_iron = {
        'format': 'paleos_api',
        'material': 'iron',
        'grid_spec': make_default_grid_iron(),
    }
    _paleos_api_mgsio3 = {
        'format': 'paleos_api',
        'material': 'mgsio3',
        'grid_spec': make_default_grid_mgsio3(),
    }
    _paleos_api_h2o = {
        'format': 'paleos_api',
        'material': 'h2o',
        'grid_spec': make_default_grid_h2o(),
        'h2o_table_path': None,
    }
    _paleos_api_2ph_mgsio3_melted = {
        'format': 'paleos_api_2phase',
        'material': 'mgsio3',
        'side': 'liquid',
        'grid_spec': make_default_grid_mgsio3(),
    }
    _paleos_api_2ph_mgsio3_solid = {
        'format': 'paleos_api_2phase',
        'material': 'mgsio3',
        'side': 'solid',
        'grid_spec': make_default_grid_mgsio3(),
    }

    return {
        # Seager2007 static
        'Seager2007:iron': {'core': _seager_iron},
        'Seager2007:MgSiO3': {'mantle': _seager_silicate},
        'Seager2007:H2O': {'ice_layer': _seager_water},
        # Wolf & Bower 2018 T-dependent
        'WolfBower2018:MgSiO3': {
            'core': _seager_iron,
            'melted_mantle': _wb_melted,
            'solid_mantle': _wb_solid,
        },
        # RTPress 100 TPa extended melt + WB2018 solid
        'RTPress100TPa:MgSiO3': {
            'core': _seager_iron,
            'melted_mantle': _rt_melted,
            'solid_mantle': _wb_solid,
        },
        # PALEOS 2-phase MgSiO3 (Zenodo 19680050; 150 pts/decade default,
        # 600 pts/decade as -highres variant for sensitivity tests).
        'PALEOS-2phase:MgSiO3': {
            'core': _seager_iron,
            'melted_mantle': _paleos2ph_melted,
            'solid_mantle': _paleos2ph_solid,
        },
        'PALEOS-2phase:MgSiO3-highres': {
            'core': _seager_iron,
            'melted_mantle': _paleos2ph_melted_highres,
            'solid_mantle': _paleos2ph_solid_highres,
        },
        # PALEOS unified
        'PALEOS:iron': _paleos_iron,
        'PALEOS:MgSiO3': _paleos_mgsio3,
        'PALEOS:H2O': _paleos_h2o,
        # PALEOS-API live-tabulated (dispatch populates eos_file on demand)
        'PALEOS-API:iron': _paleos_api_iron,
        'PALEOS-API:MgSiO3': _paleos_api_mgsio3,
        'PALEOS-API:H2O': _paleos_api_h2o,
        'PALEOS-API-2phase:MgSiO3': {
            'core': _seager_iron,
            'melted_mantle': _paleos_api_2ph_mgsio3_melted,
            'solid_mantle': _paleos_api_2ph_mgsio3_solid,
        },
        # Chabrier H/He
        'Chabrier:H': _chabrier_h,
    }

def resolve_paleos_api(key: str, entry: dict) -> None:
    """Materialise a PALEOS-API registry entry in place, once per process.

    A resolve runs git to read the paleos version, and the EOS check runs before
    every structure solve on a fresh registry, so resolved entries are kept by key.
    A kept entry whose tables were removed is resolved again, which rebuilds them.

    Parameters
    ----------
    key : str
        Registry key of the entry, e.g. ``'PALEOS-API-2phase:MgSiO3'``.
    entry : dict
        Its registry entry, updated in place with the table paths.

    Raises
    ------
    ZalmoxisMissingEOSFilesError
        When the paleos package is not installed, since the resolver's worker pool
        would then restart forever instead of failing.
    """
    cached = _PALEOS_API_RESOLVED.pop(key, None)
    if cached and all(os.path.isfile(p) for p in _table_paths(cached)):
        _PALEOS_API_RESOLVED[key] = cached
        entry.update(copy.deepcopy(cached))
        return
    if not _paleos_installed():
        raise ZalmoxisMissingEOSFilesError(
            f'PALEOS-API tables not built:\n  {key} (the paleos package is not installed)\n'
            f'{_PALEOS_API_FIX}'
        )
    from zalmoxis.eos.paleos_api_cache import resolve_registry_entry

    resolve_registry_entry(entry)
    _PALEOS_API_RESOLVED[key] = copy.deepcopy(entry)

def check_zalmoxis_eos_files(
    layer_eos_config: dict, mat_dicts: dict, paleos_companions: bool = False
) -> None:
    """Fail fast when a selected EOS table file is missing on disk.

    Walks the registry entries selected by ``layer_eos_config`` and
    collects every referenced table path that does not exist, then
    raises one actionable error. Without this check the solver emits
    one read error per shell and ends in a non-convergence failure that
    hides the real cause. PALEOS-API entries (live tabulation) are
    materialised here, so a failed build stops the run before any solve.

    Parameters
    ----------
    layer_eos_config : dict
        Per-layer EOS identifier strings (``'core'``, ``'mantle'``, ...).
    mat_dicts : dict
        EOS registry from :func:`load_zalmoxis_material_dictionaries`.
    paleos_companions : bool
        Also require the :func:`paleos_companion_keys` of the mantle.

    Raises
    ------
    ZalmoxisMissingEOSFilesError
        If any selected EOS table file is missing, naming every missing
        path and the command that downloads them.
    """
    from zalmoxis.eos.dispatch import _is_paleos_api

    selected = [
        (role, key)
        for role, identifier in layer_eos_config.items()
        for key in eos_components(identifier)
    ]
    if paleos_companions:
        companions = paleos_companion_keys(layer_eos_config.get('mantle', ''))
        selected += [('mantle', key) for key in companions]
    missing: set[str] = set()
    unbuilt: set[str] = set()
    for layer_role, key in selected:
        entry = mat_dicts.get(key)
        if entry is None:
            # Unknown identifiers fail later with a registry error.
            continue
        api = key in _PALEOS_API_RESOLVED or _is_paleos_api(entry)
        if api:
            try:
                resolve_paleos_api(key, entry)
            except ZalmoxisMissingEOSFilesError:
                unbuilt.add(f'{key} (the paleos package is not installed)')
                continue
            except Exception as exc:
                unbuilt.add(f'{key} ({type(exc).__name__}: {exc})')
                continue
        # Nested entries map layer roles to flat entries; their 'core' sub-entry
        # counts only for the core role or when it is the only sub-entry.
        if 'eos_file' in entry:
            subentries = [entry]
        else:
            has_other_roles = any(role != 'core' for role in entry)
            subentries = [
                sub
                for role, sub in entry.items()
                if role != 'core' or layer_role == 'core' or not has_other_roles
            ]
        for sub in subentries:
            if not isinstance(sub, dict):
                continue
            for field in ('eos_file', 'adiabat_grad_file'):
                path = sub.get(field)
                if path and not os.path.isfile(path) and api:
                    # A generated table removed since it was built is built again.
                    _PALEOS_API_RESOLVED.pop(key, None)
                    unbuilt.add(f'{key} ({path} not found)')
                elif path and not os.path.isfile(path):
                    missing.add(path)
    parts = []
    if missing:
        parts.append(
            f'Interior EOS table file(s) not found:\n  {"\n  ".join(sorted(missing))}\n'
            'Download them with '
            '`proteus get interiordata --config-path <config.toml>`, '
            f'or run `proteus start` once without --offline. {RELOCATE_HINT}'
        )
    if unbuilt:
        parts.append(
            f'PALEOS-API tables not built:\n  {"\n  ".join(sorted(unbuilt))}\n{_PALEOS_API_FIX}'
        )
    if parts:
        raise ZalmoxisMissingEOSFilesError('\n'.join(parts))

def energetics_entry(mantle_eos: str, mat_dicts: dict) -> tuple[str | None, dict]:
    """Energetics key of a mantle EOS and its registry entry (empty when absent)."""
    key = energetics_eos_key(mantle_eos)
    return key, mat_dicts.get(key or '', {})

def require_paleos_tables(config: Config, outdir: str) -> None:
    """Stop the run before any solve when a table of its EOS set is missing.

    Requires every table file of the core, mantle and ice-layer EOS and,
    for a mantle with a PALEOS component, the MgSiO3 2-phase pair; a resumed
    SPIDER run that keeps its P-S tables is held only to its layer tables. The
    liquidus_super initial adiabat requires the pair for any mantle. With
    ``dry_mantle = false`` the dissolved-volatile tables are also required. A
    mantle with a non-MgSiO3 component gets one WARNING that names the tables
    its energetics and its liquidus_super initial adiabat use.

    Parameters
    ----------
    config : Config
        Configuration object with struct.zalmoxis settings.
    outdir : str
        The run output directory.

    Raises
    ------
    ZalmoxisMissingEOSFilesError
        Naming every missing file and the command that downloads it.
    """
    zc = config.interior_struct.zalmoxis
    layers = {'core': zc.core_eos, 'mantle': zc.mantle_eos}
    if zc.ice_layer_eos is not None:
        layers['ice_layer'] = zc.ice_layer_eos
    if not zc.dry_mantle:
        # Dissolved volatiles join the mantle EOS during the run.
        layers['volatiles'] = '+'.join(VOLATILE_EOS_MAP.values())
    mat_dicts = load_zalmoxis_material_dictionaries()
    # Runs before the accreted mass is restored: no mass check and no report here.
    kept = config.params.resume and _resumed_ps_tables(outdir, None)
    liquidus_super = config.planet.temperature_mode == 'liquidus_super'
    if liquidus_super:
        layers['anchor'] = twophase_registry_key(zc.mantle_eos)
    # Aragog re-solves the initial condition on the pair when the mesh changes.
    check_zalmoxis_eos_files(
        layers,
        mat_dicts,
        paleos_companions=not kept or config.interior_energetics.module == 'aragog',
    )
    components = eos_components(zc.mantle_eos)
    mixture = len(set(components)) > 1
    if all(is_mgsio3(c) for c in components):
        return
    uses = []
    if config.interior_energetics.module in ('spider', 'aragog'):
        if not generates_paleos_tables(config.interior_struct):
            uses.append(
                'the energetics use the eos_dir P-S tables (by default Wolf and Bower 2018) '
                'and the melting_dir curves'
            )
        else:
            uses.append(
                'the energetics and melting curves are PALEOS MgSiO3 '
                f'({energetics_eos_key(zc.mantle_eos)}, solidus = '
                f'{zc.mushy_zone_factor:.2f} x PALEOS liquidus)'
            )
    if liquidus_super:
        uses.append('the liquidus_super initial adiabat is solved on the MgSiO3 2-phase tables')
    if uses:
        head = (
            'the non-MgSiO3 components enter only the structure density; '
            if mixture
            else 'the structure uses its density, while '
        )
        log.warning('mantle_eos=%s: %s%s', zc.mantle_eos, head, '; '.join(uses))

def resolve_2phase_mgsio3_paths(mantle_eos: str, mat_dicts: dict, required: bool = False):
    """Return the 2-phase MgSiO3 table paths, stopping on a missing one when required.

    Parameters
    ----------
    mantle_eos : str
        Configured mantle EOS.
    mat_dicts : dict
        Zalmoxis material dictionaries.
    required : bool
        Raise instead of returning None for a missing table. The initial-condition
        solves that read the pair always set it, the surface-entropy lookup sets it
        for a generated PALEOS table set, and the P-S export leaves it off.

    Returns
    -------
    tuple[str | None, str | None]
        See :func:`_twophase_mgsio3_paths`.

    Raises
    ------
    ZalmoxisMissingEOSFilesError
        If ``required`` and a table of the pair is not available.
    """
    paths = _twophase_mgsio3_paths(mantle_eos, mat_dicts)
    if required and None in paths:
        missing = [p for p, ok in zip(('solid', 'liquid'), paths) if ok is None]
        key = twophase_registry_key(mantle_eos)
        fix = (
            _PALEOS_API_FIX
            if key.startswith('PALEOS-API')
            else 'Download them with `proteus get interiordata --config-path <config.toml>`. '
            + RELOCATE_HINT
        )
        raise ZalmoxisMissingEOSFilesError(
            f'PALEOS 2-phase MgSiO3 tables {key} not available: {", ".join(missing)}. {fix}'
        )
    return paths

def load_zalmoxis_solidus_liquidus_functions(mantle_eos: str, config: Config):
    """Loads the solidus and liquidus functions for Zalmoxis based on the mantle EOS.

    Melting curves are needed for two purposes:
    1. Temperature-dependent density in the mushy zone (WolfBower2018, RTPress).
    2. phi(r) blending in VolatileProfile (any EOS with dissolved volatiles).

    The curves follow the energetics key of the mantle EOS
    (:func:`proteus.utils.helper.energetics_eos_key`), the same key as the
    energetics. For WolfBower2018/RTPress100TPa, loads SPIDER-format P-T files from FWL_DATA.
    For PALEOS unified and PALEOS-2phase, the liquidus comes from the analytic
    Belonoshko+2005 / Fei+2021 curve (Zalmoxis ``'PALEOS-liquidus'``) which is
    the basis Zalmoxis uses for MgSiO3 phase separation, and the solidus is
    derived as T_sol = T_liq * mushy_zone_factor. This keeps the curves used
    for phi-blending and 2-phase nabla_ad consistent with the unified PALEOS
    density-interpolation phase boundaries.

    Parameters
    ----------
    mantle_eos : str
        Mantle EOS string (e.g. ``"WolfBower2018:MgSiO3"``, ``"PALEOS:MgSiO3"``,
        ``"PALEOS-2phase:MgSiO3"``).
    config : Config
        PROTEUS configuration object.

    Returns
    -------
    tuple or None
        (solidus_func, liquidus_func) callable P [Pa] -> T [K], or None.
    """
    key = energetics_eos_key(mantle_eos) or ''
    if key.startswith(TDEP_EOS_PREFIXES):
        return get_zalmoxis_melting_curves(config)

    # PALEOS unified and 2-phase share the analytic PALEOS-liquidus (Belonoshko+2005 /
    # Fei+2021), with solidus = mushy_zone_factor x liquidus. The 2-phase nabla_ad needs
    # both curves, and without them the unified path falls back to phi = 0.5.
    if key.startswith(PALEOS_EOS_PREFIXES):
        try:
            from zalmoxis.melting_curves import (
                derive_solidus_from_liquidus,
                get_solidus_liquidus_functions,
            )

            _, liquidus_func = get_solidus_liquidus_functions(
                solidus_id='Stixrude14-solidus',  # required by API but unused; solidus is built below as mushy_zone_factor * liquidus (constant depression of the PALEOS liquidus)
                liquidus_id='PALEOS-liquidus',
            )
            mzf = config.interior_struct.zalmoxis.mushy_zone_factor
            solidus_func = derive_solidus_from_liquidus(liquidus_func, mzf)
            log.info(
                'PALEOS melting curves (%s): liquidus from PALEOS, '
                'solidus = liquidus * %.2f (mushy_zone_factor)',
                mantle_eos,
                mzf,
            )
            return solidus_func, liquidus_func
        except Exception as e:
            log.warning('Could not load PALEOS melting curves: %s', e)
            return None

    return None

def scale_temperature_profile_for_aragog(
    config: Config, mantle_radii: np.ndarray, mantle_temperature_profile: np.ndarray
):
    """Scales the temperature profile obtained from Zalmoxis to match the number of levels required by Aragog.
    Args:
        config (Config): The configuration object containing the configuration parameters.
        mantle_radii (np.ndarray): The radial positions of the mantle layers from Zalmoxis.
        mantle_temperature_profile (np.ndarray): The temperature profile of the mantle layers from Zalmoxis.
    Returns:
        np.ndarray: The scaled temperature profile matching the number of levels in Aragog.
    """

    # Number of levels in Aragog mesh
    mesh_grid_size = config.interior_energetics.num_levels - 1

    # Create new evenly spaced radial positions for Aragog
    radii_to_interpolate = np.linspace(mantle_radii[0], mantle_radii[-1], mesh_grid_size)

    # Cubic interpolation onto the Aragog radial mesh
    cubic_interp = interp1d(mantle_radii, mantle_temperature_profile, kind='cubic')
    return cubic_interp(radii_to_interpolate)

def write_spider_mesh_file(
    outdir: str,
    mantle_radii: np.ndarray,
    mantle_pressure: np.ndarray,
    mantle_density: np.ndarray,
    mantle_gravity: np.ndarray,
    num_basic: int,
) -> str:
    """Write an external mesh file for SPIDER from Zalmoxis mantle profiles.

    Interpolates the Zalmoxis mantle arrays onto uniformly-spaced SPIDER
    basic and staggered nodes, then writes the mesh file in the format
    expected by SPIDER's ``SetMeshFromExternalFile()``.

    Parameters
    ----------
    outdir : str
        PROTEUS output directory (file is written to ``outdir/data/``).
    mantle_radii : np.ndarray
        Radial positions from CMB to surface, ascending [m].
    mantle_pressure : np.ndarray
        Pressure at each radius [Pa].
    mantle_density : np.ndarray
        Density at each radius [kg/m^3].
    mantle_gravity : np.ndarray
        Gravity magnitude at each radius [m/s^2] (positive).
    num_basic : int
        Number of SPIDER basic nodes (shell boundaries).

    Returns
    -------
    str
        Path to the written mesh file.
    """
    num_staggered = num_basic - 1
    R_surf = float(mantle_radii[-1])
    R_cmb = float(mantle_radii[0])

    # Basic nodes: uniform spacing from surface to CMB (descending r)
    r_b = np.linspace(R_surf, R_cmb, num_basic)
    # Staggered nodes: midpoints between consecutive basic nodes
    r_s = 0.5 * (r_b[:-1] + r_b[1:])

    # Interpolate Zalmoxis profiles onto node positions
    # mantle_radii is ascending, np.interp requires ascending xp
    P_b = np.interp(r_b, mantle_radii, mantle_pressure)
    rho_b = np.interp(r_b, mantle_radii, mantle_density)
    g_b = np.interp(r_b, mantle_radii, mantle_gravity)

    P_s = np.interp(r_s, mantle_radii, mantle_pressure)
    rho_s = np.interp(r_s, mantle_radii, mantle_density)
    g_s = np.interp(r_s, mantle_radii, mantle_gravity)

    # Negate gravity for SPIDER convention (inward-pointing, negative)
    g_b = -np.abs(g_b)
    g_s = -np.abs(g_s)

    # Write mesh file
    mesh_path = os.path.join(outdir, 'data', 'spider_mesh.dat')
    with open(mesh_path, 'w') as f:
        f.write(f'# {num_basic} {num_staggered}\n')
        for i in range(num_basic):
            f.write(f'{r_b[i]:.15e} {P_b[i]:.15e} {rho_b[i]:.15e} {g_b[i]:.15e}\n')
        for i in range(num_staggered):
            f.write(f'{r_s[i]:.15e} {P_s[i]:.15e} {rho_s[i]:.15e} {g_s[i]:.15e}\n')

    log.info(
        'Wrote SPIDER mesh file: %s (%d basic + %d staggered nodes)',
        mesh_path,
        num_basic,
        num_staggered,
    )

    return mesh_path

def read_ps_cache_pointer(outdir: str) -> str | None:
    """Return the shared PS-cache table directory recorded for `outdir`.

    Reads the pointer written by :func:`_write_ps_cache_pointer`. Used by
    the resume path to locate PROTEUS_PS_CACHE_DIR tables that live outside
    the run directory.

    Parameters
    ----------
    outdir : str
        The run output directory.

    Returns
    -------
    str or None
        The recorded table directory, or None when the pointer is absent,
        unreadable, or empty.
    """

    pointer = os.path.join(outdir, 'data', PS_CACHE_POINTER_NAME)
    if not os.path.isfile(pointer):
        return None
    try:
        with open(pointer) as f:
            cache_dir = f.read().strip()
    except OSError:
        return None
    return cache_dir or None

def generate_spider_tables(config: Config, outdir: str):
    """Generate P-S EOS tables and phase boundaries from PALEOS data.

    Produces P-S lookup tables for density, temperature, heat capacity,
    thermal expansion, and adiabatic gradient, plus solidus/liquidus phase
    boundaries in S(P) format. These are consumed by the entropy-IC verify
    in Aragog (and by SPIDER if the structure module is SPIDER).

    Supports two PALEOS layouts:

    1. ``paleos_unified`` (e.g. ``PALEOS:MgSiO3``): the structural backbone is
       the single unified P-T table covering both phases plus mushy zone, while
       the per-phase property surfaces are built from the sibling two-phase
       solid + liquid tables, so the densities stay resolved across the melting-curve
       discontinuity. The solidus is derived from ``mushy_zone_factor *
       liquidus`` (default 0.8, the constant Stixrude 2014 solidus/liquidus
       ratio applied to the PALEOS liquidus); the liquidus is the analytic
       PALEOS Belonoshko+2005 / Fei+2021 curve.
    2. ``PALEOS-2phase:<solid>`` (e.g. ``PALEOS-2phase:MgSiO3``): separate
       solid + liquid PALEOS tables. The solidus is derived exactly as in the
       unified layout, ``mushy_zone_factor * liquidus`` (default 0.8); the
       two-phase tables additionally supply the latent-heat entropy gap
       between the liquid-table entropy at the liquidus and the solid-table
       entropy at the derived solidus.

    For non-PALEOS EOS types (WolfBower2018, RTPress100TPa), returns None
    and the caller is expected to fall back on pre-existing SPIDER tables.
    A PALEOS mantle with a missing table raises instead.

    Parameters
    ----------
    config : Config
        Configuration object with struct.zalmoxis settings.
    outdir : str
        Output directory. Tables are written to ``outdir/data/spider_eos/``, or,
        when the ``PROTEUS_PS_CACHE_DIR`` environment variable is set, to a
        subdirectory of it named after the sanitised :func:`_ps_cache_key` string, which
        independent runs with the same key share.

    Returns
    -------
    dict or None
        Keys ``'eos_dir'``, ``'solidus_path'``, ``'liquidus_path'`` with
        absolute paths. Returns None if the mantle EOS is not PALEOS.

    Raises
    ------
    ZalmoxisMissingEOSFilesError
        When a table of a PALEOS mantle or its MgSiO3 2-phase pair is missing.
    """
    from zalmoxis.eos_export import generate_spider_eos_tables, generate_spider_phase_boundaries
    from zalmoxis.melting_curves import (
        derive_solidus_from_liquidus,
        get_solidus_liquidus_functions,
    )

    mantle_eos = config.interior_struct.zalmoxis.mantle_eos

    # Use FWL_DATA paths (not ZALMOXIS_ROOT) for EOS file lookup
    mat_dicts = load_zalmoxis_material_dictionaries()
    key, eos_entry = energetics_entry(mantle_eos, mat_dicts)

    # A resumed run stays on the tables it started with; the key only feeds the warning.
    if config.params.resume:
        dropped = []
        resumed = _resumed_ps_tables(
            outdir,
            lambda: _ps_resume_key(config, key, eos_entry, mat_dicts),
            _ps_p_max(config),
            dropped,
        )
        if resumed is not None:
            return resumed

    if not eos_entry or key not in PALEOS_REGISTRY_KEYS:
        log.info(
            'Mantle EOS %s has no PALEOS table set; using pre-existing SPIDER tables.',
            mantle_eos,
        )
        return None
    if paleos_companion_keys(mantle_eos):
        check_zalmoxis_eos_files({'mantle': mantle_eos}, mat_dicts, paleos_companions=True)
    try:
        inputs = _ps_table_inputs(config, key, eos_entry, mat_dicts)
    except _NoPSTables as exc:
        log.log(exc.level, 'No PALEOS P-S tables: %s; using pre-existing SPIDER tables.', exc)
        return None
    eos_file, solid_eos, liquid_eos, is_twophase, P_max, cache_key = inputs
    if solid_eos and liquid_eos:
        log.info('Using PALEOS-2phase tables for entropy-IC table generation')
    nP = config.interior_struct.zalmoxis.lookup_nP
    nS = config.interior_struct.zalmoxis.lookup_nS
    # Tables dropped for a mass change were reported by _resumed_ps_tables.
    if config.params.resume and not dropped:
        log.warning(
            'Resumed run has no kept P-S entropy tables in %s or at its shared-cache '
            'pointer; it continues on the tables of the current key %s, built now if absent',
            os.path.join(outdir, 'data', 'spider_eos'),
            cache_key,
        )

    mzf = config.interior_struct.zalmoxis.mushy_zone_factor
    # Phase boundaries: PALEOS-liquidus is the analytic Belonoshko+2005 /
    # Fei+2021 Simon-Glatzel curve. The derived solidus is
    # T_solidus = T_liquidus * mushy_zone_factor for both layouts. A two-phase
    # table adds the latent-heat gap between liquid_table.s(P, T_liq) and
    # solid_table.s(P, T_sol) with T_sol = mushy_zone_factor * T_liq.
    _, liquidus_func = get_solidus_liquidus_functions(
        solidus_id='Stixrude14-solidus',  # unused, but API requires it
        liquidus_id='PALEOS-liquidus',
    )
    solidus_func = derive_solidus_from_liquidus(liquidus_func, mzf)
    if is_twophase:
        # This solidus_func also reaches Zalmoxis's own 2-phase structure
        # solve via load_zalmoxis_solidus_liquidus_functions, so mzf moves
        # nabla_ad there too, separately from the density blend below.
        log.info(
            'PALEOS-2phase phase boundaries: solidus = liquidus * %.2f '
            '(mushy_zone_factor); latent heat from 2-phase tables',
            mzf,
        )
    else:
        log.info(
            'PALEOS unified phase boundaries: solidus = liquidus * %.2f (mushy_zone_factor)',
            mzf,
        )

    # Table location: output/<run>/data/spider_eos, or with PROTEUS_PS_CACHE_DIR a
    # shared directory keyed by cache_key (fields in _ps_cache_key), so runs with
    # the same key reuse one slow full-resolution PALEOS P-S table.
    _ps_cache_root = os.environ.get('PROTEUS_PS_CACHE_DIR')
    if _ps_cache_root:
        _safe_key = cache_key.replace('.', 'p').replace('=', '-').replace('+', '')
        spider_eos_dir = os.path.join(_ps_cache_root, _safe_key)
        os.makedirs(spider_eos_dir, exist_ok=True)
        # The shared-cache tables live outside the run directory, so a resumed
        # run cannot find them via the per-run output/<run>/data/spider_eos
        # path. Leave a pointer so resume can follow the tables to the cache.
        _write_ps_cache_pointer(outdir, spider_eos_dir)
    else:
        spider_eos_dir = os.path.join(outdir, 'data', 'spider_eos')

    # Cache check: skip regeneration if tables exist and pressure range unchanged.
    # The pressure range depends on planet mass, which doesn't change during evolution.
    cache_marker = os.path.join(spider_eos_dir, '.cache_info.txt')
    if os.path.isfile(cache_marker):
        with open(cache_marker) as f:
            existing_key = f.read().strip()
        if existing_key == cache_key:
            # Tables are up to date. File names must match the writer in
            # zalmoxis.eos_export.generate_spider_phase_boundaries, which
            # emits solidus_P-S.dat / liquidus_P-S.dat.
            solidus_path = os.path.join(spider_eos_dir, 'solidus_P-S.dat')
            liquidus_path = os.path.join(spider_eos_dir, 'liquidus_P-S.dat')
            if os.path.isfile(solidus_path) and os.path.isfile(liquidus_path):
                log.info(
                    'Reusing cached PALEOS-derived P-S entropy tables (P_max=%.2e, %dx%d)',
                    P_max,
                    nP,
                    nS,
                )
                return {
                    'eos_dir': spider_eos_dir,
                    'solidus_path': solidus_path,
                    'liquidus_path': liquidus_path,
                }
        else:
            log.info(
                'Regenerating P-S entropy tables in %s: cache key %s does not match %s',
                spider_eos_dir,
                existing_key,
                cache_key,
            )

    # Choose where to generate. For a shared PROTEUS_PS_CACHE_DIR the tables are
    # written into a private staging directory on the same filesystem and then
    # published into spider_eos_dir with per-file atomic renames, so two cluster
    # runs that miss the cache marker at the same moment cannot interleave
    # partial writes into the shared table directory. For a per-run directory
    # there is no sharing, so tables are generated in place as before.
    if _ps_cache_root:
        gen_dir = tempfile.mkdtemp(prefix='.gen-', dir=os.path.dirname(spider_eos_dir))
    else:
        gen_dir = spider_eos_dir

    try:
        # Generate phase boundaries
        log.info(
            'Generating PALEOS-derived P-S phase boundaries (%d P points)...',
            nP,
        )
        generate_spider_phase_boundaries(
            solidus_func=solidus_func,
            liquidus_func=liquidus_func,
            eos_file=eos_file,
            P_range=(1e5, P_max),
            n_P=nP,
            output_dir=gen_dir,
            solid_eos_file=solid_eos,
            liquid_eos_file=liquid_eos,
        )

        # Generate full EOS tables
        log.info(
            'Generating PALEOS-derived P-S EOS tables (%d x %d)...',
            nP,
            nS,
        )
        generate_spider_eos_tables(
            eos_file=eos_file,
            solidus_func=solidus_func,
            liquidus_func=liquidus_func,
            P_range=(1e5, P_max),
            n_P=nP,
            n_S=nS,
            output_dir=gen_dir,
            solid_eos_file=solid_eos,
            liquid_eos_file=liquid_eos,
        )

        # Publish staged tables into the shared cache with atomic renames.
        if gen_dir != spider_eos_dir:
            _publish_ps_tables(gen_dir, spider_eos_dir)
    finally:
        if gen_dir != spider_eos_dir:
            shutil.rmtree(gen_dir, ignore_errors=True)

    # Write the cache marker last, atomically, so a concurrent reader only
    # trusts the directory once every table file is already in place.
    try:
        _atomic_write_text(cache_marker, cache_key)
    except OSError:
        pass

    # File names match the writer in zalmoxis.eos_export, which emits
    # solidus_P-S.dat / liquidus_P-S.dat; use the published locations.
    return {
        'eos_dir': spider_eos_dir,
        'solidus_path': os.path.join(spider_eos_dir, 'solidus_P-S.dat'),
        'liquidus_path': os.path.join(spider_eos_dir, 'liquidus_P-S.dat'),
    }