Vinet equation of state

Implements the third-order Rose-Vinet EOS used by White & Li (2025) and Boujibar et al. (2020). The pressure is expressed in terms of the Eulerian strain parameter \(f = (V/V_0)^{1/3}\) as \(P(f) = 3 K_0 f^{-2}(1-f)\exp[\eta(1-f)]\) with \(\eta = \tfrac{3}{2}(K_0' - 1)\), where \(K_0\) is the zero-pressure isothermal bulk modulus and \(K_0'\) its pressure derivative. Inversion to \(\rho(P)\) uses Brent root-finding. Material parameters cover Fe (solid: Smith+2018 with the White & Li 2025 light-element correction; liquid: Wicks+2018 Fe-7Si), MgSiO\(_3\) perovskite (Fei+2021), MgSiO\(_3\) post-perovskite (Sakai+2016), and peridotite (Stixrude & Lithgow-Bertelloni 2005).

`eos_vinet`

Rose-Vinet equation of state for planetary interior materials.

Implements the 3rd-order Vinet (Rose-Vinet) EOS used by White & Li (2025) and Boujibar et al. (2020) for computing density as a function of pressure. The Vinet EOS relates pressure to the Eulerian strain parameter f = (V/V_0)^(1/3):

P(f) = 3 K_0 f^{-2} (1 - f) exp[eta (1 - f)]

where eta = 3/2 (K_0' - 1), K_0 is the zero-pressure isothermal bulk modulus, and K_0' is its pressure derivative. Inversion to rho(P) uses Brent root-finding.

Material parameters from: - Fe (solid): Smith et al. (2018), with 12.5% thermal+light-element density reduction following White & Li (2025). - Fe (liquid): Wicks et al. (2018), Fe-7Si alloy. - MgSiO3 perovskite: Fei et al. (2021), via Boujibar et al. (2020). - MgSiO3 post-perovskite: Sakai et al. (2016). - Peridotite: Stixrude & Lithgow-Bertelloni (2005).

References

Boujibar, A., Driscoll, P. & Fei, Y. (2020). JGRP, 125, e2019JE006124. White, N. I. & Li, J. (2025). JGRP, 130, e2024JE008550. Smith, R. F. et al. (2018). Nature Astronomy, 2, 452-458. Fei, Y. et al. (2021). Phys. Rev. Lett., 127, 080501. Wicks, J. K. et al. (2018). Science Advances, 4, eaao5864. Sakai, T. et al. (2016). Geophys. Res. Lett., 43, 7648-7656.

`get_vinet_density(pressure, material_key)`

Compute density from the Rose-Vinet EOS.

Parameters:

Name	Type	Description	Default
`pressure`	`float`	Pressure in Pa.	required
`material_key`	`str`	One of the keys in VINET_MATERIALS (e.g. 'iron', 'MgSiO3').	required

Returns:

Type	Description
`float or None`	Density in kg/m³, or None for invalid input.

Raises:

Type	Description
`ValueError`	If material_key is not recognized.

Source code in src/zalmoxis/eos_vinet.py

def get_vinet_density(pressure: float, material_key: str) -> float | None:
    """Compute density from the Rose-Vinet EOS.

    Parameters
    ----------
    pressure : float
        Pressure in Pa.
    material_key : str
        One of the keys in VINET_MATERIALS (e.g. 'iron', 'MgSiO3').

    Returns
    -------
    float or None
        Density in kg/m³, or None for invalid input.

    Raises
    ------
    ValueError
        If material_key is not recognized.
    """
    if material_key not in VINET_MATERIALS:
        raise ValueError(
            f"Unknown Vinet material '{material_key}'. Valid keys: {sorted(VALID_VINET_KEYS)}"
        )

    if np.isnan(pressure):
        return None

    mat = VINET_MATERIALS[material_key]
    rho_0 = mat['rho_0']
    K_0 = mat['K_0']
    K_prime = mat['K_prime']
    thermal = mat['thermal_correction']

    if pressure <= 0:
        return float(rho_0 * thermal)

    if pressure > P_MAX_VINET:
        logger.warning(
            'Pressure %.2e Pa exceeds Vinet validity limit %.2e Pa for material %s. Clamping.',
            pressure,
            P_MAX_VINET,
            material_key,
        )
        pressure = P_MAX_VINET

    eta = 1.5 * (K_prime - 1.0)

    # Find f such that P_vinet(f) = pressure, using Brent root-finding.
    # f = 1 at P = 0, f -> 0 at P -> inf.
    # Search bracket: f in [f_min, 1 - epsilon]
    f_min = 0.01  # extreme compression (rho ~ 1e6 * rho_0)

    def residual(f):
        return _vinet_pressure(f, K_0, eta) - pressure

    # Check bracket validity
    if residual(1.0 - 1e-12) > 0:
        # P_target < P(f~1), which means P ~ 0
        return float(rho_0 * thermal)

    if residual(f_min) < 0:
        # P_target > P(f_min), need even more compression
        f_min = 0.001

    try:
        f_solution = brentq(residual, f_min, 1.0 - 1e-12, xtol=1e-14, rtol=1e-12)
    except ValueError:
        logger.warning(
            'Vinet root-finding failed for P=%.2e Pa, material=%s. Bracket: [%.4e, %.4e].',
            pressure,
            material_key,
            residual(f_min),
            residual(1.0 - 1e-12),
        )
        return None

    # rho = rho_0 / f^3, with thermal correction
    density = rho_0 / f_solution**3 * thermal

    return float(density)