Zalmoxis.melting curves

def get_melting_curve_function(curve_id):
    """Return a callable melting curve f(P [Pa]) -> T [K] for the given identifier.

    Parameters
    ----------
    curve_id : str
        Melting curve identifier. See module docstring for valid values.

    Returns
    -------
    callable
        Function accepting pressure in Pa (scalar or array) and returning
        temperature in K.

    Raises
    ------
    ValueError
        If ``curve_id`` is not recognized.
    """
    if curve_id == 'Monteux16-solidus':
        return monteux16_solidus

    elif curve_id == 'Monteux16-liquidus-A-chondritic':

        def _liq(P):
            return monteux16_liquidus(P, model='A-chondritic')

        return _liq

    elif curve_id == 'Monteux16-liquidus-F-peridotitic':

        def _liq(P):
            return monteux16_liquidus(P, model='F-peridotitic')

        return _liq

    elif curve_id == 'Stixrude14-solidus':
        return stixrude14_solidus

    elif curve_id == 'Stixrude14-liquidus':
        return stixrude14_liquidus

    elif curve_id == 'PALEOS-liquidus':
        return paleos_liquidus

    elif curve_id == 'Monteux600-solidus-tabulated':
        return _load_tabulated_curve(
            os.path.join(
                get_zalmoxis_root(), 'data', 'melting_curves_Monteux-600', 'solidus.dat'
            )
        )

    elif curve_id == 'Monteux600-liquidus-tabulated':
        return _load_tabulated_curve(
            os.path.join(
                get_zalmoxis_root(), 'data', 'melting_curves_Monteux-600', 'liquidus.dat'
            )
        )

    elif curve_id == 'Anzellini13-iron':
        return iron_melting_anzellini13

    elif curve_id == 'Sinmyo19-iron':
        return iron_melting_sinmyo19

    else:
        all_valid = sorted(VALID_SOLIDUS | VALID_LIQUIDUS | VALID_IRON_MELTING)
        raise ValueError(f"Unknown melting curve '{curve_id}'. Valid values: {all_valid}")

def get_solidus_liquidus_functions(
    solidus_id='Stixrude14-solidus',
    liquidus_id='Stixrude14-liquidus',
):
    """Load solidus and liquidus functions by config identifier.

    Parameters
    ----------
    solidus_id : str
        Solidus curve identifier.
    liquidus_id : str
        Liquidus curve identifier.

    Returns
    -------
    tuple of callable
        ``(solidus_func, liquidus_func)``
    """
    return get_melting_curve_function(solidus_id), get_melting_curve_function(liquidus_id)

def iron_melting_anzellini13(P):
    """Iron melting temperature from Anzellini et al. (2013).

    Composite Simon-Glatzel law matching the PALEOS iron phase diagram:

    - Below 98.5 GPa (gamma-Fe / liquid): Eq. 2 from Anzellini+2013
    - Above 98.5 GPa (epsilon-Fe / liquid): Eq. 3 from Anzellini+2013

    Note: the gamma-Fe branch is parameterized relative to (P0=5.2 GPa,
    T0=1991 K). Below P0, the formula extrapolates to lower temperatures
    (non-physical: melting should increase with P). For P < 1 GPa, this
    function clamps to the 1 atm melting point of pure iron (1811 K).
    The epsilon-Fe branch (P > 98.5 GPa) is the relevant one for
    planetary cores.

    Parameters
    ----------
    P : float or array-like
        Pressure [Pa].

    Returns
    -------
    float or ndarray
        Melting temperature [K].

    References
    ----------
    Anzellini, S. et al. (2013). Science, 340, 464-466.
    """
    P_arr = np.atleast_1d(np.asarray(P, dtype=float))
    P_GPa = P_arr / 1e9

    P0_GPa = 5.2  # Reference pressure [GPa]
    T0 = 1991.0  # Reference temperature [K]
    Pt_GPa = 98.5  # Triple point pressure [GPa]
    Tt = 3712.0  # Triple point temperature [K]
    T_1atm = 1811.0  # 1 atm melting point of pure iron [K]

    T = np.where(
        P_GPa < Pt_GPa,
        T0 * ((P_GPa - P0_GPa) / 27.39 + 1.0) ** (1.0 / 2.38),
        Tt * ((P_GPa - Pt_GPa) / 161.2 + 1.0) ** (1.0 / 1.72),
    )
    # Clamp to 1 atm melting point for low pressures where the
    # Simon-Glatzel extrapolation is non-physical
    T = np.maximum(T, T_1atm)
    return float(T[0]) if np.ndim(P) == 0 else T

def iron_melting_sinmyo19(P):
    """Iron melting temperature from Sinmyo et al. (2019).

    Single Simon-Glatzel law valid to ~290 GPa:
    T = T* * (P / a + 1)^b

    Parameters
    ----------
    P : float or array-like
        Pressure [Pa].

    Returns
    -------
    float or ndarray
        Melting temperature [K].

    References
    ----------
    Sinmyo, R. et al. (2019). EPSL, 510, 45-52.
    """
    P_arr = np.atleast_1d(np.asarray(P, dtype=float))

    T_star = 1811.0  # Reference temperature [K]
    a = 134.69e9  # Reference pressure [Pa]
    b = 0.93  # Exponent [-]

    T = T_star * (P_arr / a + 1.0) ** b
    return float(T[0]) if np.ndim(P) == 0 else T

def monteux16_liquidus(P, model='A-chondritic'):
    """Monteux+2016 liquidus temperature.

    Piecewise: Eq. 11 (P <= 20 GPa) and Eq. 13 (P > 20 GPa).

    Parameters
    ----------
    P : float or array-like
        Pressure in Pa (must be >= 0).
    model : str
        ``'A-chondritic'`` or ``'F-peridotitic'``.

    Returns
    -------
    float or ndarray
        Liquidus temperature in K.
    """
    P_arr = np.atleast_1d(np.asarray(P, dtype=float))
    T = np.empty_like(P_arr)
    lo = P_arr <= _P_TRANSITION
    hi = ~lo
    T[lo] = _liquidus_low(P_arr[lo])

    if model == 'A-chondritic':
        T[hi] = _liquidus_high_A(P_arr[hi])
    elif model == 'F-peridotitic':
        T[hi] = _liquidus_high_F(P_arr[hi])
    else:
        raise ValueError(f"Unknown model '{model}'. Use 'A-chondritic' or 'F-peridotitic'.")

    return float(T[0]) if np.ndim(P) == 0 else T

def monteux16_solidus(P):
    """Monteux+2016 solidus temperature.

    Piecewise: Eq. 10 (P <= 20 GPa) and Eq. 12 (P > 20 GPa).
    Both composition models share the same solidus.

    Parameters
    ----------
    P : float or array-like
        Pressure in Pa (must be >= 0).

    Returns
    -------
    float or ndarray
        Solidus temperature in K.
    """
    P_arr = np.atleast_1d(np.asarray(P, dtype=float))
    T = np.empty_like(P_arr)
    lo = P_arr <= _P_TRANSITION
    hi = ~lo
    T[lo] = _solidus_low(P_arr[lo])
    T[hi] = _solidus_high(P_arr[hi])
    return float(T[0]) if np.ndim(P) == 0 else T

def paleos_liquidus(P):
    """PALEOS MgSiO3 liquidus from Belonoshko+2005 / Fei+2021.

    Piecewise Simon-Glatzel fit:

    - P < 2.55 GPa: T = 1831 * (1 + P/4.6)^0.33  (Belonoshko+2005)
    - P >= 2.55 GPa: T = 6000 * (P/140)^0.26      (Fei+2021)

    This is the melting curve used internally by the PALEOS unified EOS
    tables. Using it for the mushy zone calculation ensures consistency
    between the table's phase boundaries and the derived solidus.

    Parameters
    ----------
    P : float or array-like
        Pressure in Pa (must be >= 0).

    Returns
    -------
    float or ndarray
        Liquidus temperature in K.
    """
    # Scalar fast path: in the numpy density Picard loop this function is
    # called once per cell per inner iteration on scalar P; the np.atleast_1d
    # / np.asarray / np.where / float() machinery dominates the trivial
    # piecewise-power-law math. Branching on scalar P first keeps the
    # per-call overhead negligible relative to a real CHILI solve.
    if isinstance(P, (int, float)) or (hasattr(P, 'ndim') and P.ndim == 0):
        Pf = float(P)
        if Pf <= 0.0:
            return 0.0
        P_GPa = Pf * 1e-9
        if P_GPa < _PALEOS_P0_GPA:
            return 1831.0 * (1.0 + P_GPa / 4.6) ** 0.33
        return 6000.0 * (P_GPa / 140.0) ** 0.26

    P_arr = np.atleast_1d(np.asarray(P, dtype=float))
    P_GPa = P_arr * 1e-9
    T = np.where(
        P_GPa < _PALEOS_P0_GPA,
        1831.0 * (1.0 + P_GPa / 4.6) ** 0.33,
        6000.0 * (P_GPa / 140.0) ** 0.26,
    )
    # Guard P=0: avoid 0^0.26 = NaN
    T = np.where(P_arr > 0, T, 0.0)
    return T

def stixrude14_liquidus(P):
    """Silicate (MgSiO3) liquidus from Stixrude (2014) Eq. 1.9.

    Simon-like power law fitted to Lindemann melting theory:
    T_rock = 5400 K * (P / 140 GPa)^0.480

    Defined for all P > 0. At P = 0, returns 0 K (extrapolation).

    Parameters
    ----------
    P : float or array-like
        Pressure in Pa (must be >= 0).

    Returns
    -------
    float or ndarray
        Liquidus temperature in K.
    """
    # Scalar fast path. Skip np.atleast_1d / np.asarray / np.where
    # dispatch overhead for scalar P. The solver's RHS calls this per
    # cell and per step (millions of calls per solve); the array branch
    # below is dispatched only when P is array-like. Scalar math is
    # identical to the array np.where branch for P > 0 and the P == 0
    # branch for P == 0.
    if isinstance(P, (int, float, np.floating, np.integer)):
        return _STIX14_T_REF * (P / _STIX14_P_REF) ** _STIX14_EXPONENT if P > 0 else 0.0
    P_arr = np.atleast_1d(np.asarray(P, dtype=float))
    # Guard P=0: the power law gives 0 K at P=0
    T = np.where(
        P_arr > 0,
        _STIX14_T_REF * (P_arr / _STIX14_P_REF) ** _STIX14_EXPONENT,
        0.0,
    )
    return float(T[0]) if np.ndim(P) == 0 else T

def stixrude14_solidus(P):
    """Silicate solidus from Stixrude (2014) Eqs. 1.9 + 1.10.

    The solidus is the pure-substance liquidus (Eq. 1.9) depressed by
    the cryoscopic equation (Eq. 1.10) with x0 = 0.79 (mole fraction
    of pure MgSiO3 in an Earth-like mantle composition):

    T_solidus = T_liquidus * (1 - ln(x0))^{-1}

    With x0 = 0.79 this gives ~4370 K at 140 GPa, approximately
    consistent with the experimentally measured solidus of an Earth-like
    mantle composition (~4100 K at 140 GPa, Fiquet et al. 2010). The
    ~270 K offset is within the uncertainty of the cryoscopic model.

    Parameters
    ----------
    P : float or array-like
        Pressure in Pa (must be >= 0).

    Returns
    -------
    float or ndarray
        Solidus temperature in K.
    """
    return stixrude14_liquidus(P) * _STIX14_CRYO_FACTOR