[cig-commits] [commit] split_mt: Alpha of MT split into isothermal and thermal (58450e1)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Dec 31 07:50:26 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : split_mt
Link : https://github.com/geodynamics/burnman/compare/0000000000000000000000000000000000000000...58450e1a99551f71cdd4dde6146b43ad18d66e5c
>---------------------------------------------------------------
commit 58450e1a99551f71cdd4dde6146b43ad18d66e5c
Author: Bob Myhill <myhill.bob at gmail.com>
Date: Wed Dec 31 15:49:49 2014 +0000
Alpha of MT split into isothermal and thermal
>---------------------------------------------------------------
58450e1a99551f71cdd4dde6146b43ad18d66e5c
burnman/eos/__init__.py | 1 +
burnman/eos/helper.py | 5 +-
burnman/eos/{modified_tait.py => hp.py} | 76 +++----
burnman/eos/modified_tait.py | 352 +++++++-------------------------
4 files changed, 104 insertions(+), 330 deletions(-)
diff --git a/burnman/eos/__init__.py b/burnman/eos/__init__.py
index 2009223..00b4cb2 100644
--- a/burnman/eos/__init__.py
+++ b/burnman/eos/__init__.py
@@ -14,6 +14,7 @@ from birch_murnaghan import BM2, BM3
from mie_grueneisen_debye import MGD2, MGD3
from slb import SLB2, SLB3
from modified_tait import MT
+from hp import HPMT
from cork import CORK
from helper import create
diff --git a/burnman/eos/helper.py b/burnman/eos/helper.py
index 54b6c4a..115c8f3 100644
--- a/burnman/eos/helper.py
+++ b/burnman/eos/helper.py
@@ -7,6 +7,7 @@ import slb
import mie_grueneisen_debye as mgd
import birch_murnaghan as bm
import modified_tait as mt
+import hp
import cork
from equation_of_state import EquationOfState
@@ -29,8 +30,10 @@ def create(method):
return bm.BM2()
elif method == "bm3":
return bm.BM3()
- elif method == "mtait":
+ elif method == "isomt":
return mt.MT()
+ elif method == "mtait":
+ return hp.HPMT()
elif method == "cork":
return cork.CORK()
else:
diff --git a/burnman/eos/modified_tait.py b/burnman/eos/hp.py
similarity index 88%
copy from burnman/eos/modified_tait.py
copy to burnman/eos/hp.py
index 752be96..afe2a79 100644
--- a/burnman/eos/modified_tait.py
+++ b/burnman/eos/hp.py
@@ -3,24 +3,43 @@
# Released under GPL v2 or later.
import numpy as np
+import warnings
+import modified_tait as mt
import equation_of_state as eos
+
import einstein
+
from burnman.endmemberdisorder import *
-import burnman.constants as constants
-import warnings
+
T_0=298.15 # Standard temperature = 25 C
P_0=1.e5 # Standard pressure = 1.e5 Pa
-class MT(eos.EquationOfState):
+class HPMT(eos.EquationOfState):
"""
Base class for a generic modified Tait equation of state.
References for this can be found in Huang and Chow (1974)
and Holland and Powell (2011; followed here).
"""
+ def volume(self, pressure,temperature,params):
+ """
+ Returns volume [m^3] as a function of pressure [Pa] and temperature [K]
+ EQ 12
+ """
+ Pth=self.__relative_thermal_pressure(temperature,params)
+ return mt.volume(pressure-Pth, params)
+
+ def pressure(self, temperature, volume, params):
+ """
+ Returns pressure [Pa] as a function of temperature [K] and volume[m^3]
+ EQ B7
+ """
+ Pth=self.__relative_thermal_pressure(temperature,params)
+ return mt.modified_tait(params['V_0']/volume, params) + Pth
+
def grueneisen_parameter(self, pressure, temperature, volume, params):
"""
Returns grueneisen parameter [unitless] as a function of pressure,
@@ -31,25 +50,13 @@ class MT(eos.EquationOfState):
C_V = self.heat_capacity_v( pressure, temperature, volume, params)
return alpha * K_T * volume / C_V
- def volume(self, pressure,temperature,params):
- """
- Returns volume [m^3] as a function of pressure [Pa] and temperature [K]
- EQ 12
- """
- a, b, c = self.__tait_constants(params)
- Pth=self.__relative_thermal_pressure(temperature,params)
- x = 1 - a*( 1. - np.power(( 1. + b*(pressure-Pth)), -1.0*c))
- return x*params['V_0']
-
def isothermal_bulk_modulus(self, pressure,temperature,volume, params):
"""
Returns isothermal bulk modulus [Pa] as a function of pressure [Pa],
temperature [K], and volume [m^3]. EQ 13+2
"""
- a, b, c = self.__tait_constants(params)
Pth=self.__relative_thermal_pressure(temperature,params)
- psubpth=pressure-Pth
- return params['K_0']*(1. + b*(psubpth))*(a + (1.-a)*np.power((1. + b*(psubpth)), c))
+ return mt.bulk_modulus(pressure-Pth, params)
#calculate the shear modulus as a function of P, V, and T
def shear_modulus(self, pressure, temperature, volume, params):
@@ -76,7 +83,7 @@ class MT(eos.EquationOfState):
Returns thermal expansivity at the pressure, temperature, and volume [1/K]
Replace -Pth in EQ 13+1 with P-Pth for non-ambient temperature
"""
- a, b, c = self.__tait_constants(params)
+ a, b, c = mt.tait_constants(params)
Pth=self.__relative_thermal_pressure(temperature,params)
psubpth=pressure-Pth
einstein_T=self.__einstein_temperature(params['S_0'], params['n'])
@@ -124,7 +131,7 @@ class MT(eos.EquationOfState):
and temperature [K].
"""
# Calculate temperature and pressure integrals
- a, b, c = self.__tait_constants(params)
+ a, b, c = mt.tait_constants(params)
Pth=self.__relative_thermal_pressure(temperature,params)
psubpth=pressure-Pth
@@ -149,7 +156,7 @@ class MT(eos.EquationOfState):
Returns the entropy [J/K/mol] as a function of pressure [Pa]
and temperature [K].
"""
- a, b, c = self.__tait_constants(params)
+ a, b, c = mt.tait_constants(params)
Pth=self.__relative_thermal_pressure(temperature,params)
einstein_T=self.__einstein_temperature(params['S_0'], params['n'])
@@ -192,7 +199,7 @@ class MT(eos.EquationOfState):
Returns the heat capacity [J/K/mol] as a function of pressure [Pa]
and temperature [K].
"""
- a, b, c = self.__tait_constants(params)
+ a, b, c = mt.tait_constants(params)
Pth=self.__relative_thermal_pressure(temperature,params)
einstein_T=self.__einstein_temperature(params['S_0'], params['n'])
@@ -208,15 +215,6 @@ class MT(eos.EquationOfState):
return self.heat_capacity_p0(temperature,params) + temperature*dSdT + Cpdisord
- # calculate P = P(T0) + Pth
- def pressure(self, temperature, volume, params):
- """
- Returns pressure [Pa] as a function of temperature [K] and volume[m^3]
- EQ B7
- """
- return self.__modified_tait(params['V_0']/volume, params) + \
- self.__relative_thermal_pressure(temperature, params)
-
def __einstein_temperature(self, S, n):
"""
@@ -225,27 +223,7 @@ class MT(eos.EquationOfState):
"""
return 10636./(S/n + 6.44)
- def __tait_constants(self, params):
- """
- returns parameters for the modified Tait equation of state
- derived from K_T and its two first pressure derivatives
- EQ 4 from Holland and Powell, 2011
- """
- a=(1.+params['Kprime_0'])/(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])
- b=params['Kprime_0']/params['K_0'] - params['Kdprime_0']/(1. + params['Kprime_0'])
- c=(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])/(params['Kprime_0']*params['Kprime_0'] + params['Kprime_0'] - params['K_0']*params['Kdprime_0'])
- return a, b, c
-
- def __modified_tait(self, x, params):
- """
- equation for the modified Tait equation of state, returns
- pressure in the same units that are supplied for the reference bulk
- modulus (params['K_0'])
- EQ 2 from Holland and Powell, 2011
- """
- a, b, c = self.__tait_constants(params)
- return (np.power((x + a - 1.) / a, -1./c) - 1.)/b
def __thermal_pressure(self,T,params):
"""
diff --git a/burnman/eos/modified_tait.py b/burnman/eos/modified_tait.py
index 752be96..00464f1 100644
--- a/burnman/eos/modified_tait.py
+++ b/burnman/eos/modified_tait.py
@@ -3,16 +3,49 @@
# Released under GPL v2 or later.
import numpy as np
-
import equation_of_state as eos
-import einstein
-from burnman.endmemberdisorder import *
-import burnman.constants as constants
import warnings
-T_0=298.15 # Standard temperature = 25 C
P_0=1.e5 # Standard pressure = 1.e5 Pa
+def tait_constants(params):
+ """
+ returns parameters for the modified Tait equation of state
+ derived from K_T and its two first pressure derivatives
+ EQ 4 from Holland and Powell, 2011
+ """
+ a=(1.+params['Kprime_0'])/(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])
+ b=params['Kprime_0']/params['K_0'] - params['Kdprime_0']/(1. + params['Kprime_0'])
+ c=(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])/(params['Kprime_0']*params['Kprime_0'] + params['Kprime_0'] - params['K_0']*params['Kdprime_0'])
+ return a, b, c
+
+def modified_tait(x, params):
+ """
+ equation for the modified Tait equation of state, returns
+ pressure in the same units that are supplied for the reference bulk
+ modulus (params['K_0'])
+ EQ 2 from Holland and Powell, 2011
+ """
+ a, b, c = tait_constants(params)
+ return (np.power((x + a - 1.) / a, -1./c) - 1.)/b
+
+def volume(pressure,params):
+ """
+ Returns volume [m^3] as a function of pressure [Pa] and temperature [K]
+ EQ 12
+ """
+ a, b, c = tait_constants(params)
+ x = 1 - a*( 1. - np.power(( 1. + b*(pressure)), -1.0*c))
+ return x*params['V_0']
+
+def bulk_modulus(pressure, params):
+ """
+ Returns isothermal bulk modulus :math:`K_T` of the mineral. :math:`[Pa]`.
+ EQ 13+2
+ """
+ a, b, c = tait_constants(params)
+ return params['K_0']*(1. + b*(pressure))*(a + (1.-a)*np.power((1. + b*(pressure)), c))
+
class MT(eos.EquationOfState):
"""
@@ -21,294 +54,71 @@ class MT(eos.EquationOfState):
and Holland and Powell (2011; followed here).
"""
- def grueneisen_parameter(self, pressure, temperature, volume, params):
- """
- Returns grueneisen parameter [unitless] as a function of pressure,
- temperature, and volume.
- """
- alpha = self.thermal_expansivity (pressure, temperature, volume, params)
- K_T = self.isothermal_bulk_modulus (pressure, temperature, volume, params)
- C_V = self.heat_capacity_v( pressure, temperature, volume, params)
- return alpha * K_T * volume / C_V
-
def volume(self, pressure,temperature,params):
"""
- Returns volume [m^3] as a function of pressure [Pa] and temperature [K]
- EQ 12
- """
- a, b, c = self.__tait_constants(params)
- Pth=self.__relative_thermal_pressure(temperature,params)
- x = 1 - a*( 1. - np.power(( 1. + b*(pressure-Pth)), -1.0*c))
- return x*params['V_0']
-
- def isothermal_bulk_modulus(self, pressure,temperature,volume, params):
- """
- Returns isothermal bulk modulus [Pa] as a function of pressure [Pa],
- temperature [K], and volume [m^3]. EQ 13+2
- """
- a, b, c = self.__tait_constants(params)
- Pth=self.__relative_thermal_pressure(temperature,params)
- psubpth=pressure-Pth
- return params['K_0']*(1. + b*(psubpth))*(a + (1.-a)*np.power((1. + b*(psubpth)), c))
-
- #calculate the shear modulus as a function of P, V, and T
- def shear_modulus(self, pressure, temperature, volume, params):
+ Returns volume :math:`[m^3]` as a function of pressure :math:`[Pa]`.
"""
- Not implemented.
- Returns 0.
- Could potentially apply a fixed Poissons ratio as a rough estimate.
- """
- return 0.
-
- # Cv, heat capacity at constant volume
- def heat_capacity_v(self, pressure, temperature, volume, params):
- """
- Returns heat capacity at constant volume at the pressure, temperature, and volume [J/K/mol].
- """
- C_p=self.heat_capacity_p(pressure, temperature, volume, params)
- V=self.volume(pressure,temperature,params)
- alpha=self.thermal_expansivity(pressure, temperature, volume , params)
- K_T=self.isothermal_bulk_modulus(pressure,temperature,volume, params)
- return C_p - V*temperature*alpha*alpha*K_T
-
- def thermal_expansivity(self, pressure, temperature, volume , params):
- """
- Returns thermal expansivity at the pressure, temperature, and volume [1/K]
- Replace -Pth in EQ 13+1 with P-Pth for non-ambient temperature
- """
- a, b, c = self.__tait_constants(params)
- Pth=self.__relative_thermal_pressure(temperature,params)
- psubpth=pressure-Pth
- einstein_T=self.__einstein_temperature(params['S_0'], params['n'])
- C_V0 = einstein.heat_capacity_v( T_0, einstein_T, params['n'] )
- C_V = einstein.heat_capacity_v(temperature, einstein_T,params['n'])
- alpha = params['a_0'] * (C_V/C_V0) *1./((1.+b*psubpth)*(a + (1.-a)*np.power((1+b*psubpth), c)))
-
- return alpha
-
- def heat_capacity_p0(self,temperature,params):
- """
- Returns heat capacity at ambient pressure as a function of temperature [J/K/mol]
- Cp = a + bT + cT^-2 + dT^-0.5 in Holland and Powell, 2011
- """
- Cp = params['Cp'][0] + params['Cp'][1]*temperature + params['Cp'][2]*np.power(temperature,-2.) + params['Cp'][3]*np.power(temperature,-0.5)
- return Cp
-
- def heat_capacity_p_einstein(self, pressure, temperature, volume, params):
- """
- Returns heat capacity at constant pressure at the pressure, temperature, and volume, using the C_v and Einstein model [J/K/mol]
- WARNING: Only for comparison with internally self-consistent C_p
- """
- alpha = self.thermal_expansivity(pressure, temperature, volume, params)
- gr = self.grueneisen_parameter(pressure, temperature, volume, params)
- C_v = self.heat_capacity_v(pressure, temperature, volume, params)
- C_p = C_v*(1. + gr * alpha * temperature)
- return C_p
-
-
- def adiabatic_bulk_modulus(self,pressure,temperature,volume,params):
- """
- Returns adiabatic bulk modulus [Pa] as a function of pressure [Pa],
- temperature [K], and volume [m^3].
- """
- K_T= self.isothermal_bulk_modulus(pressure,temperature,volume,params)
- alpha = self.thermal_expansivity(pressure,temperature,volume,params)
- C_p = self.heat_capacity_p(pressure, temperature, volume, params)
- C_v = self.heat_capacity_v(pressure, temperature, volume, params)
- K_S = K_T*C_p/C_v
- return K_S
-
- def gibbs_free_energy(self,pressure,temperature, volume, params):
- """
- Returns the gibbs free energy [J/mol] as a function of pressure [Pa]
- and temperature [K].
- """
- # Calculate temperature and pressure integrals
- a, b, c = self.__tait_constants(params)
- Pth=self.__relative_thermal_pressure(temperature,params)
-
- psubpth=pressure-Pth
+ return volume(pressure,params)
- # EQ 13
- intVdP = pressure*params['V_0']*(1. - a + (a*(np.power((1.-b*Pth), 1.-c) - np.power((1. + b*(pressure-Pth)), 1.-c))/(b*(c-1.)*pressure)))
-
- # Add order-disorder terms if required
- if params.has_key('landau_Tc'): # For a phase transition described by Landau term
- Gdisord=gibbs_disorder_Landau(pressure, temperature, params)
- else:
- if params.has_key('BW_deltaH'): # Add Bragg-Williams disordering
- Gdisord=gibbs_disorder_BW(pressure, temperature, params) - gibbs_disorder_BW(P_0, T_0, params)
- else:
- Gdisord=0.0
-
- return params['H_0'] + self.__intCpdT(temperature, params) - temperature*(params['S_0'] + self.__intCpoverTdT(temperature, params)) + intVdP + Gdisord
-
-
- def entropy(self,pressure,temperature, volume, params):
- """
- Returns the entropy [J/K/mol] as a function of pressure [Pa]
- and temperature [K].
- """
- a, b, c = self.__tait_constants(params)
- Pth=self.__relative_thermal_pressure(temperature,params)
-
- einstein_T=self.__einstein_temperature(params['S_0'], params['n'])
- ksi_over_ksi_0=einstein.heat_capacity_v( temperature, einstein_T, params['n'] )/einstein.heat_capacity_v( T_0, einstein_T, params['n'] )
-
- dintVdpdx=(params['V_0']*params['a_0']*params['K_0']*a*ksi_over_ksi_0)*(np.power((1.+b*(pressure-Pth)), 0.-c) - np.power((1.-b*Pth), 0.-c))
-
- # Add order-disorder terms if required
- if params.has_key('landau_Tc'): # For a phase transition described by Landau term
- Sdisord=entropy_disorder_Landau(pressure, temperature, params)
- else:
- if params.has_key('BW_deltaH'): # Add Bragg-Williams disordering
- Sdisord=entropy_disorder_BW(pressure, temperature, params) - entropy_disorder_BW(P_0, T_0, params)
- else:
- Sdisord=0.0
-
- return params['S_0'] + self.__intCpoverTdT(temperature, params) + dintVdpdx + Sdisord
-
- def enthalpy(self, pressure, temperature, volume, params):
- """
- Returns the enthalpy [J/mol] as a function of pressure [Pa]
- and temperature [K].
- """
- gibbs=self.gibbs_free_energy(pressure,temperature,volume, params)
- entropy=self.entropy(pressure,temperature,volume, params)
-
- # Add order-disorder terms if required
- if params.has_key('landau_Tc'): # For a phase transition described by Landau term
- Hdisord=enthalpy_disorder_Landau(pressure, temperature, params)
- else:
- if params.has_key('BW_deltaH'): # Add Bragg-Williams disordering
- Hdisord=enthalpy_disorder_BW(pressure, temperature, params) - enthalpy_disorder_BW(P_0, T_0, params)
- else:
- Hdisord=0.0
-
- return gibbs + temperature*entropy + Hdisord
-
- def heat_capacity_p(self, pressure, temperature, volume, params):
- """
- Returns the heat capacity [J/K/mol] as a function of pressure [Pa]
- and temperature [K].
- """
- a, b, c = self.__tait_constants(params)
- Pth=self.__relative_thermal_pressure(temperature,params)
-
- einstein_T=self.__einstein_temperature(params['S_0'], params['n'])
- ksi_over_ksi_0=einstein.heat_capacity_v( temperature, einstein_T, params['n'] )/einstein.heat_capacity_v( T_0, einstein_T, params['n'] )
-
- dSdT=params['V_0']*params['K_0']*np.power((ksi_over_ksi_0*params['a_0']),2.0)*(np.power((1.+b*(pressure-Pth)), -1.-c) - np.power((1.-b*Pth), -1.-c))
-
- # Add order-disorder terms if required
- if params.has_key('landau_Tc'): # For a phase transition described by Landau term
- Cpdisord=heat_capacity_p_disorder_Landau(pressure, temperature, params)
- else:
- Cpdisord=0.0
-
- return self.heat_capacity_p0(temperature,params) + temperature*dSdT + Cpdisord
-
- # calculate P = P(T0) + Pth
def pressure(self, temperature, volume, params):
"""
Returns pressure [Pa] as a function of temperature [K] and volume[m^3]
- EQ B7
- """
- return self.__modified_tait(params['V_0']/volume, params) + \
- self.__relative_thermal_pressure(temperature, params)
-
-
- def __einstein_temperature(self, S, n):
- """
- Empirical Einstein temperature
- Holland and Powell, 2011; base of p.346, para.1
"""
- return 10636./(S/n + 6.44)
+ return modified_tait(params['V_0']/volume, params)
- def __tait_constants(self, params):
+ def isothermal_bulk_modulus(self, pressure,temperature,volume, params):
"""
- returns parameters for the modified Tait equation of state
- derived from K_T and its two first pressure derivatives
- EQ 4 from Holland and Powell, 2011
+ Returns isothermal bulk modulus :math:`K_T` of the mineral. :math:`[Pa]`.
"""
- a=(1.+params['Kprime_0'])/(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])
- b=params['Kprime_0']/params['K_0'] - params['Kdprime_0']/(1. + params['Kprime_0'])
- c=(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])/(params['Kprime_0']*params['Kprime_0'] + params['Kprime_0'] - params['K_0']*params['Kdprime_0'])
- return a, b, c
+ return bulk_modulus(pressure, params)
- def __modified_tait(self, x, params):
+ def adiabatic_bulk_modulus(self,pressure, temperature, volume, params):
"""
- equation for the modified Tait equation of state, returns
- pressure in the same units that are supplied for the reference bulk
- modulus (params['K_0'])
- EQ 2 from Holland and Powell, 2011
+ Since this equation of state does not contain temperature effects, simply return a very large number. :math:`[Pa]`
"""
- a, b, c = self.__tait_constants(params)
-
- return (np.power((x + a - 1.) / a, -1./c) - 1.)/b
+ return 1.e99
- def __thermal_pressure(self,T,params):
+ def shear_modulus(self, pressure, temperature, volume, params):
"""
- Returns thermal pressure [Pa] as a function of T [K]
- EQ 12 - 1 of Holland and Powell, 2011
+ Not implemented in the Modified Tait EoS. :math:`[Pa]`
+ Returns 0.
+ Could potentially apply a fixed Poissons ratio as a rough estimate.
"""
+ return 0.
- # This is basically the mie-gruneisen equation of state for thermal
- # pressure using an Einstein model for heat capacity. The additional
- # assumption that they make is that alpha*K/Cv, (or gamma / V) is
- # constant over a wide range of compressions.
-
- # Note that the xi function in HP2011 is just the Einstein heat capacity
- # divided by 3nR. I don't know why they don't use that, but anyhow...
-
- einstein_T=self.__einstein_temperature(params['S_0'],params['n'])
- E_th = einstein.thermal_energy( T, einstein_T, params['n'] )
- C_V0 = einstein.heat_capacity_v( T_0, einstein_T, params['n'] )
- P_th = params['a_0']*params['K_0'] / C_V0 * E_th
- return P_th
-
- def __relative_thermal_pressure( self, T, params):
+ def heat_capacity_v(self, pressure, temperature, volume, params):
"""
- Returns relative thermal pressure [Pa] as a function of T-T_0 [K]
- EQ 12 - 1 of Holland and Powell, 2011
+ Since this equation of state does not contain temperature effects, simply return a very large number. :math:`[J/K/mol]`
"""
- return self.__thermal_pressure(T, params) - \
- self.__thermal_pressure(T_0, params)
+ return 1.e99
- def __intCpdT (self, temperature, params):
+ def thermal_expansivity(self,pressure, temperature, volume, params):
"""
- Returns the thermal addition to the standard state enthalpy [J/mol]
- at ambient pressure [Pa]
+ Since this equation of state does not contain temperature effects, simply return zero. :math:`[1/K]`
"""
- return (params['Cp'][0]*temperature + 0.5*params['Cp'][1]*np.power(temperature,2.) - params['Cp'][2]/temperature + 2.*params['Cp'][3]*np.sqrt(temperature)) - (params['Cp'][0]*T_0 + 0.5*params['Cp'][1]*T_0*T_0 - params['Cp'][2]/T_0 + 2.0*params['Cp'][3]*np.sqrt(T_0))
+ return 0.
- def __intCpoverTdT (self, temperature, params):
+ def grueneisen_parameter(self,pressure,temperature,volume,params):
"""
- Returns the thermal addition to the standard state entropy [J/K/mol]
- at ambient pressure [Pa]
+ Since this equation of state does not contain temperature effects, simply return zero. :math:`[unitless]`
"""
- return (params['Cp'][0]*np.log(temperature) + params['Cp'][1]*temperature - 0.5*params['Cp'][2]/np.power(temperature,2.) - 2.0*params['Cp'][3]/np.sqrt(temperature)) - (params['Cp'][0]*np.log(T_0) + params['Cp'][1]*T_0 - 0.5*params['Cp'][2]/(T_0*T_0) - 2.0*params['Cp'][3]/np.sqrt(T_0))
-
+ return 0.
def validate_parameters(self, params):
"""
Check for existence and validity of the parameters
"""
- #if G and Gprime are not included this is presumably deliberate,
- #as we can model density and bulk modulus just fine without them,
- #so just add them to the dictionary as nans
- if 'H_0' not in params:
- params['H_0'] = float('nan')
- if 'S_0' not in params:
- params['S_0'] = float('nan')
+ # G and Gprime are not defined in this equation of state,
+ # We can model density and bulk modulus just fine without them,
+ # so just add them to the dictionary as nans
if 'G_0' not in params:
params['G_0'] = float('nan')
if 'Gprime_0' not in params:
params['Gprime_0'] = float('nan')
#check that all the required keys are in the dictionary
- expected_keys = ['H_0', 'S_0', 'V_0', 'Cp', 'a_0', 'K_0', 'Kprime_0', 'Kdprime_0', 'n', 'molar_mass']
+ expected_keys = ['V_0', 'K_0', 'Kprime_0', 'Kdprime_0', 'G_0', 'Gprime_0']
for k in expected_keys:
if k not in params:
raise KeyError('params object missing parameter : ' + k)
@@ -316,32 +126,14 @@ class MT(eos.EquationOfState):
#now check that the values are reasonable. I mostly just
#made up these values from experience, and we are only
#raising a warning. Better way to do this? [IR]
- if params['G_0'] is not float('nan') and (params['G_0'] < 0. or params['G_0'] > 1.e13):
- warnings.warn( 'Unusual value for G_0', stacklevel=2 )
- if params['Gprime_0'] is not float('nan') and (params['Gprime_0'] < -5. or params['Gprime_0'] > 10.):
- warnings.warn( 'Unusual value for Gprime_0', stacklevel=2 )
-
- # no test for H_0
- if params['S_0'] is not float('nan') and params['S_0'] < 0.:
- warnings.warn( 'Unusual value for S_0', stacklevel=2 )
- if params['V_0'] < 1.e-7 or params['V_0'] > 1.e-2:
+ if params['V_0'] < 1.e-7 or params['V_0'] > 1.e-3:
warnings.warn( 'Unusual value for V_0', stacklevel=2 )
-
-
- if self.heat_capacity_p0(T_0,params) < 0.:
- warnings.warn( 'Negative heat capacity at T_0', stacklevel=2 )
- if self.heat_capacity_p0(2000.,params) < 0.:
- warnings.warn( 'Negative heat capacity at 2000K', stacklevel=2 )
-
- if params['a_0'] < 0. or params['a_0'] > 1.e-3:
- warnings.warn( 'Unusual value for a_0', stacklevel=2 )
if params['K_0'] < 1.e9 or params['K_0'] > 1.e13:
- warnings.warn( 'Unusual value for K_0', stacklevel=2 )
- if params['Kprime_0'] < 0. or params['Kprime_0'] > 10.:
+ warnings.warn( 'Unusual value for K_0' , stacklevel=2)
+ if params['Kprime_0'] < -5. or params['Kprime_0'] > 10.:
warnings.warn( 'Unusual value for Kprime_0', stacklevel=2 )
- # no test for Kdprime_0
+ if params['G_0'] < 0.0 or params['G_0'] > 1.e13:
+ warnings.warn( 'Unusual value for G_0', stacklevel=2 )
+ if params['Gprime_0'] < -5. or params['Gprime_0'] > 10.:
+ warnings.warn( 'Unusual value for Gprime_0', stacklevel=2 )
- if params['n'] < 1. or params['n'] > 1000.:
- warnings.warn( 'Unusual value for n', stacklevel=2 )
- if params['molar_mass'] < 0.001 or params['molar_mass'] > 10.:
- warnings.warn( 'Unusual value for molar_mass', stacklevel=2 )
More information about the CIG-COMMITS
mailing list