[cig-commits] [commit] add_gibbs_energy: Implemented sample_pressure attribute for experimental data inversion (6744559)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Dec 11 17:10:26 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_gibbs_energy
Link : https://github.com/geodynamics/burnman/compare/0000000000000000000000000000000000000000...2148b324d3e8aa7b527f831eb397590942563008
>---------------------------------------------------------------
commit 67445590d5692ee5de14409cae2e0129fd35aacc
Author: Bob Myhill <myhill.bob at gmail.com>
Date: Tue Aug 26 14:38:01 2014 +0200
Implemented sample_pressure attribute for experimental data inversion
>---------------------------------------------------------------
67445590d5692ee5de14409cae2e0129fd35aacc
burnman/P-V_test.py | 26 ++++++++++++++++++++++++++
burnman/mineral.py | 6 +++++-
burnman/slb.py | 13 +++++++++++++
3 files changed, 44 insertions(+), 1 deletion(-)
diff --git a/burnman/P-V_test.py b/burnman/P-V_test.py
new file mode 100644
index 0000000..a31e6e9
--- /dev/null
+++ b/burnman/P-V_test.py
@@ -0,0 +1,26 @@
+import os, sys
+if not os.path.exists('burnman') and os.path.exists('../burnman'):
+ sys.path.insert(1,os.path.abspath('..'))
+
+import burnman
+from burnman import minerals
+
+P0=2.e10
+T=2273.15
+
+wad=minerals.other.Katsura_2009_wadsleyite()
+wad.set_method("slb3")
+
+# Here we find the volume at the P, T of interest
+wad.set_state(P0,T)
+print 'V @', T, 'K,', P0, 'Pa:'
+print wad.V, 'm^3/mol'
+print ''
+
+# Here we take the V from the previous step and calculate the pressure. Hopefully it's the same as our original input!!
+V=wad.V
+
+P1=wad.sample_pressure(T,V)
+print 'Retrieved pressure:', P1, 'Pa'
+print 'Pressure difference:', P1-P0, 'Pa'
+print 'Fractional difference:', (P1-P0)/P0
diff --git a/burnman/mineral.py b/burnman/mineral.py
index ced6490..c560e11 100644
--- a/burnman/mineral.py
+++ b/burnman/mineral.py
@@ -110,6 +110,9 @@ class Mineral(Material):
def unroll(self):
return ([1.0],[self])
+ def sample_pressure(self, temperature, volume):
+ return self.method.sample_pressure(temperature, volume, self.params)
+
def set_state(self, pressure, temperature):
"""
Update the material to the given pressure [Pa] and temperature [K].
@@ -138,7 +141,8 @@ class Mineral(Material):
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)
- self.gibbs = self.method.gibbs(self.pressure, self.temperature, self.params)
+ if (self.params.has_key('H_0') and self.params.has_key('S_0')):
+ self.gibbs = self.method.gibbs(self.pressure, self.temperature, 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)
diff --git a/burnman/slb.py b/burnman/slb.py
index 3c2fab0..2b494cd 100644
--- a/burnman/slb.py
+++ b/burnman/slb.py
@@ -16,6 +16,7 @@ class SLBBase(eos.EquationOfState):
in Stixrude and Lithgow-Bertelloni (2005). 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]
@@ -53,6 +54,18 @@ class SLBBase(eos.EquationOfState):
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
+ def sample_pressure(self, temperature, volume, params):
+ return bm.birch_murnaghan(params['V_0']/volume, params) + \
+ self.__thermal_pressure(temperature,volume, params) - \
+ self.__thermal_pressure(300.,volume, params)
+
+ #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., T, V, params) # P not important
+ P_th = gr * debye.thermal_energy(T,Debye_T, params['n'])/V
+ return P_th
def volume(self, pressure, temperature, params):
"""
More information about the CIG-COMMITS
mailing list