Source code for burnman.eos.slb

# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit
# for the Earth and Planetary Sciences
# Copyright (C) 2012 - 2017 by the BurnMan team, released under the GNU
# GPL v2 or later.

from __future__ import absolute_import

import numpy as np
import scipy.optimize as opt
import warnings

# 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


from . import birch_murnaghan as bm
from . import debye
from . import equation_of_state as eos
from ..utils.math import bracket


@jit(nopython=True)
def _grueneisen_parameter_fast(V_0, volume, gruen_0, q_0):
    """global function with plain parameters so jit will work"""
    x = V_0 / volume
    f = 1.0 / 2.0 * (pow(x, 2.0 / 3.0) - 1.0)
    a1_ii = 6.0 * gruen_0  # EQ 47
    a2_iikk = -12.0 * gruen_0 + 36.0 * gruen_0 * gruen_0 - 18.0 * q_0 * gruen_0  # EQ 47
    nu_o_nu0_sq = 1.0 + a1_ii * f + (1.0 / 2.0) * a2_iikk * f * f  # EQ 41
    return 1.0 / 6.0 / nu_o_nu0_sq * (2.0 * f + 1.0) * (a1_ii + a2_iikk * f)


@jit(nopython=True)
def _delta_pressure(
    x, pressure, temperature, V_0, T_0, Debye_0, n, a1_ii, a2_iikk, b_iikk, b_iikkmm
):
    f = 0.5 * (pow(V_0 / x, 2.0 / 3.0) - 1.0)
    nu_o_nu0_sq = 1.0 + a1_ii * f + 1.0 / 2.0 * a2_iikk * f * f
    debye_temperature = Debye_0 * np.sqrt(nu_o_nu0_sq)
    E_th = debye.thermal_energy(
        temperature, debye_temperature, n
    )  # thermal energy at temperature T
    E_th_ref = debye.thermal_energy(
        T_0, debye_temperature, n
    )  # thermal energy at reference temperature
    nu_o_nu0_sq = 1.0 + a1_ii * f + (1.0 / 2.0) * a2_iikk * f * f  # EQ 41
    gr = 1.0 / 6.0 / nu_o_nu0_sq * (2.0 * f + 1.0) * (a1_ii + a2_iikk * f)

    return (
        (1.0 / 3.0)
        * (pow(1.0 + 2.0 * f, 5.0 / 2.0))
        * ((b_iikk * f) + (0.5 * b_iikkmm * f * f))
        + gr * (E_th - E_th_ref) / x
        - pressure
    )  # EQ 21


[docs] class SLBBase(eos.EquationOfState): """ Base class for the finite strain-Mie-Grueneiesen-Debye equation of state detailed in :cite:`Stixrude2005`. For the most part the equations are all third order in strain, but see further the :class:`burnman.slb.SLB2` and :class:`burnman.slb.SLB3` classes. """ def _debye_temperature(self, x, params): """ Finite strain approximation for Debye Temperature [K] x = ref_vol/vol """ f = 1.0 / 2.0 * (pow(x, 2.0 / 3.0) - 1.0) a1_ii = 6.0 * params["grueneisen_0"] # EQ 47 a2_iikk = ( -12.0 * params["grueneisen_0"] + 36.0 * pow(params["grueneisen_0"], 2.0) - 18.0 * params["q_0"] * params["grueneisen_0"] ) # EQ 47 nu_o_nu0_sq = 1.0 + a1_ii * f + 1.0 / 2.0 * a2_iikk * f * f if nu_o_nu0_sq > 0.0: return params["Debye_0"] * np.sqrt(nu_o_nu0_sq) else: raise Exception( f"This volume (V = {1./x:.2f}*V_0) exceeds the " "valid range of the thermal " "part of the slb equation of state." )
[docs] def volume_dependent_q(self, x, params): """ Finite strain approximation for :math:`q`, the isotropic volume strain derivative of the grueneisen parameter. """ f = 1.0 / 2.0 * (pow(x, 2.0 / 3.0) - 1.0) a1_ii = 6.0 * params["grueneisen_0"] # EQ 47 a2_iikk = ( -12.0 * params["grueneisen_0"] + 36.0 * pow(params["grueneisen_0"], 2.0) - 18.0 * params["q_0"] * params["grueneisen_0"] ) # EQ 47 nu_o_nu0_sq = 1.0 + a1_ii * f + (1.0 / 2.0) * a2_iikk * f * f # EQ 41 gr = 1.0 / 6.0 / nu_o_nu0_sq * (2.0 * f + 1.0) * (a1_ii + a2_iikk * f) # avoids divide by zero if grueneisen_0 = 0. if np.abs(params["grueneisen_0"]) < 1.0e-10: q = 1.0 / 9.0 * (18.0 * gr - 6.0) else: q = ( 1.0 / 9.0 * ( 18.0 * gr - 6.0 - 1.0 / 2.0 / nu_o_nu0_sq * (2.0 * f + 1.0) * (2.0 * f + 1.0) * a2_iikk / gr ) ) return q
def _isotropic_eta_s(self, x, params): """ Finite strain approximation for :math:`eta_{s0}`, the isotropic shear strain derivative of the grueneisen parameter. """ f = 1.0 / 2.0 * (pow(x, 2.0 / 3.0) - 1.0) a2_s = -2.0 * params["grueneisen_0"] - 2.0 * params["eta_s_0"] # EQ 47 a1_ii = 6.0 * params["grueneisen_0"] # EQ 47 a2_iikk = ( -12.0 * params["grueneisen_0"] + 36.0 * pow(params["grueneisen_0"], 2.0) - 18.0 * params["q_0"] * params["grueneisen_0"] ) # EQ 47 nu_o_nu0_sq = 1.0 + a1_ii * f + (1.0 / 2.0) * a2_iikk * pow(f, 2.0) # EQ 41 gr = 1.0 / 6.0 / nu_o_nu0_sq * (2.0 * f + 1.0) * (a1_ii + a2_iikk * f) # EQ 46 NOTE the typo from Stixrude 2005: eta_s = -gr - ( 1.0 / 2.0 * pow(nu_o_nu0_sq, -1.0) * pow((2.0 * f) + 1.0, 2.0) * a2_s ) return eta_s # calculate isotropic thermal pressure, see # Matas et. al. (2007) eq B4 def _thermal_pressure(self, T, V, params): Debye_T = self._debye_temperature(params["V_0"] / V, params) gr = self.grueneisen_parameter(0.0, T, V, params) # P not important P_th = gr * debye.thermal_energy(T, Debye_T, params["n"]) / V return P_th
[docs] def volume(self, pressure, temperature, params): """ Returns molar volume. :math:`[m^3]` """ T_0 = params["T_0"] Debye_0 = params["Debye_0"] V_0 = params["V_0"] dV = 1.0e-2 * params["V_0"] n = params["n"] a1_ii = 6.0 * params["grueneisen_0"] # EQ 47 a2_iikk = ( -12.0 * params["grueneisen_0"] + 36.0 * pow(params["grueneisen_0"], 2.0) - 18.0 * params["q_0"] * params["grueneisen_0"] ) # EQ 47 b_iikk = 9.0 * params["K_0"] # EQ 28 b_iikkmm = 27.0 * params["K_0"] * (params["Kprime_0"] - 4.0) # EQ 29z # Finding the volume at a given pressure requires a # root-finding scheme. Here we use brentq to find the root. # Root-finding using brentq requires bounds to be specified. # We do this using a bracketing function. args = ( pressure, temperature, V_0, T_0, Debye_0, n, a1_ii, a2_iikk, b_iikk, b_iikkmm, ) try: # The first attempt to find a bracket for # root finding uses V_0 as a starting point sol = bracket(_delta_pressure, V_0, dV, args) except Exception: # At high temperature, the naive bracketing above may # try a volume guess that exceeds the point at which the # bulk modulus goes negative at that temperature. # In this case, we try a more nuanced approach by # first finding the volume at which the bulk modulus goes # negative, and then either (a) raising an exception if the # desired pressure is less than the pressure at that volume, # or (b) using that pressure to create a better bracket for # brentq. def _K_T(V, T, params): return self.isothermal_bulk_modulus_reuss(0.0, T, V, params) sol_K_T = bracket(_K_T, V_0, dV, args=(temperature, params)) V_crit = opt.brentq( _K_T, sol_K_T[0], sol_K_T[1], args=(temperature, params) ) P_min = self.pressure(temperature, V_crit, params) if P_min > pressure: raise Exception( "The desired pressure is not achievable " "at this temperature. The minimum pressure " f"achievable is {P_min:.2e} Pa." ) else: try: sol = bracket(_delta_pressure, V_crit - dV, dV, args) except Exception: raise Exception( "Cannot find a volume, perhaps you are " "outside of the range of validity for " "the equation of state?" ) return opt.brentq(_delta_pressure, sol[0], sol[1], args=args)
[docs] def pressure(self, temperature, volume, params): """ Returns the pressure of the mineral at a given temperature and volume [Pa] """ debye_T = self._debye_temperature(params["V_0"] / volume, params) gr = self.grueneisen_parameter( 0.0, temperature, volume, params ) # does not depend on pressure # thermal energy at temperature T E_th = debye.thermal_energy(temperature, debye_T, params["n"]) # thermal energy at reference temperature E_th_ref = debye.thermal_energy(params["T_0"], debye_T, params["n"]) b_iikk = 9.0 * params["K_0"] # EQ 28 b_iikkmm = 27.0 * params["K_0"] * (params["Kprime_0"] - 4.0) # EQ 29 f = 0.5 * (pow(params["V_0"] / volume, 2.0 / 3.0) - 1.0) # EQ 24 P = (1.0 / 3.0) * (pow(1.0 + 2.0 * f, 5.0 / 2.0)) * ( (b_iikk * f) + (0.5 * b_iikkmm * pow(f, 2.0)) ) + gr * ( E_th - E_th_ref ) / volume # EQ 21 return P
[docs] def grueneisen_parameter(self, pressure, temperature, volume, params): """ Returns grueneisen parameter :math:`[unitless]` """ return _grueneisen_parameter_fast( params["V_0"], volume, params["grueneisen_0"], params["q_0"] )
[docs] def isothermal_bulk_modulus_reuss(self, pressure, temperature, volume, params): """ Returns isothermal bulk modulus :math:`[Pa]` """ T_0 = params["T_0"] debye_T = self._debye_temperature(params["V_0"] / volume, params) gr = self.grueneisen_parameter(pressure, temperature, volume, params) # thermal energy at temperature T E_th = debye.thermal_energy(temperature, debye_T, params["n"]) # thermal energy at reference temperature E_th_ref = debye.thermal_energy(T_0, debye_T, params["n"]) # heat capacity at temperature T C_v = debye.molar_heat_capacity_v(temperature, debye_T, params["n"]) # heat capacity at reference temperature C_v_ref = debye.molar_heat_capacity_v(T_0, debye_T, params["n"]) q = self.volume_dependent_q(params["V_0"] / volume, params) K = ( bm.bulk_modulus(volume, params) + (gr + 1.0 - q) * (gr / volume) * (E_th - E_th_ref) - (pow(gr, 2.0) / volume) * (C_v * temperature - C_v_ref * T_0) ) return K
[docs] def isentropic_bulk_modulus_reuss(self, pressure, temperature, volume, params): """ Returns adiabatic bulk modulus. :math:`[Pa]` """ K_T = self.isothermal_bulk_modulus_reuss(pressure, temperature, volume, params) alpha = self.thermal_expansivity(pressure, temperature, volume, params) gr = self.grueneisen_parameter(pressure, temperature, volume, params) K_S = K_T * (1.0 + gr * alpha * temperature) return K_S
[docs] def shear_modulus(self, pressure, temperature, volume, params): """ Returns shear modulus. :math:`[Pa]` """ T_0 = params["T_0"] debye_T = self._debye_temperature(params["V_0"] / volume, params) eta_s = self._isotropic_eta_s(params["V_0"] / volume, params) E_th = debye.thermal_energy(temperature, debye_T, params["n"]) E_th_ref = debye.thermal_energy(T_0, debye_T, params["n"]) if self.order == 2: return ( bm.shear_modulus_second_order(volume, params) - eta_s * (E_th - E_th_ref) / volume ) elif self.order == 3: return ( bm.shear_modulus_third_order(volume, params) - eta_s * (E_th - E_th_ref) / volume ) else: raise NotImplementedError("")
[docs] def molar_heat_capacity_v(self, pressure, temperature, volume, params): """ Returns heat capacity at constant volume. :math:`[J/K/mol]` """ debye_T = self._debye_temperature(params["V_0"] / volume, params) return debye.molar_heat_capacity_v(temperature, debye_T, params["n"])
[docs] def molar_heat_capacity_p(self, pressure, temperature, volume, params): """ Returns heat capacity at constant pressure. :math:`[J/K/mol]` """ alpha = self.thermal_expansivity(pressure, temperature, volume, params) gr = self.grueneisen_parameter(pressure, temperature, volume, params) C_v = self.molar_heat_capacity_v(pressure, temperature, volume, params) C_p = C_v * (1.0 + gr * alpha * temperature) return C_p
[docs] def thermal_expansivity(self, pressure, temperature, volume, params): """ Returns thermal expansivity. :math:`[1/K]` """ C_v = self.molar_heat_capacity_v(pressure, temperature, volume, params) gr = self.grueneisen_parameter(pressure, temperature, volume, params) K = self.isothermal_bulk_modulus_reuss(pressure, temperature, volume, params) alpha = gr * C_v / K / volume return alpha
[docs] def gibbs_free_energy(self, pressure, temperature, volume, params): """ Returns the Gibbs free energy at the pressure and temperature of the mineral [J/mol] """ G = ( self.helmholtz_free_energy(pressure, temperature, volume, params) + pressure * volume ) return G
[docs] def molar_internal_energy(self, pressure, temperature, volume, params): """ Returns the internal energy at the pressure and temperature of the mineral [J/mol] """ return self.helmholtz_free_energy( pressure, temperature, volume, params ) + temperature * self.entropy(pressure, temperature, volume, params)
[docs] def entropy(self, pressure, temperature, volume, params): """ Returns the entropy at the pressure and temperature of the mineral [J/K/mol] """ Debye_T = self._debye_temperature(params["V_0"] / volume, params) S = debye.entropy(temperature, Debye_T, params["n"]) return S
[docs] def enthalpy(self, pressure, temperature, volume, params): """ Returns the enthalpy at the pressure and temperature of the mineral [J/mol] """ return ( self.helmholtz_free_energy(pressure, temperature, volume, params) + temperature * self.entropy(pressure, temperature, volume, params) + pressure * volume )
[docs] def helmholtz_free_energy(self, pressure, temperature, volume, params): """ Returns the Helmholtz free energy at the pressure and temperature of the mineral [J/mol] """ x = params["V_0"] / volume f = 1.0 / 2.0 * (pow(x, 2.0 / 3.0) - 1.0) Debye_T = self._debye_temperature(params["V_0"] / volume, params) F_quasiharmonic = debye.helmholtz_free_energy( temperature, Debye_T, params["n"] ) - debye.helmholtz_free_energy(params["T_0"], Debye_T, params["n"]) b_iikk = 9.0 * params["K_0"] # EQ 28 b_iikkmm = 27.0 * params["K_0"] * (params["Kprime_0"] - 4.0) # EQ 29 F = ( params["F_0"] + 0.5 * b_iikk * f * f * params["V_0"] + (1.0 / 6.0) * params["V_0"] * b_iikkmm * f * f * f + F_quasiharmonic ) return F
[docs] def validate_parameters(self, params): """ Check for existence and validity of the parameters """ if "T_0" not in params: params["T_0"] = 300.0 # If eta_s_0 is not included this is presumably deliberate, # as we can model density and bulk modulus just fine without it, # so just add it to the dictionary as nan # The same goes for the standard state Helmholtz free energy if "eta_s_0" not in params: params["eta_s_0"] = float("nan") if "F_0" not in params: params["F_0"] = float("nan") # First, let's check the EoS parameters for Tref bm.BirchMurnaghanBase.validate_parameters(bm.BirchMurnaghanBase(), params) # Now check all the required keys for the # thermal part of the EoS are in the dictionary expected_keys = ["molar_mass", "n", "Debye_0", "grueneisen_0", "q_0", "eta_s_0"] 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["T_0"] < 0.0: warnings.warn("Unusual value for T_0", stacklevel=2) if params["molar_mass"] < 0.001 or params["molar_mass"] > 10.0: warnings.warn("Unusual value for molar_mass", stacklevel=2) if params["n"] < 1.0 or params["n"] > 1000.0: warnings.warn("Unusual value for n", stacklevel=2) if params["Debye_0"] < 1.0 or params["Debye_0"] > 10000.0: warnings.warn("Unusual value for Debye_0", stacklevel=2) if params["grueneisen_0"] < -1.0 or params["grueneisen_0"] > 10.0: warnings.warn("Unusual value for grueneisen_0", stacklevel=2) if params["q_0"] < -20.0 or params["q_0"] > 20.0: warnings.warn("Unusual value for q_0", stacklevel=2) if params["eta_s_0"] < -10.0 or params["eta_s_0"] > 10.0: warnings.warn("Unusual value for eta_s_0", stacklevel=2)
[docs] class SLB3(SLBBase): """ SLB equation of state with third order finite strain expansion for the shear modulus (this should be preferred, as it is more thermodynamically consistent.) """ def __init__(self): self.order = 3
[docs] class SLB2(SLBBase): """ SLB equation of state with second order finite strain expansion for the shear modulus. In general, this should not be used, but sometimes shear modulus data is fit to a second order equation of state. In that case, you should use this. The moral is, be careful! """ def __init__(self): self.order = 2