added SPOCK EoS

burnman/eos/__init__.py

@@ -27,6 +27,7 @@
 from .morse_potential import Morse
 from .reciprocal_kprime import RKprime
 from .macaw import MACAW
+from .spock import SPOCK
 from .dks_liquid import DKS_L
 from .dks_solid import DKS_S
 from .aa import AA

burnman/eos/helper.py

@@ -12,6 +12,7 @@
 from . import birch_murnaghan_4th as bm4
 from . import modified_tait as mt
 from . import macaw
+from . import spock
 from . import dks_liquid
 from . import dks_solid
 from . import hp
@@ -73,6 +74,8 @@ def create(method):
             return mt.MT()
         elif method == "macaw":
             return macaw.MACAW()
+        elif method == "spock":
+            return spock.SPOCK()
         elif method == "hp98":
             return hp.HP98()
         elif method == "hp_tmt":

burnman/eos/spock.py

@@ -0,0 +1,229 @@
+from __future__ import absolute_import
+
+# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit
+# for the Earth and Planetary Sciences.
+# Copyright (C) 2012 - 2024 by the BurnMan team, released under the GNU
+# GPL v2 or later.
+
+
+import scipy.optimize as opt
+from scipy.special import gamma, gammainc, exp1
+from . import equation_of_state as eos
+import warnings
+import numpy as np
+
+
+# Try to import the jit from numba.  If it is
+# not available, just go with the standard
+# python interpreter
+try:
+    from numba import jit
+except ImportError:
+
+    def jit(fn):
+        return fn
+
+
+def gammaincc(a, x):
+    """
+    An implementation of the non-regularised upper incomplete gamma
+    function. Computed using the relationship with the regularised
+    lower incomplete gamma function (scipy.special.gammainc).
+    Uses the recurrence relation wherever a<0.
+    """
+    n = int(-np.floor(a))
+    if n > 0:
+        a = a + n
+        u_gamma = exp1(x) if a == 0 else (1.0 - gammainc(a, x)) * gamma(a)
+
+        for _ in range(n):
+            a = a - 1.0
+            u_gamma = (u_gamma - np.power(x, a) * np.exp(-x)) / a
+        return u_gamma
+    else:
+        return (1.0 - gammainc(a, x)) * gamma(a)
+
+
+@jit(nopython=True)
+def make_params(K_0, Kp_0, Kp_inf, Kdp_0):
+    dKpdlnV_zero = -Kdp_0 * K_0
+    c = Kp_inf
+    a = dKpdlnV_zero / (Kp_0 - Kp_inf)
+    b = (Kp_0 - Kp_inf) / a
+    return a, b, c
+
+
+class SPOCK(eos.EquationOfState):
+    """
+    Class for the Scaled Power Of Compression K-prime equation of state.
+    This equation is derived from the assumption that K' = b*(V/V_0)^a.
+
+    This equation of state has no temperature dependence.
+    """
+
+    def isothermal_bulk_modulus_reuss(self, pressure, temperature, volume, params):
+        """
+        Returns isothermal bulk modulus :math:`K_T` :math:`[Pa]` as a function of pressure :math:`[Pa]`,
+        temperature :math:`[K]` and volume :math:`[m^3]`.
+        """
+        ai, bi, ci = make_params(
+            params["K_0"], params["Kprime_0"], params["Kprime_inf"], params["Kdprime_0"]
+        )
+
+        lnVrel = np.log(volume / params["V_0"])
+        return params["K_0"] * np.exp(-bi * (np.exp(ai * lnVrel) - 1.0) - ci * lnVrel)
+
+    def volume(self, pressure, temperature, params):
+        """
+        Get the Vinet volume at a reference temperature for a given
+        pressure :math:`[Pa]`. Returns molar volume in :math:`[m^3]`
+        """
+
+        def delta_pressure(x):
+            return self.pressure(0.0, x, params) - pressure
+
+        V = opt.brentq(delta_pressure, 0.1 * params["V_0"], 1.5 * params["V_0"])
+        return V
+
+    def pressure(self, temperature, volume, params):
+        """
+        Returns pressure :math:`[Pa]` as a function of volume :math:`[m^3]`.
+        """
+        ai, bi, ci = make_params(
+            params["K_0"], params["Kprime_0"], params["Kprime_inf"], params["Kdprime_0"]
+        )
+        lnVrel = np.log(volume / params["V_0"])
+        return params["P_0"] + (
+            params["K_0"]
+            * np.exp(bi)
+            / ai
+            * np.power(bi, ci / ai)
+            * (
+                gammaincc(
+                    -ci / ai,
+                    bi * np.exp(ai * lnVrel),
+                )
+                - gammaincc(-ci / ai, bi)
+            )
+        )
+
+    def molar_internal_energy(self, pressure, temperature, volume, params):
+        """
+        Returns the internal energy :math:`\\mathcal{E}` of the mineral. :math:`[J/mol]`
+        """
+        ai, bi, ci = make_params(
+            params["K_0"], params["Kprime_0"], params["Kprime_inf"], params["Kdprime_0"]
+        )
+        lnVrel = np.log(volume / params["V_0"])
+        f = (
+            -params["V_0"]
+            * params["K_0"]
+            * np.exp(bi)
+            / ai
+            * np.power(bi, (ci - 1.0) / ai)
+        )
+
+        Vrel = np.exp(lnVrel)
+        I1 = (
+            np.power(bi, 1.0 / ai)
+            * Vrel
+            * (
+                gammaincc(
+                    -ci / ai,
+                    bi * np.exp(ai * lnVrel),
+                )
+                - gammaincc(-ci / ai, bi)
+            )
+        )
+        I2 = gammaincc(
+            (1.0 - ci) / ai,
+            bi * np.exp(ai * lnVrel),
+        ) - gammaincc((1.0 - ci) / ai, bi)
+
+        return params["E_0"] + params["P_0"] * (volume - params["V_0"]) + f * (I1 - I2)
+
+    def gibbs_free_energy(self, pressure, temperature, volume, params):
+        """
+        Returns the Gibbs free energy :math:`\\mathcal{G}` of the mineral. :math:`[J/mol]`
+        """
+        return (
+            self.molar_internal_energy(pressure, temperature, volume, params)
+            + pressure * volume
+        )
+
+    def isentropic_bulk_modulus_reuss(self, pressure, temperature, volume, params):
+        """
+        Returns adiabatic bulk modulus :math:`K_s` of the mineral. :math:`[Pa]`.
+        """
+        return self.isothermal_bulk_modulus_reuss(pressure, temperature, volume, params)
+
+    def shear_modulus(self, pressure, temperature, volume, params):
+        """
+        Returns shear modulus :math:`G` of the mineral. :math:`[Pa]`
+        """
+        return 1.0e99
+
+    def entropy(self, pressure, temperature, volume, params):
+        """
+        Returns the molar entropy :math:`\\mathcal{S}` of the mineral. :math:`[J/K/mol]`
+        """
+        return 0.0
+
+    def molar_heat_capacity_v(self, pressure, temperature, volume, params):
+        """
+        Since this equation of state does not contain temperature effects, return a very small number. :math:`[J/K/mol]`
+        """
+        return 1.0e-99
+
+    def molar_heat_capacity_p(self, pressure, temperature, volume, params):
+        """
+        Since this equation of state does not contain temperature effects, return a very small number. :math:`[J/K/mol]`
+        """
+        return 1.0e-99
+
+    def thermal_expansivity(self, pressure, temperature, volume, params):
+        """
+        Since this equation of state does not contain temperature effects, return zero. :math:`[1/K]`
+        """
+        return 0.0
+
+    def grueneisen_parameter(self, pressure, temperature, volume, params):
+        """
+        Since this equation of state does not contain temperature effects, return zero. :math:`[unitless]`
+        """
+        return 0.0
+
+    def validate_parameters(self, params):
+        """
+        Check for existence and validity of the parameters.
+        The value for :math:`K'_{\\infty}` is thermodynamically bounded
+        between 5/3 and :math:`K'_0` :cite:`StaceyDavis2004`.
+        """
+
+        if "E_0" not in params:
+            params["E_0"] = 0.0
+        if "P_0" not in params:
+            params["P_0"] = 1.0e5
+
+        # Check that all the required keys are in the dictionary
+        expected_keys = ["V_0", "K_0", "Kprime_0", "Kdprime_0", "Kprime_inf"]
+        for k in expected_keys:
+            if k not in params:
+                raise KeyError("params object missing parameter : " + k)
+
+        # Finally, check that the values are reasonable.
+        if params["P_0"] < 0.0:
+            warnings.warn("Unusual value for P_0", stacklevel=2)
+        if params["V_0"] < 1.0e-7 or params["V_0"] > 1.0e-3:
+            warnings.warn("Unusual value for V_0", stacklevel=2)
+        if params["K_0"] < 1.0e9 or params["K_0"] > 1.0e13:
+            warnings.warn("Unusual value for K_0", stacklevel=2)
+        if params["Kprime_0"] < 0.0 or params["Kprime_0"] > 10.0:
+            warnings.warn("Unusual value for Kprime_0", stacklevel=2)
+        if params["Kdprime_0"] > 0.0:
+            warnings.warn("Unusual value for Kdprime_0", stacklevel=2)
+        if (
+            params["Kprime_inf"] < 5.0 / 3.0
+            or params["Kprime_inf"] > params["Kprime_0"]
+        ):
+            warnings.warn("Unusual value for Kprime_inf", stacklevel=2)

misc/benchmarks/spock_benchmarks.py

@@ -0,0 +1,56 @@
+from __future__ import absolute_import
+
+# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit
+# for the Earth and Planetary Sciences.
+# Copyright (C) 2012 - 2024 by the BurnMan team, released under the GNU
+# GPL v2 or later.
+
+from burnman.tools.eos import check_eos_consistency
+from burnman import Mineral
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.image as mpimg
+
+HMX_params = {
+    "P_0": 1.0e5,
+    "V_0": 1.0e-6,  # arbitrary value
+    "K_0": 15.22e9,
+    "Kprime_0": 7.54,
+    "Kdprime_0": -7.54 / 15.22e9,
+    "Kprime_inf": 2.63,
+    "molar_mass": 0.296155,
+    "equation_of_state": "spock",
+}
+
+HMX = Mineral(HMX_params)
+
+if check_eos_consistency(HMX, tol=1.0e-5, including_shear_properties=False):
+    print("The SPOCK EoS is internally consistent.\n")
+
+pressures = np.linspace(0.0, 100.0e9, 6)
+temperatures = 0.0 + 0.0 * pressures
+V, K_T = HMX.evaluate(["V", "K_T"], pressures, temperatures)
+
+for i in range(6):
+    print(
+        f"{pressures[i]/1.e9:3.0f} GPa: "
+        f"V/V_0 = {V[i]/HMX_params['V_0']:.3f}, "
+        f"K_T = {K_T[i]/1.e9:6.2f} GPa"
+    )
+
+
+fig1 = mpimg.imread("figures/Lozano_Aslam_2022_Fig6c_HMX.png")
+plt.imshow(fig1, extent=[0.0, 100.0, 0, 500.0], aspect="auto")
+
+for a in [1.0, 1.5, 2.0]:
+    HMX_params["Kdprime_0"] = -a * HMX_params["Kprime_0"] / HMX_params["K_0"]
+
+    pressures = np.linspace(0.0, 100.0e9, 101)
+    temperatures = 0.0 + 0.0 * pressures
+    K_T = HMX.evaluate(["K_T"], pressures, temperatures)[0]
+
+    plt.plot(pressures / 1.0e9, K_T / 1.0e9, linestyle=":", label=f"SPOCK {a}")
+
+plt.ylim(0.0, 500.0)
+plt.legend(loc=(0.025, 0.5))
+plt.show()

misc/ref/spock_benchmarks.py.out

@@ -0,0 +1,8 @@
+The SPOCK EoS is internally consistent.
+
+  0 GPa: V/V_0 = 1.000, K_T =  15.22 GPa
+ 20 GPa: V/V_0 = 0.711, K_T = 137.70 GPa
+ 40 GPa: V/V_0 = 0.638, K_T = 243.48 GPa
+ 60 GPa: V/V_0 = 0.596, K_T = 342.73 GPa
+ 80 GPa: V/V_0 = 0.566, K_T = 437.90 GPa
+100 GPa: V/V_0 = 0.543, K_T = 530.17 GPa