From 93ec10bc21d50094b2b5c30c859bc086410d7e0b Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Mon, 25 Nov 2024 12:28:51 +0000 Subject: [PATCH] added SPOCK EoS --- burnman/eos/__init__.py | 1 + burnman/eos/helper.py | 3 + burnman/eos/spock.py | 229 ++++++++++++++++++++++++++++ misc/benchmarks/spock_benchmarks.py | 56 +++++++ misc/ref/spock_benchmarks.py.out | 8 + 5 files changed, 297 insertions(+) create mode 100644 burnman/eos/spock.py create mode 100644 misc/benchmarks/spock_benchmarks.py create mode 100644 misc/ref/spock_benchmarks.py.out diff --git a/burnman/eos/__init__.py b/burnman/eos/__init__.py index 16fecdac..ff59474c 100644 --- a/burnman/eos/__init__.py +++ b/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 diff --git a/burnman/eos/helper.py b/burnman/eos/helper.py index 59972ee2..f60afd42 100644 --- a/burnman/eos/helper.py +++ b/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": diff --git a/burnman/eos/spock.py b/burnman/eos/spock.py new file mode 100644 index 00000000..98c6b580 --- /dev/null +++ b/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) diff --git a/misc/benchmarks/spock_benchmarks.py b/misc/benchmarks/spock_benchmarks.py new file mode 100644 index 00000000..84413183 --- /dev/null +++ b/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() diff --git a/misc/ref/spock_benchmarks.py.out b/misc/ref/spock_benchmarks.py.out new file mode 100644 index 00000000..a67f2a6b --- /dev/null +++ b/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