[cig-commits] [commit] add_thermodynamic_potentials: Implemented sample_pressure attribute for experimental data inversion (777fd1d)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/burnman

On branch  : add_thermodynamic_potentials
Link       : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51

>---------------------------------------------------------------

commit 777fd1d6a9540e6992346617bade2831b06f4174
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


>---------------------------------------------------------------

777fd1d6a9540e6992346617bade2831b06f4174
 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 8d2a8ed..08c6458 100644
--- a/burnman/mineral.py
+++ b/burnman/mineral.py
@@ -97,6 +97,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].
@@ -122,7 +125,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 184ab7d..5b0ad81 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):
         """