aragog.eos

def dTdPs(self, P: npt.NDArray | float, S: npt.NDArray | float) -> npt.NDArray:
    """Adiabatic temperature gradient dT/dP|_S (P, S) [K/Pa].

    This is the SPIDER convention: dT/dP along the adiabat, NOT
    the dimensionless nabla_ad = d ln T / d ln P.
    """
    return self._lookup_phase_weighted('dTdPs', P, S)

def density(self, P: npt.NDArray | float, S: npt.NDArray | float) -> npt.NDArray:
    """Density rho(P, S) [kg/m^3].

    In the mushy zone (0 < phi < 1): harmonic mean of end-member
    densities evaluated at phase boundary entropies (SPIDER
    eos_composite.c:236-237).

    Outside the mushy zone (phi=0 or phi=1): single-phase table
    evaluated at the actual S (clamped to table range), matching
    SPIDER's combine_matprop(smth=0, mixed, single) path.
    """
    P = np.asarray(P, dtype=float)
    S = np.asarray(S, dtype=float)
    phi = self._melt_fraction(P, S)
    mushy = (phi > 0) & (phi < 1)

    solid_table = self._tables['density_solid']
    melt_table = self._tables['density_melt']

    # Mushy zone: evaluate at phase boundaries (Lever Rule)
    rho_sol_boundary = self._lookup_at_phase_boundary('density', P, 'solid')
    rho_liq_boundary = self._lookup_at_phase_boundary('density', P, 'melt')
    inv_rho_mushy = phi / np.maximum(rho_liq_boundary, 1.0) + (1.0 - phi) / np.maximum(
        rho_sol_boundary, 1.0
    )
    rho_mushy = 1.0 / np.maximum(inv_rho_mushy, 1e-30)

    # Single-phase: evaluate at actual S (clamped to table range),
    # selected below by phi >= 0.5. Mushy points use rho_mushy, and
    # each non-mushy point uses only the branch its selector picks,
    # so mask each range check to the points that use that branch.
    melt_selected = phi >= 0.5
    solid_used = ~mushy & ~melt_selected
    melt_used = ~mushy & melt_selected
    self._check_entropy_range(
        np.where(solid_used, S, solid_table['S'][0]),
        solid_table['S'][0],
        solid_table['S'][-1],
        'density (solid table lookup)',
    )
    self._check_entropy_range(
        np.where(melt_used, S, melt_table['S'][0]),
        melt_table['S'][0],
        melt_table['S'][-1],
        'density (melt table lookup)',
    )
    S_solid_clamped = np.clip(S, solid_table['S'][0], solid_table['S'][-1])
    S_melt_clamped = np.clip(S, melt_table['S'][0], melt_table['S'][-1])
    P_solid_clamped = np.clip(P, solid_table['P'][0], solid_table['P'][-1])
    P_melt_clamped = np.clip(P, melt_table['P'][0], melt_table['P'][-1])
    pts_solid = np.column_stack([P_solid_clamped.ravel(), S_solid_clamped.ravel()])
    pts_melt = np.column_stack([P_melt_clamped.ravel(), S_melt_clamped.ravel()])
    rho_solid_single = solid_table['interp'](pts_solid).reshape(P.shape)
    rho_melt_single = melt_table['interp'](pts_melt).reshape(P.shape)
    # 0.5 here is a binary table-selector for the non-mushy fallback
    # (phi is essentially 0 or 1 outside the mushy band); it is NOT
    # the rheological critical melt fraction. The RCMF is
    # ``EntropyPhaseEvaluator._phi_rheo`` and drives the viscosity
    # tanh blend separately.
    rho_single = np.where(phi >= 0.5, rho_melt_single, rho_solid_single)

    return np.where(mushy, rho_mushy, rho_single)

def entropy_at_temperature(self, P: float, T: float) -> float:
    """Entropy [J/kg/K] at pressure ``P`` where ``temperature(P, S) = T``.

    Brent's method over the tables' entropy range, which keeps its bracket where ``T(S)`` is
    flat; the range edge where ``T`` lies beyond it, and NaN for a NaN ``P`` or ``T`` or where
    the method does not converge.
    """
    if np.isnan(P + T):
        return np.nan
    solid, melt = self._tables['temperature_solid'], self._tables['temperature_melt']
    lo, hi = min(solid['S'][0], melt['S'][0]), max(solid['S'][-1], melt['S'][-1])

    def f(s):
        return self.temperature_scalar(P, s) - T

    f_lo, f_hi = f(lo), f(hi)
    if f_lo >= 0.0 or f_hi <= 0.0:
        return lo if f_lo >= 0.0 else hi
    # next to a flat T(S) Brent takes more than its default 100 iterations
    s, result = brentq(f, lo, hi, maxiter=500, full_output=True, disp=False)
    return s if result.converged else np.nan

def heat_capacity(self, P: npt.NDArray | float, S: npt.NDArray | float) -> npt.NDArray:
    """Specific heat capacity Cp(P, S) [J/kg/K].

    Linear blend of solid and melt P-S tables by melt fraction
    (the ``cp_blend = 'linear'`` convention). Use
    ``heat_capacity_latent_blend`` for SPIDER-parity mushy-zone
    Cp that includes the latent-heat contribution.
    """
    return self._lookup_phase_weighted('heat_capacity', P, S)

def heat_capacity_latent_blend(
    self, P: npt.NDArray | float, S: npt.NDArray | float, width: float = 0.01
) -> npt.NDArray:
    """Latent-heat-augmented Cp(P, S), matching SPIDER convention.

    Mirrors SPIDER's ``eos_composite.c:226-232`` formula:

        Cp_mix = (S_liq - S_sol) / (T_liq - T_sol)
               * (T_sol + 0.5 * (T_liq - T_sol))

    which captures the latent-heat contribution to apparent heat
    capacity in the mushy zone. Outside the phase boundary the
    result reduces to the pure-phase Cp via tanh smoothing.

    The plain linear blend used by ``heat_capacity`` above misses
    the latent-heat term and underestimates Cp by up to a factor
    of 6 in the deep mushy regime relative to SPIDER's internal
    cp_s. This formula closes that gap.

    Parameters
    ----------
    P : array or float
        Pressure [Pa].
    S : array or float
        Entropy [J/kg/K].
    width : float
        Tanh smoothing width across the phase boundary, in units
        of melt fraction. SPIDER default is 0.01.

    Returns
    -------
    ndarray
        Heat capacity [J/kg/K] with the same shape as the inputs.
    """
    P = np.asarray(P, dtype=float)
    S = np.asarray(S, dtype=float)

    # Phase boundaries at this pressure
    S_sol = np.asarray(self.solidus_entropy(P)).reshape(P.shape)
    S_liq = np.asarray(self.liquidus_entropy(P)).reshape(P.shape)
    T_sol = np.asarray(self._lookup_at_phase_boundary('temperature', P, 'solid')).reshape(
        P.shape
    )
    T_liq = np.asarray(self._lookup_at_phase_boundary('temperature', P, 'melt')).reshape(
        P.shape
    )

    dS = S_liq - S_sol
    dT = T_liq - T_sol

    # Mixed-phase Cp from the latent-heat formula
    # SPIDER: Cp_mix = (S_liq - S_sol) / (T_liq - T_sol) * T_mid
    safe_dT = np.where(np.abs(dT) > 1e-10, dT, 1e-10)
    T_mid = T_sol + 0.5 * dT
    Cp_mix = (dS / safe_dT) * T_mid

    # Pure-phase Cp from the linear-blend tables
    Cp_pure = self._lookup_phase_weighted('heat_capacity', P, S)

    # Un-truncated phase fraction (SPIDER's gphi convention)
    safe_dS = np.where(np.abs(dS) > 1e-10, dS, 1e-10)
    gphi = (S - S_sol) / safe_dS

    # SPIDER-parity smoothing (util.c:245-270, get_smoothing).
    # Bell shape centred on the mushy zone: rises at gphi=0 with
    # tanh(gphi/w), falls at gphi=1 with tanh((1-gphi)/w). The
    # width parameter w controls the shoulder into the single-
    # phase regime.
    #
    # IMPORTANT: SPIDER uses tanh(x/w), not expit(x/w) which
    # equals 0.5*(1+tanh(x/(2w))). The previous expit formula
    # had an effective width of w/2, making the transition twice
    # as narrow and producing 16% less latent-heat Cp at
    # crystallization onset (gphi=0.994).
    w = max(width, 1e-6)
    smooth = np.where(
        gphi > 0.5,
        0.5 * (1.0 + np.tanh((1.0 - gphi) / w)),  # fall at liquidus
        0.5 * (1.0 + np.tanh(gphi / w)),  # rise at solidus
    )

    return smooth * Cp_mix + (1.0 - smooth) * Cp_pure

def invert_temperature(self, P: float, T_target: float) -> float:
    """Find entropy S such that T(P, S) = T_target.

    Uses Brent root-finding with pure-Python bilinear interpolation
    on the P-S temperature tables. All arithmetic inside the brentq
    callback uses only Python floats, avoiding numpy scalar
    conversions that break on numpy >= 2.4 (where ndim > 0 arrays
    cannot be implicitly converted to scalars).

    Parameters
    ----------
    P : float
        Pressure [Pa].
    T_target : float
        Target temperature [K].

    Returns
    -------
    float
        Entropy [J/kg/K] such that T(P, S) ~ T_target.

    Raises
    ------
    ValueError
        If T_target is outside the range of T(P, S) for this P.
    """
    from scipy.optimize import brentq

    # All values must be pure Python floats for numpy 2.4 safety.
    P_clamped = max(self.P_min, min(float(P), self.P_max))
    T_tgt = float(T_target)

    def residual(S_cand: float) -> float:
        # Pure-Python path: no numpy arrays created or converted.
        return self.temperature_scalar(P_clamped, S_cand) - T_tgt

    S_lo = self.S_min  # already Python float from __init__
    S_hi = self.S_max
    f_lo = residual(S_lo)
    f_hi = residual(S_hi)

    if f_lo * f_hi > 0:
        raise ValueError(
            f'Cannot invert T={T_target:.1f} K at P={P:.2e} Pa: '
            f'T(S_min={S_lo:.0f})={T_tgt + f_lo:.0f} K, '
            f'T(S_max={S_hi:.0f})={T_tgt + f_hi:.0f} K. '
            f'Target outside table range.'
        )

    S_root = brentq(residual, S_lo, S_hi, xtol=0.1, rtol=1e-10)
    # brentq returns a Python float (from C double), but be safe
    return float(S_root)

def latent_heat(self, P: npt.NDArray | float) -> npt.NDArray:
    """Latent heat L(P) = T_fus × (S_liq - S_sol) [J/kg].

    P-dependent, following SPIDER convention. T_fus is the average
    of solidus and liquidus temperatures at the given pressure.
    """
    P = np.asarray(P, dtype=float)
    S_sol = self.solidus_entropy(P)
    S_liq = self.liquidus_entropy(P)
    T_sol = self._lookup_at_phase_boundary('temperature', P, 'solid')
    T_liq = self._lookup_at_phase_boundary('temperature', P, 'melt')
    T_fus = 0.5 * (T_sol + T_liq)
    return T_fus * np.maximum(S_liq - S_sol, 1.0)

def liquidus_entropy(self, P: npt.NDArray | float) -> npt.NDArray:
    """Liquidus entropy S_liq(P) [J/kg/K]."""
    return np.asarray(self._liquidus['interp'](P), dtype=float)

def liquidus_entropy_dP(self, P: npt.NDArray | float) -> npt.NDArray:
    """dS_liq/dP at the given pressure(s), in J/(kg·K·Pa)."""
    return np.asarray(self._liquidus['dinterp'](P), dtype=float)

def material_expansivity(
    self, P: npt.NDArray | float, S: npt.NDArray | float
) -> npt.NDArray:
    """Thermal expansivity [1/K] of the phases, without the phase-change term.

    The volume-weighted blend ``rho sum_i x_i alpha_i / rho_i`` of each phase's
    expansivity, with ``x_i`` the Lever-rule mass fractions and ``rho`` the harmonic
    mixture density of ``density``. Each ``alpha_i`` is from the thermal_exp tables when
    present, otherwise ``rho cp |dT/dP|_S / T`` of the phase at its own entropy.
    """
    phi = self._melt_fraction(P, S)
    melt = phi / self._lookup_at_phase_boundary('density', P, 'melt')
    volume = melt / (
        melt + (1.0 - phi) / self._lookup_at_phase_boundary('density', P, 'solid')
    )
    if self._has_alpha_tables:
        return self._lookup_phase_weighted('thermal_exp', P, S, melt_weight=volume)
    return self._lookup_phase_weighted(
        ('density', 'heat_capacity', 'dTdPs', 'temperature'),
        P,
        S,
        per_phase=lambda rho, cp, dTdPs, T: rho * cp * np.abs(dTdPs) / np.maximum(T, 1.0),
        melt_weight=volume,
    )

def melt_fraction(self, P: npt.NDArray | float, S: npt.NDArray | float) -> npt.NDArray:
    """Melt fraction phi from position between solidus and liquidus.

    phi = 0 for S <= S_sol, phi = 1 for S >= S_liq, linear between.
    """
    S = np.asarray(S, dtype=float)
    # phi feeds the live CVODE melt-fraction step cap directly, so flag a
    # non-finite or far-out-of-range S: the clip in the core maps +-inf to
    # a valid-looking 0 or 1 and passes NaN straight through.
    self._check_entropy_range(S, self.S_min, self.S_max, 'melt_fraction')
    return self._melt_fraction(P, S)

def solidus_entropy(self, P: npt.NDArray | float) -> npt.NDArray:
    """Solidus entropy S_sol(P) [J/kg/K]."""
    return np.asarray(self._solidus['interp'](P), dtype=float)

def solidus_entropy_dP(self, P: npt.NDArray | float) -> npt.NDArray:
    """dS_sol/dP at the given pressure(s), in J/(kg·K·Pa).

    Needed by the SPIDER-parity mixing-flux formula in
    ``entropy_state.update``. SPIDER computes Jmix via a bracket
    expression ``dS/dr − [φ dS_liq/dP + (1−φ) dS_sol/dP] dP/dr``
    (``energy.c::GetMixingHeatFlux`` lines 307-309). Exposing the
    phase-boundary P-derivatives here lets Aragog match SPIDER's
    formula exactly without resorting to un-truncated gphi
    arithmetic.
    """
    return np.asarray(self._solidus['dinterp'](P), dtype=float)

def specific_enthalpy(
    self,
    P: npt.NDArray | float,
    S: npt.NDArray | float,
) -> npt.NDArray:
    """EOS-consistent specific enthalpy h(P, S) [J/kg].

    Bilinear interpolation on the precomputed table built in
    ``_build_enthalpy_table``. Inputs outside the table P or S
    range are clamped to the nearest table edge and logged as a
    warning; pass ``strict_range=True`` to the constructor to raise
    a ``RuntimeError`` instead of returning the clamped value.

    Parameters
    ----------
    P : array or float
        Pressure [Pa].
    S : array or float
        Entropy [J/kg/K].

    Returns
    -------
    ndarray
        Specific enthalpy [J/kg], same shape as broadcast(P, S).
    """
    P = np.asarray(P, dtype=float)
    S = np.asarray(S, dtype=float)
    self._check_entropy_range(
        S,
        self._h_grid_S[0],
        self._h_grid_S[-1],
        'specific_enthalpy',
    )
    P_clamped = np.clip(P, self._h_grid_P[0], self._h_grid_P[-1])
    S_clamped = np.clip(S, self._h_grid_S[0], self._h_grid_S[-1])
    pts = np.column_stack([P_clamped.ravel(), S_clamped.ravel()])
    return self._h_interp(pts).reshape(P_clamped.shape)

def temperature(self, P: npt.NDArray | float, S: npt.NDArray | float) -> npt.NDArray:
    """Temperature T(P, S) [K]."""
    return self._lookup_phase_weighted('temperature', P, S)

def temperature_scalar(self, P: float, S: float) -> float:
    """Temperature T(P, S) for a single point, pure Python.

    Uses the same Lever Rule blending as ``temperature()`` but
    with pure-Python arithmetic throughout. Safe for use inside
    scipy.optimize callbacks on numpy >= 2.4.

    Parameters
    ----------
    P : float
        Pressure [Pa].
    S : float
        Entropy [J/kg/K].

    Returns
    -------
    float
        Temperature [K].
    """
    phi = self._melt_fraction_scalar(P, S)

    solid_table = self._tables['temperature_solid']
    melt_table = self._tables['temperature_melt']

    # Lever Rule: in the mushy zone, evaluate at phase boundary
    # entropies; outside, evaluate at actual S.
    if 0.0 < phi < 1.0:
        S_for_solid = self._solidus_entropy_scalar(P)
        S_for_melt = self._liquidus_entropy_scalar(P)
    else:
        S_for_solid = S
        S_for_melt = S

    # Only check a branch's S where it actually has nonzero weight
    # in the blend below (phi < 1.0 for solid, phi > 0.0 for melt).
    # Written as "not phi >= ..." so a NaN phi (comparisons always
    # False) still passes through the check instead of being masked.
    solid_used = not (phi >= 1.0)
    melt_used = not (phi <= 0.0)
    self._check_entropy_range(
        S_for_solid if solid_used else solid_table['S_list'][0],
        solid_table['S_list'][0],
        solid_table['S_list'][-1],
        'temperature_scalar (solid table lookup)',
    )
    self._check_entropy_range(
        S_for_melt if melt_used else melt_table['S_list'][0],
        melt_table['S_list'][0],
        melt_table['S_list'][-1],
        'temperature_scalar (melt table lookup)',
    )

    # Clamp to each table's own range
    S_sol_clamped = max(
        solid_table['S_list'][0], min(S_for_solid, solid_table['S_list'][-1])
    )
    S_melt_clamped = max(melt_table['S_list'][0], min(S_for_melt, melt_table['S_list'][-1]))
    P_sol_clamped = max(solid_table['P_list'][0], min(P, solid_table['P_list'][-1]))
    P_melt_clamped = max(melt_table['P_list'][0], min(P, melt_table['P_list'][-1]))

    val_solid = _interp_bilinear_scalar(
        solid_table['P_list'],
        solid_table['S_list'],
        solid_table['values'],
        P_sol_clamped,
        S_sol_clamped,
    )
    val_melt = _interp_bilinear_scalar(
        melt_table['P_list'],
        melt_table['S_list'],
        melt_table['values'],
        P_melt_clamped,
        S_melt_clamped,
    )

    # NaN-safe phase-weighted blend. phi > 0.0 and phi < 1.0 are both
    # False when phi is NaN, so without this check a NaN phi would
    # silently fall through to result = 0.0 instead of NaN.
    if math.isnan(phi):
        return math.nan
    result = 0.0
    if phi > 0.0:
        result += phi * val_melt
    if phi < 1.0:
        result += (1.0 - phi) * val_solid
    return result

def thermal_expansivity(
    self, P: npt.NDArray | float, S: npt.NDArray | float
) -> npt.NDArray:
    """Thermal expansivity alpha(P, S) [1/K].

    When P-S thermal_exp tables are available (generated alongside
    SPIDER tables), uses them directly for exact parity with SPIDER.
    Otherwise, derives alpha from the thermodynamic identity:
    alpha = rho * Cp * |dT/dP|_S| / T.
    """
    if self._has_alpha_tables:
        return self._lookup_phase_weighted('thermal_exp', P, S)

    T = self.temperature(P, S)
    rho = self.density(P, S)
    Cp = self.heat_capacity(P, S)
    dTdPs_val = self.dTdPs(P, S)
    alpha = rho * Cp * np.abs(dTdPs_val) / np.maximum(T, 1.0)
    return alpha

def thermal_expansivity_composite_blend(
    self,
    P: npt.NDArray | float,
    S: npt.NDArray | float,
    width: float = 0.01,
) -> npt.NDArray:
    """Composite two-phase thermal expansivity, matching SPIDER.

    Mirrors SPIDER's ``eos_composite.c:246`` formula:

        alpha_mix = (rho_solid - rho_melt) / (T_liq - T_sol) / rho

    where ``rho`` is the harmonic-mean composite density. This
    effective alpha captures the density change from the phase
    transition (Clausius-Clapeyron contribution), which dominates
    over the single-phase alpha by a factor of ~5 in the deep
    mushy zone. Outside the mushy band the result reduces to the
    pure-phase alpha from the P-S tables via tanh smoothing
    (SPIDER ``eos_composite.c:279``).

    Parameters
    ----------
    P : array or float
        Pressure [Pa].
    S : array or float
        Entropy [J/kg/K].
    width : float
        Tanh smoothing width across the phase boundary.

    Returns
    -------
    ndarray
        Thermal expansivity [1/K].
    """
    P = np.asarray(P, dtype=float)
    S = np.asarray(S, dtype=float)

    # Phase-boundary quantities at this pressure
    S_sol = np.asarray(self.solidus_entropy(P)).reshape(P.shape)
    S_liq = np.asarray(self.liquidus_entropy(P)).reshape(P.shape)
    T_sol = np.asarray(self._lookup_at_phase_boundary('temperature', P, 'solid')).reshape(
        P.shape
    )
    T_liq = np.asarray(self._lookup_at_phase_boundary('temperature', P, 'melt')).reshape(
        P.shape
    )
    rho_sol = np.asarray(self._lookup_at_phase_boundary('density', P, 'solid')).reshape(
        P.shape
    )
    rho_liq = np.asarray(self._lookup_at_phase_boundary('density', P, 'melt')).reshape(
        P.shape
    )

    # Composite density (harmonic mean, SPIDER eos_composite.c:236)
    phi = self._melt_fraction(P, S)
    inv_rho = phi / np.maximum(rho_liq, 1.0) + (1.0 - phi) / np.maximum(rho_sol, 1.0)
    rho_comp = 1.0 / np.maximum(inv_rho, 1e-30)

    # Mixed-phase alpha (SPIDER eos_composite.c:246)
    dT = T_liq - T_sol
    safe_dT = np.where(np.abs(dT) > 1e-10, dT, 1e-10)
    alpha_mix = (rho_sol - rho_liq) / safe_dT / rho_comp

    # Pure-phase alpha from the linear-blend tables
    alpha_pure = self._lookup_phase_weighted('thermal_exp', P, S)

    # SPIDER-parity tanh smoothing (same as heat_capacity_latent_blend)
    dS = S_liq - S_sol
    safe_dS = np.where(np.abs(dS) > 1e-10, dS, 1e-10)
    gphi = (S - S_sol) / safe_dS

    w = max(width, 1e-6)
    smooth = np.where(
        gphi > 0.5,
        0.5 * (1.0 + np.tanh((1.0 - gphi) / w)),
        0.5 * (1.0 + np.tanh(gphi / w)),
    )

    return smooth * alpha_mix + (1.0 - smooth) * alpha_pure

def capacitance(self) -> npt.NDArray:
    """Capacitance for the entropy equation: rho * T [kg K / m^3].

    The entropy equation is rho * T * dS/dt = div(F).
    """
    return self._density * self._temperature

def dTdPs(self) -> npt.NDArray:
    """Adiabatic temperature gradient dT/dP|_S [K/Pa]."""
    return self._dTdPs_val

def dTdrs(self) -> npt.NDArray:
    """Adiabatic temperature gradient dT/dr|_S [K/m].

    dT/dr|_S = -g * alpha * T / Cp
    """
    return -self._g * self._thermal_expansivity * self._temperature / self._heat_capacity

def delta_specific_volume(self) -> FloatOrArray:
    """Specific volume difference between solid and liquid [m^3/kg]."""
    if self._const_properties:
        return np.zeros_like(self._density)
    rho_s = self._eos._lookup_at_phase_boundary('density', self.pressure, 'solid')
    rho_l = self._eos._lookup_at_phase_boundary('density', self.pressure, 'melt')
    return 1.0 / np.maximum(rho_l, 1.0) - 1.0 / np.maximum(rho_s, 1.0)

def latent_heat(self) -> FloatOrArray:
    """Latent heat L(P) = T_fus × (S_liq - S_sol) [J/kg].

    P-dependent from EOS tables, matching SPIDER convention.
    """
    return self._latent_heat_val

def phase_boundary_densities(self) -> tuple[FloatOrArray, FloatOrArray]:
    """Solid and melt densities at the phase boundary at the node pressure [kg/m^3].

    Under ``const_properties`` there is no phase contrast and both equal
    the node density.
    """
    if self._const_properties:
        return self._density.copy(), self._density.copy()
    rho_s = self._eos._lookup_at_phase_boundary('density', self.pressure, 'solid')
    rho_l = self._eos._lookup_at_phase_boundary('density', self.pressure, 'melt')
    return rho_s, rho_l

def porosity(self) -> FloatOrArray:
    """Density-derived porosity ``(rho_s - rho) / (rho_s - rho_l)`` [1].

    This is the porosity :meth:`relative_velocity` uses, except under
    ``const_properties``, where there is no phase contrast, the mantle counts as
    fully liquid so this is one everywhere, and :meth:`relative_velocity` is zero.
    It is soft-clipped to [0, 1] with a width of 1e-3: its range is [-2.5e-7, 1]
    and it is 4.9975e-4 where the node density equals ``rho_s``. It is not the
    melt fraction: a solid node whose density is below ``rho_s`` can give any
    value, and where ``rho_s <= rho_l`` a mush is denser than ``rho_s``, so this
    is about 0 although the node holds melt. Do not use it for melt drainage
    where ``rho_s <= rho_l``: test the sign of the contrast first.
    """
    if self._const_properties:
        return np.ones_like(self._density)
    return self._porosity_from(*self.phase_boundary_densities())

def relative_velocity(self) -> FloatOrArray:
    """Melt-solid relative velocity for gravitational separation [m/s].

    Returns zero in const_properties mode (no phase contrast).

    Uses the Abe (1995) gravitational separation velocity with the
    three-regime permeability of Bower et al. (2018) section 2.1, based
    on the density-derived porosity (see :meth:`porosity`), matching
    SPIDER's GetGravitationalHeatFlux in energy.c.

    Regimes (F = K/porosity, the quantity multiplying delta_rho*g/eta):
    1. Blake-Kozeny-Carman (low porosity): F = d^2 por^2 / ((1-por)^2 * 1000)
    2. Rumpf-Gupte (intermediate): F = d^2 por^4.5 * (5/7)
    3. Stokes settling (high porosity): F = d^2 * 2/9

    The drag viscosity is set by ``separation_viscosity``: ``'melt'``
    uses the fixed single-phase liquid viscosity, matching SPIDER's
    ``GetGravitationalHeatFlux``; ``'mixture'`` uses the
    rheological-transition-blended bulk viscosity, so settling locks
    up through the same phi_rheo transition that sets the bulk
    rheology.

    The porosity denominator ``rho_s - rho_l`` is floored at about 1 kg/m^3,
    which acts where the tabulated densities cross (the melt locally denser
    or lighter by less than 1 kg/m^3). There the porosity is
    ``(rho_s - rho) / (1 kg/m^3)`` clipped to [0, 1]: about 0 where the node
    density is at or above ``rho_s``, 1 where it is more than 1 kg/m^3
    below, and intermediate in between. The separation direction is fixed by the
    non-negative contrast ``|rho_l - rho_s|`` in ``v_rel`` and cannot
    reverse.
    """
    if self._const_properties:
        return np.zeros_like(self._density)
    rho_s, rho_l = self.phase_boundary_densities()
    delta_rho = rho_l - rho_s  # typically negative (melt lighter)
    g = self._g
    # Explicit dispatch: a third mode added to SEPARATION_VISCOSITY_MODES
    # must fail here rather than silently fall back to 'melt'.
    if self._separation_viscosity == 'mixture':
        eta_l = self._viscosity_val
    elif self._separation_viscosity == 'melt':
        eta_l = self._visc_liquid
    else:
        raise ValueError(f'unhandled separation_viscosity {self._separation_viscosity!r}')

    porosity = self._porosity_from(rho_s, rho_l)
    F = mobility_function(porosity, self._grain_size)

    # Relative velocity: v = |delta_rho| * g * F / eta_liquid
    # Sign convention: positive = outward (melt rising, solid sinking)
    # Smoothed |delta_rho| via sqrt(x^2 + eps^2); eps tiny compared
    # to physical delta_rho ~ -500 kg/m^3, so bulk result unchanged.
    abs_drho = np.sqrt(delta_rho * delta_rho + 1.0e-12)
    v_rel = abs_drho * g * F / np.maximum(eta_l, 1e-10)

    return v_rel

def set_entropy(self, entropy: npt.NDArray) -> None:
    """Set the entropy profile."""
    self.entropy = np.asarray(entropy, dtype=float)

def set_pressure(self, pressure: npt.NDArray) -> None:
    """Set the pressure profile."""
    self.pressure = np.asarray(pressure, dtype=float)

def set_temperature(self, temperature: npt.NDArray) -> None:
    """Not used in entropy mode. Use set_entropy() instead."""
    raise NotImplementedError(
        'EntropyPhaseEvaluator uses set_entropy(), not set_temperature(). '
        'This evaluator cannot be used with the T-based solver.'
    )

def temperature(self) -> FloatOrArray:
    """Temperature from EOS lookup (not a state variable)."""
    return self._temperature

def update(self) -> None:
    """Recompute all cached properties from current (P, S).

    Two modes:
    - const_properties=True: analytical T(S) with constant rho, Cp,
      alpha, k, visc. Matches SPIDER's -use_const_properties. No EOS
      table lookups. phi=1 always (no phase transitions).
    - const_properties=False (default): single-pass EOS evaluation
      mirroring SPIDER's EOSEval_Composite_TwoPhase.
    """
    if self._const_properties:
        return self._update_const()
    return self._update_eos()