[cig-commits] [commit] inversion, master, validate_MT_params: Adding non-functional solidsolution.py (a979bde)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Fri Dec 12 18:24:17 PST 2014
Repository : https://github.com/geodynamics/burnman
On branches: inversion,master,validate_MT_params
Link : https://github.com/geodynamics/burnman/compare/80c2a295c42dfdb38f83f6c1334bf7d8f97a8463...409647ff05dfad6a686198cac1481bd46b5e2e62
>---------------------------------------------------------------
commit a979bde5d7672d014b78bea93920a3e985c910e7
Author: Bob Myhill <myhill.bob at gmail.com>
Date: Tue Jul 29 00:24:54 2014 -0700
Adding non-functional solidsolution.py
>---------------------------------------------------------------
a979bde5d7672d014b78bea93920a3e985c910e7
.../{mie_grueneisen_debye.py => modified_tait.py} | 29 ++-
burnman/solidsolution.py | 237 +++++++++++++++++++++
2 files changed, 249 insertions(+), 17 deletions(-)
diff --git a/burnman/mie_grueneisen_debye.py b/burnman/modified_tait.py
similarity index 85%
copy from burnman/mie_grueneisen_debye.py
copy to burnman/modified_tait.py
index fd606b2..cc6786b 100644
--- a/burnman/mie_grueneisen_debye.py
+++ b/burnman/modified_tait.py
@@ -1,5 +1,5 @@
# BurnMan - a lower mantle toolkit
-# Copyright (C) 2012, 2013, Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
+# Copyright (C) 2012, 2013, Myhill, R., Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
# Released under GPL v2 or later.
import numpy as np
@@ -8,14 +8,13 @@ import scipy.optimize as opt
import burnman.equation_of_state as eos
import burnman.birch_murnaghan as bm
import burnman.debye as debye
-import constants
+
class MGDBase(eos.EquationOfState):
"""
- Base class for a generic finite-strain-mie-grueneisen-debye
- equation of state. References for this can be found in many
- places, such as Shim, Duffym and Kenichi (2002) and Jackson and Rigedn
- (1996). Here we mostly follow the appendices of Matas et al (2007)
+ 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 __init__(self):
@@ -33,10 +32,9 @@ class MGDBase(eos.EquationOfState):
Returns volume [m^3] as a function of pressure [Pa] and temperature [K]
EQ B7
"""
- T_0 = self.reference_temperature( params )
func = lambda x: bm.birch_murnaghan(params['V_0']/x, params) + \
self.__thermal_pressure(temperature, x, params) - \
- self.__thermal_pressure(T_0, x, params) - pressure
+ self.__thermal_pressure(300., x, params) - pressure
V = opt.brentq(func, 0.5*params['V_0'], 1.5*params['V_0'])
return V
@@ -45,10 +43,9 @@ class MGDBase(eos.EquationOfState):
Returns isothermal bulk modulus [Pa] as a function of pressure [Pa],
temperature [K], and volume [m^3]. EQ B8
"""
- T_0 = self.reference_temperature( params )
K_T = bm.bulk_modulus(volume, params) + \
self.__thermal_bulk_modulus(temperature,volume, params) - \
- self.__thermal_bulk_modulus(T_0,volume, params) #EQB13
+ self.__thermal_bulk_modulus(300.,volume, params) #EQB13
return K_T
#calculate the mgd shear modulus as a function of P, V, and T
@@ -57,15 +54,14 @@ class MGDBase(eos.EquationOfState):
Returns shear modulus [Pa] as a function of pressure [Pa],
temperature [K], and volume [m^3]. EQ B11
"""
- T_0 = self.reference_temperature( params )
if self.order==2:
return bm.shear_modulus_second_order(volume,params) + \
self.__thermal_shear_modulus(temperature,volume, params) - \
- self.__thermal_shear_modulus(T_0,volume, params) # EQ B11
+ self.__thermal_shear_modulus(300.,volume, params) # EQ B11
elif self.order==3:
return bm.shear_modulus_third_order(volume,params) + \
self.__thermal_shear_modulus(temperature,volume, params) - \
- self.__thermal_shear_modulus(T_0,volume, params) # EQ B11
+ self.__thermal_shear_modulus(300.,volume, params) # EQ B11
else:
raise NotImplementedError("")
@@ -115,17 +111,16 @@ class MGDBase(eos.EquationOfState):
Returns pressure [Pa] as a function of temperature [K] and volume[m^3]
EQ B7
"""
- T_0 = self.reference_temperature( params )
return bm.birch_murnaghan(params['V_0']/volume, params) + \
self.__thermal_pressure(temperature,volume, params) - \
- self.__thermal_pressure(T_0,volume, params)
+ self.__thermal_pressure(300.,volume, params)
#calculate the thermal correction to the shear modulus as a function of V, T
def __thermal_shear_modulus(self, T, V, params):
gr = self.__grueneisen_parameter(params['V_0']/V, params)
Debye_T = self.__debye_temperature(params['V_0']/V, params)
G_th= 3./5. * ( self.__thermal_bulk_modulus(T,V,params) - \
- 6*constants.gas_constant*T*params['n']/V * gr * debye.debye_fn(Debye_T/T) ) # EQ B10
+ 6*debye.R*T*params['n']/V * gr * debye.debye_fn(Debye_T/T) ) # EQ B10
return G_th
#compute the Debye temperature in K. Takes the
@@ -156,7 +151,7 @@ class MGDBase(eos.EquationOfState):
def __thermal_bulk_modulus(self, T,V,params):
gr = self.__grueneisen_parameter(params['V_0']/V, params)
Debye_T = self.__debye_temperature(params['V_0']/V, params)
- K_th = 3.*params['n']*constants.gas_constant*T/V * gr * \
+ K_th = 3.*params['n']*debye.R*T/V * gr * \
((1. - params['q_0'] - 3.*gr)*debye.debye_fn(Debye_T/T)+3.*gr*(Debye_T/T)/(np.exp(Debye_T/T) - 1.)) # EQ B5
return K_th
diff --git a/burnman/solidsolution.py b/burnman/solidsolution.py
new file mode 100644
index 0000000..ed8bf33
--- /dev/null
+++ b/burnman/solidsolution.py
@@ -0,0 +1,237 @@
+# BurnMan - a lower mantle toolkit
+# Copyright (C) 2012, 2013, Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
+# Released under GPL v2 or later.
+
+import warnings
+
+import numpy as np
+
+# Kronecker delta
+def kd(x, y):
+ if x == y:
+ return 1
+ else:
+ return 0
+
+import burnman.symmetric_model as sym
+import burnman.asymmetric_model as asym
+
+class SolidSolution(Mineral):
+ """
+ This is the base class for all solid solutions.
+ States of the solid solution can only be queried after setting
+ the pressure and temperature using set_state().
+ The method for computing properties of
+ the solution is set using set_method(), which should be done
+ once after creating the material.
+
+ This class is available as ``burnman.SolidSolution``.
+
+ If deriving from this class, set the properties in self.params
+ to the desired values. For more complicated materials you
+ can overwrite set_state(), change the params and then call
+ set_state() from this class.
+
+ All the solid solution parameters are expected to be in SI units. This
+ means that the interaction parameters should be in J/mol, with the T
+ and P derivatives in J/K/mol and m^3/(mol molecule).
+ """
+
+ def __init__(self, base_material, molar_fraction, site_occupancy, interaction_parameter, van_laar_parameter):
+ self.base_material = base_material
+ self.molar_fraction = molar_fraction
+ assert(len(base_material) == len(molar_fraction))
+ assert(len(base_material) == len(van_laar_parameter))
+ assert(sum(molar_fraction) > 0.9999)
+ assert(sum(molar_fraction) < 1.0001)
+
+ # make sure the number of atoms in each endmember is the same
+ # (when we start adding minerals with vacancies
+ # this will need to be changed).
+ for m in base_material:
+ if(base_material[0].params.has_key('n')):
+ assert(m.params['n'] == base_material[0].params['n'])
+
+ # Set the state of all the endmembers
+ def set_state(self, pressure, temperature):
+ for mat in self.base_materials:
+ mat.method = self.method
+ mat.set_state(pressure, temperature)
+
+ itrange = range(0, len(self.base_material))
+
+ self.params = {}
+
+ # NEEDS CHANGING FROM HERE!!!
+ # Description of RTlngamma terms given by the asymmetric formalism
+ # is included in Holland and Powell (2003)
+
+ # the volume is calculated in the same way as Gex (p.493),
+ # replacing Wij with Wvij
+ Vex=-2.*sum([ sum([ self.molar_fraction[i]*self.van_laar_parameter[i]*self.molar_fraction[j]*self.van_laar_parameter[j]*self.interaction_parameter[i][2] / (self.van_laar_parameter[i] + self.van_laar_parameter[j]) for j in range(i, len(self.base_material)) ]) for i in range(0, len(self.base_material)-1) ])/sum([ self.van_laar_parameter[l]*self.molar_fraction[l] for l in itrange ])
+ V= sum([ self.base_materials[i].molar_volume() * self.molar_fraction[i] for i in itrange ]) + Vex
+
+ # the volume is calculated in the same way as Gex (p.493),
+ # replacing Wij with Wvij
+ # NEED TO CORRECT INDICES FOR INT PARAM
+ # W[i][j] == W[(itrange-1)*i+(j-1)]
+ Vex=-2.*sum([ sum([ self.molar_fraction[i]*self.van_laar_parameter[i]*self.molar_fraction[j]*self.van_laar_parameter[j]*self.interaction_parameter[(itrange-1)*i+(j-1)][2] / (self.van_laar_parameter[i] + self.van_laar_parameter[j]) for j in range(i, len(self.base_material)) ]) for i in range(0, len(self.base_material)-1) ])/sum([ self.van_laar_parameter[l]*self.molar_fraction[l] for l in itrange ])
+ V= sum([ self.base_materials[i].molar_volume() * self.molar_fraction[i] for i in itrange ]) + Vex
+
+
+ # Same for the Gibbs free energy
+ # NEED TO CORRECT INDICES FOR INT PARAM
+ Gex=-2.*sum([ sum([ self.molar_fraction[i]*self.van_laar_parameter[i]*self.molar_fraction[j]*self.van_laar_parameter[j]*(pressure*self.interaction_parameter[(itrange-1)*i+(j-1)][0] + temperature*self.interaction_parameter[(itrange-1)*i+(j-1)][1] + pressure*self.interaction_parameter[(itrange-1)*i+(j-1)][2]) / (self.van_laar_parameter[i] + self.van_laar_parameter[j]) for j in range(i, len(self.base_material)) ]) for i in range(0, len(self.base_material)-1) ])/sum([ self.van_laar_parameter[l]*self.molar_fraction[l] for l in itrange ])
+ G= sum([ self.base_materials[i].molar_gibbs() * self.molar_fraction[i] for i in itrange ]) + Gex
+
+
+
+
+
+ for prop in self.base_materials[0].params:
+ try:
+ self.params[prop] = sum([ self.base_materials[i].params[prop] * self.molar_fraction[i] for i in itrange ])
+ except TypeError:
+ #if there is a type error, it is probably a string. Just go with the value of the first base_material.
+ self.params[prop] = self.base_materials[0].params[prop]
+ Mineral.set_state(self, pressure, temperature)
+
+ # BELOW IS FROM THE MINERAL CLASS.
+ def set_method(self, method):
+ """
+ Set the equation of state to be used for this mineral.
+ Takes a string corresponding to any of the predefined
+ equations of state: 'bm2', 'bm3', 'mgd2', 'mgd3', 'slb2',
+ or 'slb3'. Alternatively, you can pass a user defined
+ class which derives from the equation_of_state base class.
+ """
+ if( isinstance(method, basestring)):
+ if (method == "slb2"):
+ self.method = slb.SLB2()
+ elif (method == "mgd2"):
+ self.method = mgd.MGD2()
+ elif (method == "mgd3"):
+ self.method = mgd.MGD3()
+ elif (method == "slb3"):
+ self.method = slb.SLB3()
+ elif (method == "bm2"):
+ self.method = bm.BM2()
+ elif (method == "bm3"):
+ self.method = bm.BM3()
+ else:
+ raise Exception("unsupported material method " + method)
+ elif ( issubclass(method, eos.EquationOfState) ):
+ self.method = method()
+ else:
+ raise Exception("unsupported material method " + method.__class__.__name__ )
+
+ def to_string(self):
+ """
+ Returns the name of the mineral class
+ """
+ return "'" + self.__class__.__module__.replace(".minlib_",".") + "." + self.__class__.__name__ + "'"
+
+ def unroll(self):
+ return ([1.0],[self])
+
+ def set_state(self, pressure, temperature):
+ """
+ Update the material to the given pressure [Pa] and temperature [K].
+
+ 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.V = self.method.volume(self.pressure, self.temperature, self.params)
+ 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)
+ self.C_v = self.method.heat_capacity_v(self.pressure, self.temperature, self.V, self.params)
+ 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)
+
+ if (self.params.has_key('G_0') and self.params.has_key('Gprime_0')):
+ self.G = self.method.shear_modulus(self.pressure, self.temperature, self.V, self.params)
+ else:
+ self.G = float('nan') #nan if there is no G, this should propagate through calculations to the end
+ warnings.warn(('Warning: G_0 and or Gprime_0 are undefined for ' + self.to_string()))
+
+ def molar_mass(self):
+ """
+ Returns molar mass of the mineral [kg/mol]
+ """
+ return self.params['molar_mass']
+
+ def density(self):
+ """
+ Returns density of the mineral [kg/m^3]
+ """
+ return self.params['molar_mass'] / self.V
+
+ def molar_volume(self):
+ """
+ Returns molar volume of the mineral [m^3/mol]
+ """
+ return self.V
+ def grueneisen_parameter(self):
+ """
+ Returns grueneisen parameter of the mineral [unitless]
+ """
+ return self.gr
+ def isothermal_bulk_modulus(self):
+ """
+ Returns isothermal bulk modulus of the mineral [Pa]
+ """
+ return self.K_T
+ def adiabatic_bulk_modulus(self):
+ """
+ Returns adiabatic bulk modulus of the mineral [Pa]
+ """
+ return self.K_S
+ def shear_modulus(self):
+ """
+ Returns shear modulus of the mineral [Pa]
+ """
+ return self.G
+ def thermal_expansivity(self):
+ """
+ Returns thermal expansion coefficient of the mineral [1/K]
+ """
+ return self.alpha
+ def heat_capacity_v(self):
+ """
+ Returns heat capacity at constant volume of the mineral [J/K/mol]
+ """
+ return self.C_v
+ def heat_capacity_p(self):
+ """
+ Returns heat capacity at constant pressure of the mineral [J/K/mol]
+ """
+ return self.C_p
+ def v_s(self):
+ """
+ Returns shear wave speed of the mineral [m/s]
+ """
+ return np.sqrt(self.shear_modulus() / \
+ self.density())
+ def v_p(self):
+ """
+ Returns P wave speed of the mineral [m/s]
+ """
+ return np.sqrt((self.adiabatic_bulk_modulus() + 4. / 3. * \
+ self.shear_modulus()) / self.density())
+ def v_phi(self):
+ """
+ Returns bulk sound speed of the mineral [m/s]
+ """
+ return np.sqrt(self.adiabatic_bulk_modulus() / self.density())
More information about the CIG-COMMITS
mailing list