[cig-commits] [commit] jit_optimizations: eos restructuring for performance reasons (ed9edf1)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Sun Dec 14 18:07:32 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : jit_optimizations
Link : https://github.com/geodynamics/burnman/compare/841a4b72425b7731caeeeb7f0a992a8b9b97230a...ed9edf1a9e99173d286d95de225645f19cce153f
>---------------------------------------------------------------
commit ed9edf1a9e99173d286d95de225645f19cce153f
Author: Timo Heister <timo.heister at gmail.com>
Date: Sat Dec 13 14:42:59 2014 -0800
eos restructuring for performance reasons
Only SLB is working for now.
>---------------------------------------------------------------
ed9edf1a9e99173d286d95de225645f19cce153f
burnman/eos/equation_of_state.py | 26 +++---
burnman/eos/slb.py | 180 ++++++++++++++++++++++-----------------
burnman/mineral.py | 36 ++++++--
3 files changed, 146 insertions(+), 96 deletions(-)
diff --git a/burnman/eos/equation_of_state.py b/burnman/eos/equation_of_state.py
index 82f4a5a..34b6300 100644
--- a/burnman/eos/equation_of_state.py
+++ b/burnman/eos/equation_of_state.py
@@ -80,7 +80,7 @@ class EquationOfState(object):
"""
return params["molar_mass"] / self.volume(pressure, temperature, params)
- def grueneisen_parameter(self, pressure, temperature, volume, params):
+ def grueneisen_parameter(self, mineral):
"""
Parameters
----------
@@ -101,7 +101,7 @@ class EquationOfState(object):
"""
raise NotImplementedError("")
- def isothermal_bulk_modulus(self, pressure, temperature, volume, params):
+ def isothermal_bulk_modulus(self, mineral):
"""
Parameters
----------
@@ -122,7 +122,7 @@ class EquationOfState(object):
"""
raise NotImplementedError("")
- def adiabatic_bulk_modulus(self, pressure, temperature, volume, params):
+ def adiabatic_bulk_modulus(self, mineral):
"""
Parameters
----------
@@ -143,7 +143,7 @@ class EquationOfState(object):
"""
raise NotImplementedError("")
- def shear_modulus(self, pressure, temperature, volume, params):
+ def shear_modulus(self, mineral):
"""
Parameters
----------
@@ -164,7 +164,7 @@ class EquationOfState(object):
"""
raise NotImplementedError("")
- def heat_capacity_v(self, pressure, temperature, volume, params):
+ def heat_capacity_v(self, mineral):
"""
Parameters
----------
@@ -185,7 +185,7 @@ class EquationOfState(object):
"""
raise NotImplementedError("")
- def heat_capacity_p(self, pressure, temperature, volume, params):
+ def heat_capacity_p(self, mineral):
"""
Parameters
----------
@@ -206,7 +206,7 @@ class EquationOfState(object):
"""
raise NotImplementedError("")
- def thermal_expansivity(self, pressure, temperature, volume, params):
+ def thermal_expansivity(self, mineral):
"""
Parameters
----------
@@ -227,7 +227,7 @@ class EquationOfState(object):
"""
raise NotImplementedError("")
- def reference_temperature( self, params ):
+ def reference_temperature( self, params):
"""
Parameters
----------
@@ -245,7 +245,7 @@ class EquationOfState(object):
else:
return 300.0
- def gibbs_free_energy( self, pressure, temperature, volume, params ):
+ def gibbs_free_energy(self, mineral):
"""
Parameters
----------
@@ -266,7 +266,7 @@ class EquationOfState(object):
"""
raise NotImplementedError("")
- def helmholtz_free_energy( self, pressure, temperature, volume, params ):
+ def helmholtz_free_energy(self, mineral):
"""
Parameters
----------
@@ -285,14 +285,14 @@ class EquationOfState(object):
"""
raise NotImplementedError("")
- def entropy( self, pressure, temperature, volume, params):
+ def entropy( self, mineral):
"""
Returns the entropy at the pressure and temperature of the mineral [J/K/mol]
"""
raise NotImplementedError("")
- def enthalpy( self, pressure, temperature, volume, params ):
+ def enthalpy( self, mineral):
"""
Parameters
----------
@@ -310,7 +310,7 @@ class EquationOfState(object):
"""
raise NotImplementedError("")
- def internal_energy( self, pressure, temperature, volume, params ):
+ def internal_energy( self, mineral):
"""
Parameters
----------
diff --git a/burnman/eos/slb.py b/burnman/eos/slb.py
index cabfa06..f1c6521 100644
--- a/burnman/eos/slb.py
+++ b/burnman/eos/slb.py
@@ -66,32 +66,50 @@ class SLBBase(eos.EquationOfState):
a2_iikk = -12.*params['grueneisen_0']+36.*pow(params['grueneisen_0'],2.) - 18.*params['q_0']*params['grueneisen_0'] # EQ 47
return params['Debye_0'] * np.sqrt(1. + a1_ii * f + 1./2. * a2_iikk*f*f)
- def volume_dependent_q(self, x, params):
+ def _volume_dependent_q(self, mineral):
"""
Finite strain approximation for :math:`q`, the isotropic volume strain
derivative of the grueneisen parameter.
"""
- f = 1./2. * (pow(x, 2./3.) - 1.)
- a1_ii = 6. * params['grueneisen_0'] # EQ 47
- a2_iikk = -12.*params['grueneisen_0']+36.*pow(params['grueneisen_0'],2.) - 18.*params['q_0']*params['grueneisen_0'] # EQ 47
- nu_o_nu0_sq = 1.+ a1_ii*f + (1./2.)*a2_iikk * f*f # EQ 41
- gr = 1./6./nu_o_nu0_sq * (2.*f+1.) * ( a1_ii + a2_iikk*f )
- q = 1./9.*(18.*gr - 6. - 1./2. / nu_o_nu0_sq * (2.*f+1.)*(2.*f+1.)*a2_iikk/gr)
- return q
-
- def __isotropic_eta_s(self, x, params):
+ if 'q' not in mineral.eos_data:
+ x = mineral.params['V_0'] / mineral.molar_volume()
+ grueneisen_0 = mineral.params['grueneisen_0']
+ q_0 = mineral.params['q_0']
+ f = 1./2. * (pow(x, 2./3.) - 1.)
+ a1_ii = 6. * grueneisen_0 # EQ 47
+ a2_iikk = -12.*grueneisen_0+36.*pow(grueneisen_0, 2.) - 18.*q_0*grueneisen_0 # EQ 47
+ nu_o_nu0_sq = 1.+ a1_ii*f + (1./2.)*a2_iikk * f*f # EQ 41
+ gr = mineral.grueneisen_parameter()
+ mineral.eos_data['q'] = 1./9.*(18.*gr - 6. - 1./2. / nu_o_nu0_sq * (2.*f+1.)*(2.*f+1.)*a2_iikk/gr)
+ return mineral.eos_data['q']
+
+
+ def _debye_T(self, mineral):
+ if 'debye_T' not in mineral.eos_data:
+ grueneisen_0 = mineral.params['grueneisen_0']
+ q_0 = mineral.params['q_0']
+ Debye_0 = mineral.params['Debye_0']
+ mineral.eos_data['debye_T'] = _debye_temperature_fast(mineral.params['V_0']/mineral.molar_volume(),
+ grueneisen_0, Debye_0, q_0)
+ return mineral.eos_data['debye_T']
+
+ def _isotropic_eta_s(self, mineral):
"""
Finite strain approximation for :math:`eta_{s0}`, the isotropic shear
strain derivative of the grueneisen parameter.
"""
- f = 1./2. * (pow(x, 2./3.) - 1.)
- a2_s = -2.*params['grueneisen_0'] - 2.*params['eta_s_0'] # EQ 47
- a1_ii = 6. * params['grueneisen_0'] # EQ 47
- a2_iikk = -12.*params['grueneisen_0']+36.*pow(params['grueneisen_0'],2.) - 18.*params['q_0']*params['grueneisen_0'] # EQ 47
- nu_o_nu0_sq = 1.+ a1_ii*f + (1./2.)*a2_iikk * pow(f,2.) # EQ 41
- gr = 1./6./nu_o_nu0_sq * (2.*f+1.) * ( a1_ii + a2_iikk*f )
- eta_s = - gr - (1./2. * pow(nu_o_nu0_sq,-1.) * pow((2.*f)+1.,2.)*a2_s) # EQ 46 NOTE the typo from Stixrude 2005
- return eta_s
+ if 'eta_s' not in mineral.eos_data:
+ x = mineral.params['V_0']/mineral.molar_volume()
+ f = 1./2. * (pow(x, 2./3.) - 1.)
+ grueneisen_0 = mineral.params['grueneisen_0']
+ q_0 = mineral.params['q_0']
+ a2_s = -2.*grueneisen_0 - 2.*mineral.params['eta_s_0'] # EQ 47
+ a1_ii = 6. * grueneisen_0 # EQ 47
+ a2_iikk = -12.*grueneisen_0+36.*pow(grueneisen_0, 2.) - 18.*q_0*grueneisen_0 # EQ 47
+ nu_o_nu0_sq = 1.+ a1_ii*f + (1./2.)*a2_iikk * pow(f,2.) # EQ 41
+ gr = mineral.grueneisen_parameter()
+ mineral.eos_data['eta_s'] = - gr - (1./2. * pow(nu_o_nu0_sq, -1.) * pow((2.*f)+1., 2.)*a2_s) # EQ 46 NOTE the typo from Stixrude 2005
+ return mineral.eos_data['eta_s']
def pressure(self, temperature, volume, params):
return bm.birch_murnaghan(params['V_0']/volume, params) + \
@@ -125,10 +143,10 @@ class SLBBase(eos.EquationOfState):
Returns molar volume. :math:`[m^3]`
"""
T_0 = self.reference_temperature( params )
- debye_T = lambda x : self.__debye_temperature(params['V_0']/x, params)
- gr = lambda x : self.grueneisen_parameter(pressure, temperature, x, params)
- E_th = lambda x : debye.thermal_energy(temperature, debye_T(x), params['n']) #thermal energy at temperature T
- E_th_ref = lambda x : debye.thermal_energy(T_0, debye_T(x), params['n']) #thermal energy at reference temperature
+ #debye_T = lambda x : self.__debye_temperature(params['V_0']/x, params)
+ #gr = lambda x : self.grueneisen_parameter(pressure, temperature, x, params)
+ #E_th = lambda x : debye.thermal_energy(temperature, debye_T(x), params['n']) #thermal energy at temperature T
+ #E_th_ref = lambda x : debye.thermal_energy(T_0, debye_T(x), params['n']) #thermal energy at reference temperature
b_iikk= 9.*params['K_0'] # EQ 28
b_iikkmm= 27.*params['K_0']*(params['Kprime_0']-4.) # EQ 29
@@ -179,24 +197,28 @@ class SLBBase(eos.EquationOfState):
return P
+ def grueneisen_parameter(self, mineral):
+ gruen_0 = mineral.params['grueneisen_0']
+ V_0 = mineral.params['V_0']
+ q_0 = mineral.params['q_0']
+ return _grueneisen_parameter_fast(V_0, mineral.V, gruen_0, q_0)
- def grueneisen_parameter(self, pressure, temperature, volume, params):
- """
- Returns grueneisen parameter :math:`[unitless]`
- """
- gruen_0 = params['grueneisen_0']
- V_0 = params['V_0']
- q_0 = params['q_0']
- return _grueneisen_parameter_fast(V_0, volume, gruen_0, q_0)
-
- def isothermal_bulk_modulus(self, pressure,temperature, volume, params):
+ def isothermal_bulk_modulus(self, mineral):
"""
Returns isothermal bulk modulus :math:`[Pa]`
"""
+ temperature = mineral.temperature
+ params = mineral.params
+ volume = mineral.molar_volume()
T_0 = self.reference_temperature( params )
- debye_T = self.__debye_temperature(params['V_0']/volume, params)
- gr = self.grueneisen_parameter(pressure, temperature, volume, params)
+ gruen_0 = mineral.params['grueneisen_0']
+ V_0 = mineral.params['V_0']
+ q_0 = mineral.params['q_0']
+ Debye_0 = params['Debye_0']
+
+ debye_T = self._debye_T(mineral)
+ gr = mineral.grueneisen_parameter()
E_th = debye.thermal_energy(temperature, debye_T, params['n']) #thermal energy at temperature T
E_th_ref = debye.thermal_energy(T_0,debye_T, params['n']) #thermal energy at reference temperature
@@ -204,7 +226,7 @@ class SLBBase(eos.EquationOfState):
C_v = debye.heat_capacity_v(temperature, debye_T, params['n']) #heat capacity at temperature T
C_v_ref = debye.heat_capacity_v(T_0,debye_T, params['n']) #heat capacity at reference temperature
- q = self.volume_dependent_q(params['V_0']/volume, params)
+ q = self._volume_dependent_q(mineral)
K = bm.bulk_modulus(volume, params) \
+ (gr + 1.-q)* ( gr / volume ) * (E_th - E_th_ref) \
@@ -212,26 +234,29 @@ class SLBBase(eos.EquationOfState):
return K
- def adiabatic_bulk_modulus(self, pressure, temperature, volume, params):
+ def adiabatic_bulk_modulus(self, mineral):
"""
Returns adiabatic bulk modulus. :math:`[Pa]`
"""
- K_T=self.isothermal_bulk_modulus(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. + gr * alpha * temperature)
+ K_T = mineral.isothermal_bulk_modulus()
+ alpha = mineral.thermal_expansivity()
+ gr = mineral.grueneisen_parameter()
+ K_S = K_T*(1. + gr * alpha * mineral.temperature)
return K_S
- def shear_modulus(self, pressure, temperature, volume, params):
+ def shear_modulus(self, mineral):
"""
Returns shear modulus. :math:`[Pa]`
"""
+ params = mineral.params
+ volume = mineral.molar_volume()
T_0 = self.reference_temperature( params )
- debye_T = self.__debye_temperature(params['V_0']/volume, params)
- eta_s = self.__isotropic_eta_s(params['V_0']/volume, params)
+ debye_T = self._debye_T(mineral)
+
+ eta_s = self._isotropic_eta_s(mineral) # TODO
- E_th = debye.thermal_energy(temperature ,debye_T, params['n'])
- E_th_ref = debye.thermal_energy(T_0,debye_T, params['n'])
+ E_th = debye.thermal_energy(mineral.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
@@ -240,74 +265,73 @@ class SLBBase(eos.EquationOfState):
else:
raise NotImplementedError("")
- def heat_capacity_v(self, pressure, temperature, volume, params):
+ def heat_capacity_v(self, mineral):
"""
Returns heat capacity at constant volume. :math:`[J/K/mol]`
"""
- debye_T = self.__debye_temperature(params['V_0']/volume, params)
- return debye.heat_capacity_v(temperature, debye_T,params['n'])
+ debye_T = self._debye_T(mineral)
+ return debye.heat_capacity_v(mineral.temperature, debye_T, mineral.params['n'])
- def heat_capacity_p(self, pressure, temperature, volume, params):
+ def heat_capacity_p(self, mineral):
"""
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.heat_capacity_v(pressure, temperature, volume, params)
- C_p = C_v*(1. + gr * alpha * temperature)
+ alpha = mineral.thermal_expansivity()
+ gr = mineral.grueneisen_parameter()
+ C_v = mineral.heat_capacity_v()
+ C_p = C_v*(1. + gr * alpha * mineral.temperature)
return C_p
- def thermal_expansivity(self, pressure, temperature, volume, params):
+ def thermal_expansivity(self, mineral):
"""
Returns thermal expansivity. :math:`[1/K]`
"""
- C_v = self.heat_capacity_v(pressure, temperature, volume, params)
- gr = self.grueneisen_parameter(pressure, temperature, volume, params)
- K = self.isothermal_bulk_modulus(pressure, temperature, volume, params)
- alpha = gr * C_v / K / volume
+ C_v = mineral.heat_capacity_v()
+ gr = mineral.grueneisen_parameter()
+ K = mineral.isothermal_bulk_modulus()
+ alpha = gr * C_v / K / mineral.molar_volume()
return alpha
- def gibbs_free_energy( self, pressure, temperature, volume, params):
+ def gibbs_free_energy(self, mineral):
"""
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
+ F = mineral.molar_helmholtz(mineral)
+ G = F + mineral.pressure * mineral.molar_volume()
return G
- def entropy( self, pressure, temperature, volume, params):
+ def entropy( self, mineral):
"""
Returns the entropy at the pressure and temperature of the mineral [J/K/mol]
"""
- x = params['V_0'] / volume
- f = 1./2. * (pow(x, 2./3.) - 1.)
- Debye_T = self.__debye_temperature(params['V_0']/volume, params)
- S = debye.entropy( temperature, Debye_T, params['n'] )
+ Debye_T = self._debye_T(mineral)
+ S = debye.entropy(mineral.temperature, Debye_T, mineral.params['n'] )
return S
- def enthalpy( self, pressure, temperature, volume, params):
+ def enthalpy(self, mineral):
"""
Returns the enthalpy at the pressure and temperature of the mineral [J/mol]
"""
+ F = mineral.molar_helmholtz(mineral)
+ entropy = mineral.molar_entropy()
+ return F + mineral.temperature * entropy
- return self.helmholtz_free_energy( pressure, temperature, volume, params) + \
- temperature * self.entropy( pressure, temperature, volume, params)
-
- def helmholtz_free_energy( self, pressure, temperature, volume, params):
+ def helmholtz_free_energy(self, mineral):
"""
Returns the Helmholtz free energy at the pressure and temperature of the mineral [J/mol]
"""
- x = params['V_0'] / volume
+ x = mineral.params['V_0'] / mineral.molar_volume()
f = 1./2. * (pow(x, 2./3.) - 1.)
- Debye_T = self.__debye_temperature(params['V_0']/volume, params)
+ Debye_T = self._debye_T(mineral)
- F_quasiharmonic = debye.helmholtz_free_energy( temperature, Debye_T, params['n'] ) - \
- debye.helmholtz_free_energy( 300., Debye_T, params['n'] )
+ F_quasiharmonic = debye.helmholtz_free_energy(mineral.temperature, Debye_T, mineral.params['n'] ) - \
+ debye.helmholtz_free_energy( 300., Debye_T, mineral.params['n'] )
- b_iikk= 9.*params['K_0'] # EQ 28
- b_iikkmm= 27.*params['K_0']*(params['Kprime_0']-4.) # EQ 29
+ b_iikk= 9.*mineral.params['K_0'] # EQ 28
+ b_iikkmm= 27.*mineral.params['K_0']*(mineral.params['Kprime_0']-4.) # EQ 29
- F = params['F_0'] + \
- 0.5*b_iikk*f*f*params['V_0'] + (1./6.)*params['V_0']*b_iikkmm*f*f*f +\
+ F = mineral.params['F_0'] + \
+ 0.5*b_iikk*f*f*mineral.params['V_0'] + (1./6.)*mineral.params['V_0']*b_iikkmm*f*f*f +\
F_quasiharmonic
return F
diff --git a/burnman/mineral.py b/burnman/mineral.py
index f1ece62..44861b1 100644
--- a/burnman/mineral.py
+++ b/burnman/mineral.py
@@ -42,6 +42,7 @@ class Mineral(Material):
self.method = None
if 'equation_of_state' in self.params:
self.set_method(self.params['equation_of_state'])
+ self.eos_data = {}
def set_method(self, equation_of_state):
"""
@@ -95,22 +96,27 @@ class Mineral(Material):
This updates the other properties of this class (v_s, v_p, ...).
"""
-
+ """
#in an effort to avoid additional work, don't do all the calculations if nothing has changed
try:
if self.pressure == pressure and self.temperature == temperature and self.old_params == self.params:
return
except AttributeError:
pass #do nothing
-
+ """
self.pressure = pressure
self.temperature = temperature
- self.old_params = self.params
+ #self.old_params = self.params
+ self.eos_data = {}
if self.method is None:
raise AttributeError, "no method set for mineral, or equation_of_state given in mineral.params"
self.V = self.method.volume(self.pressure, self.temperature, self.params)
+ self.gr = self.K_T = self.K_S = self.G = self.C_v = self.C_p = self.alpha = None
+ self.gibbs = self.helmholtz = self.S = self.H = None
+
+ """
self.gr = self.method.grueneisen_parameter(self.pressure, self.temperature, self.V, self.params)
self.K_T = self.method.isothermal_bulk_modulus(self.pressure, self.temperature, self.V, self.params)
self.K_S = self.method.adiabatic_bulk_modulus(self.pressure, self.temperature, self.V, self.params)
@@ -119,6 +125,7 @@ class Mineral(Material):
self.C_p = self.method.heat_capacity_p(self.pressure, self.temperature, self.V, self.params)
self.alpha = self.method.thermal_expansivity(self.pressure, self.temperature, self.V, self.params)
+
# Attempt to calculate the gibbs free energy and helmholtz free energy, but don't complain if the
# equation of state does not calculate it, or if the mineral params do not have the requisite entries.
try:
@@ -136,7 +143,7 @@ class Mineral(Material):
try:
self.H = self.method.enthalpy(self.pressure, self.temperature, self.V, self.params)
except (KeyError, NotImplementedError):
- self.H = float('nan')
+ self.H = float('nan')"""
# The following gibbs function avoids having to calculate a bunch of unnecessary parameters over P-T space. This will be useful for gibbs minimisation.
@@ -164,42 +171,57 @@ class Mineral(Material):
"""
Returns grueneisen parameter of the mineral [unitless]
"""
+ self.gr = self.gr or self.method.grueneisen_parameter(self)
return self.gr
+
def isothermal_bulk_modulus(self):
"""
Returns isothermal bulk modulus of the mineral [Pa]
"""
+ self.K_T = self.K_T or self.method.isothermal_bulk_modulus(self)
return self.K_T
+
def compressibility(self):
"""
Returns isothermal bulk modulus of the mineral [Pa]
"""
- return 1./self.K_T
+ return 1./self.isothermal_bulk_modulus()
+
def adiabatic_bulk_modulus(self):
"""
Returns adiabatic bulk modulus of the mineral [Pa]
"""
+ self.K_S = self.K_S or self.method.adiabatic_bulk_modulus(self)
return self.K_S
+
def shear_modulus(self):
"""
Returns shear modulus of the mineral [Pa]
"""
+ self.G = self.G or self.method.shear_modulus(self)
return self.G
+
def thermal_expansivity(self):
"""
Returns thermal expansion coefficient of the mineral [1/K]
"""
+ self.alpha = self.alpha or self.method.thermal_expansivity(self)
return self.alpha
+
def heat_capacity_v(self):
"""
Returns heat capacity at constant volume of the mineral [J/K/mol]
"""
+ self.C_v = self.C_v or self.method.heat_capacity_v(self)
return self.C_v
+
def heat_capacity_p(self):
"""
Returns heat capacity at constant pressure of the mineral [J/K/mol]
"""
+ self.C_p = self.C_p or self.method.heat_capacity_p(self)
return self.C_p
+
def v_s(self):
"""
Returns shear wave speed of the mineral [m/s]
@@ -222,22 +244,26 @@ class Mineral(Material):
"""
Returns Gibbs free energy of the mineral [J]
"""
+ self.gibbs = self.gibbs or self.method.gibbs_free_energy(self)
return self.gibbs
def molar_helmholtz(self):
"""
Returns Helmholtz free energy of the mineral [J]
"""
+ self.helmholtz = self.helmholtz or self.method.helmholtz_free_energy(self)
return self.helmholtz
def molar_enthalpy(self):
"""
Returns enthalpy of the mineral [J]
"""
+ self.H = self.H or self.method.enthalpy(self)
return self.H
def molar_entropy(self):
"""
Returns enthalpy of the mineral [J]
"""
+ self.S = self.S or self.method.entropy(self)
return self.S
More information about the CIG-COMMITS
mailing list