Squashed commit of the following:

commit 2454445d3bfb8f08d4e97a02217b2dce82436c9d Merge: 74be206 cb7e671 Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Sat Feb 27 00:04:40 2016 -0500 Merge remote-tracking branch 'upstream/master' into Planets commit 74be206 Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Thu Feb 25 16:13:43 2016 -0500 Cleaned up tests, more to do commit b47c22f Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Wed Feb 24 14:29:51 2016 -0500 Fixed default method in example. Reformatted BM4 to f notation commit cb8e698 Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Wed Feb 24 11:08:25 2016 -0500 Fixed comments commit a83a921 Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Mon Feb 22 19:43:16 2016 -0500 Added tests and benchmarks for Vinet and BM4 commit e151d90 Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Fri Feb 12 12:21:40 2016 -0500 Added vinet to example_user_input commit 8e0a7df Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Thu Jan 7 15:04:30 2016 -0500 Cosmetic fixes commit 167e261 Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Tue Jan 5 14:14:54 2016 -0500 Added K back in. commit 5c515f7 Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Tue Jan 5 14:11:09 2016 -0500 Redundant Code removal commit dded24a Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Tue Jan 5 14:06:41 2016 -0500 Removed K, bracketed V commit be7f28b Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Tue Jan 5 13:59:56 2016 -0500 Cleanup commit cfcd96c Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Tue Jan 5 13:46:29 2016 -0500 Split BM into 4th order commit 2c96d8a Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Tue Jan 5 05:22:09 2016 -0500 Fixed decmials commit 3359fe0 Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Tue Jan 5 05:18:47 2016 -0500 Removed unnecessary Fe mineral commit ec8c5d5 Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Tue Jan 5 05:16:57 2016 -0500 Removing duplicate mineral commit ca412fa Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Tue Jan 5 04:48:48 2016 -0500 Cleaned up failed test commit bea8270 Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Tue Jan 5 04:44:34 2016 -0500 Added a liquid Fe EOS for BM4 Calculations commit 7ba35ed Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Tue Jan 5 04:41:41 2016 -0500 Fixed decimals commit d579344 Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Tue Jan 5 04:35:31 2016 -0500 Added New line Timo Requested commit b0af75f Merge: a8325da d5d3d03 Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Tue Jan 5 04:24:34 2016 -0500 Merge remote-tracking branch 'geodynamics/master' into Planets commit a8325da Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Sun Dec 13 14:09:58 2015 -0800 Fixed uppercase/lowercase issue commit 9de0593 Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Sat Dec 12 23:02:23 2015 -0800 Fixed typo causing test failure commit cf492ab Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Sat Dec 12 22:27:12 2015 -0800 Lowercase V in vines and updated test table commit 200f3b6 Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Sat Dec 12 22:14:58 2015 -0800 Revert "Split Bm4 into it’s own EOS and reverted BM.py" This reverts commit 29a8b20. commit 018d5cb Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Sat Dec 12 19:18:30 2015 -0800 Split Bm4 into it’s own EOS and reverted BM.py commit bb7f6d8 Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Sat Dec 12 18:15:27 2015 -0800 Fixed no new line issue commit fa7efd5 Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Sat Dec 12 18:06:08 2015 -0800 Fixed an import commit 3e830f3 Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Sat Dec 12 18:02:53 2015 -0800 Included the 4th order Birch murnaghan EOS. Added 2 Epsilon Fe minerals in minerals.other and a liquid Fe EOS to use the 4th order BM EOS commit d274159 Author: Cayman Unterborn <kmanunterborn@gmail.com> Date: Thu Dec 10 18:54:27 2015 -0800 Added the Vinet EOS Fixed if to elif
geodynamics · Feb 27, 2016 · 8dac81b · 8dac81b
1 parent cb7e671
commit 8dac81b
Show file tree

Hide file tree

Showing 11 changed files with 450 additions and 10 deletions.
diff --git a/burnman/data/input_figures/Ahmad.png b/burnman/data/input_figures/Ahmad.png
diff --git a/burnman/data/input_figures/Dewaele.png b/burnman/data/input_figures/Dewaele.png
diff --git a/burnman/eos/__init__.py b/burnman/eos/__init__.py
@@ -17,6 +17,8 @@
 from .modified_tait import MT
 from .hp import HP_TMT
 from .cork import CORK
+from .vinet import Vinet
+from .birch_murnaghan_4th import BM4
 
 from .property_modifiers import calculate_property_modifications
 

diff --git a/burnman/eos/birch_murnaghan.py b/burnman/eos/birch_murnaghan.py
@@ -32,7 +32,7 @@ def birch_murnaghan(x, params):
     """
 
     return 3.*params['K_0']/2. * (pow(x, 7./3.) - pow(x, 5./3.)) \
-    * (1 - .75*(4-params['Kprime_0'] )*(pow(x, 2./3.) - 1)) + params['P_0']
+    * (1. - .75*(4.-params['Kprime_0'] )*(pow(x, 2./3.) - 1.)) + params['P_0']
 
 
 def volume(pressure, params):
@@ -68,14 +68,14 @@ def shear_modulus_third_order(volume, params):
 
     x = params['V_0']/volume
     f = 0.5*(pow(x, 2./3.) - 1.0)
-    G = pow((1. + 2*f), 5./2.)*(params['G_0']+(3.*params['K_0']*params['Gprime_0'] - 5.*params['G_0'])*f + (6.*params['K_0']*params['Gprime_0']-24.*params['K_0']-14.*params['G_0']+9./2. * params['K_0']*params['Kprime_0'])*f*f)
+    G = pow((1. + 2.*f), 5./2.)*(params['G_0']+(3.*params['K_0']*params['Gprime_0'] - 5.*params['G_0'])*f + (6.*params['K_0']*params['Gprime_0']-24.*params['K_0']-14.*params['G_0']+9./2. * params['K_0']*params['Kprime_0'])*f*f)
     return G
 
 
 class BirchMurnaghanBase(eos.EquationOfState):
     """
     Base class for the isothermal Birch Murnaghan equation of state.  This is third order in strain, and
-    has no temperature dependence.  However, the shear modulus is sometimes fit to a second order 
+    has no temperature dependence.  However, the shear modulus is sometimes fit to a second order
     function, so if this is the case, you should use that.  For more see :class:`burnman.birch_murnaghan.BM2` and :class:`burnman.birch_murnaghan.BM3`.
     """
     def volume(self,pressure, temperature, params):
@@ -90,7 +90,7 @@ def pressure(self, temperature, volume, params):
     def isothermal_bulk_modulus(self,pressure,temperature, volume, params):
         """
         Returns isothermal bulk modulus :math:`K_T` :math:`[Pa]` as a function of pressure :math:`[Pa]`,
-        temperature :math:`[K]` and volume :math:`[m^3]`. 
+        temperature :math:`[K]` and volume :math:`[m^3]`.
         """
         return bulk_modulus(volume, params)
 
@@ -137,7 +137,7 @@ def validate_parameters(self, params):
         """
         Check for existence and validity of the parameters
         """
-     
+
         if 'P_0' not in params:
             params['P_0'] = 0.
 
@@ -148,13 +148,13 @@ def validate_parameters(self, 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 = ['V_0', 'K_0', 'Kprime_0', 'G_0', 'Gprime_0']
         for k in expected_keys:
             if k not in params:
                 raise KeyError('params object missing parameter : ' + k)
-        
+
         # Finally, check that the values are reasonable.
         if params['P_0'] < 0.:
             warnings.warn( 'Unusual value for P_0', stacklevel=2 )
@@ -173,7 +173,7 @@ def validate_parameters(self, params):
 
 class BM3(BirchMurnaghanBase):
     """
-    Third order Birch Murnaghan isothermal equation of state.  
+    Third order Birch Murnaghan isothermal equation of state.
     This uses the third order expansion for shear modulus.
     """
     def __init__(self):
@@ -182,7 +182,7 @@ def __init__(self):
 
 class BM2(BirchMurnaghanBase):
     """
-    Third order Birch Murnaghan isothermal equation of state.  
+    Third order Birch Murnaghan isothermal equation of state.
     This uses the third order expansion for shear modulus.
     """
     def __init__(self):

diff --git a/burnman/eos/birch_murnaghan_4th.py b/burnman/eos/birch_murnaghan_4th.py
@@ -0,0 +1,142 @@
+from __future__ import absolute_import
+# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for the Earth and Planetary Sciences
+# Copyright (C) 2012 - 2015 by the BurnMan team, released under the GNU GPL v2 or later.
+
+
+import scipy.optimize as opt
+from . import equation_of_state as eos
+from ..tools import bracket
+import warnings
+
+def bulk_modulus_fourth(volume, params):
+    """
+    compute the bulk modulus as per the fourth order
+    birch-murnaghan equation of state.  Returns bulk
+    modulus in the same units as the reference bulk
+    modulus.  Pressure must be in :math:`[Pa]`.
+    """
+
+    x = params['V_0']/volume
+    f = 0.5*(pow(x, 2./3.) - 1.0)
+
+    Xi = (3./4.)*(4.-params['Kprime_0'])
+    Zeta = (3./8.)*((params['K_0']*params['Kprime_prime_0'])+params['Kprime_0']*(params['Kprime_0']-7.)+143./9.)
+
+    K = (5.*f*pow((1.+2.*f),5./2.)*params['K_0']*(1.-(2.*Xi*f)+(4.*Zeta*pow(f,2.)))) + \
+        (pow(1.+(2.*f),7./2.)*params['K_0']*(1.-(4.*Xi*f)+(12.*Zeta*pow(f,2.))))
+
+
+    return K
+
+def volume_fourth_order(pressure,params):
+    func = lambda x: birch_murnaghan_fourth(params['V_0']/x, params) - pressure
+    try:
+        sol = bracket(func, params['V_0'], 1.e-2*params['V_0'])
+    except:
+        raise ValueError('Cannot find a volume, perhaps you are outside of the range of validity for the equation of state?')
+    return opt.brentq(func, sol[0], sol[1])
+
+def birch_murnaghan_fourth(x, params):
+    """
+    equation for the fourth order birch-murnaghan equation of state, returns
+    pressure in the same units that are supplied for the reference bulk
+    modulus (params['K_0'])
+    """
+
+    f = 0.5*(pow(x, 2./3.) - 1.0)
+    Xi = (3./4.)*(4.-params['Kprime_0'])
+    Zeta = (3./8.)*((params['K_0']*params['Kprime_prime_0'])+params['Kprime_0']*(params['Kprime_0']-7.)+143./9.)
+
+    return 3.*f*pow(1.+2.*f,5./2.)*params['K_0']*(1.-(2.*Xi*f)+(4.*Zeta*pow(f,2.)))
+
+
+class BM4(eos.EquationOfState):
+    """
+    Base class for the isothermal Birch Murnaghan equation of state.  This is fourth order in strain, and
+    has no temperature dependence.
+    """
+    def volume(self,pressure, temperature, params):
+        """
+        Returns volume :math:`[m^3]` as a function of pressure :math:`[Pa]`.
+        """
+        return volume_fourth_order(pressure,params)
+
+
+    def pressure(self, temperature, volume, params):
+        return birch_murnaghan_fourth(volume/params['V_0'], params)
+
+
+    def isothermal_bulk_modulus(self,pressure,temperature, volume, params):
+        """
+        Returns isothermal bulk modulus :math:`K_T` :math:`[Pa]` as a function of pressure :math:`[Pa]`,
+        temperature :math:`[K]` and volume :math:`[m^3]`.
+        """
+        return bulk_modulus_fourth(volume,params)
+    def adiabatic_bulk_modulus(self,pressure, temperature, volume, params):
+        """
+        Returns adiabatic bulk modulus :math:`K_s` of the mineral. :math:`[Pa]`.
+        """
+        return bulk_modulus_fourth(volume,params)
+
+    def shear_modulus(self,pressure, temperature, volume, params):
+        """
+        Returns shear modulus :math:`G` of the mineral. :math:`[Pa]`
+        """
+        return 0.
+    def heat_capacity_v(self,pressure, temperature, volume, params):
+        """
+        Since this equation of state does not contain temperature effects, simply return a very large number. :math:`[J/K/mol]`
+        """
+        return 1.e99
+
+    def heat_capacity_p(self,pressure, temperature, volume, params):
+        """
+        Since this equation of state does not contain temperature effects, simply return a very large number. :math:`[J/K/mol]`
+        """
+        return 1.e99
+
+    def thermal_expansivity(self,pressure, temperature, volume, params):
+        """
+        Since this equation of state does not contain temperature effects, simply return zero. :math:`[1/K]`
+        """
+        return 0.
+
+    def grueneisen_parameter(self,pressure,temperature,volume,params):
+        """
+        Since this equation of state does not contain temperature effects, simply return zero. :math:`[unitless]`
+        """
+        return 0.
+
+    def validate_parameters(self, params):
+        """
+        Check for existence and validity of the parameters
+        """
+
+        if 'P_0' not in params:
+            params['P_0'] = 0.
+
+        # 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 '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 = ['V_0', 'K_0', 'Kprime_0']
+        for k in expected_keys:
+            if k not in params:
+                raise KeyError('params object missing parameter : ' + k)
+
+        # Finally, check that the values are reasonable.
+        if params['P_0'] < 0.:
+            warnings.warn( 'Unusual value for P_0', stacklevel=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 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 Kprime_0', stacklevel=2 )
+        if params['Kprime_prime_0'] > 0. or params['Kprime_prime_0'] < -10.:
+            warnings.warn( 'Unusual value for Kprime_prime_0', stacklevel=2 )
diff --git a/burnman/eos/helper.py b/burnman/eos/helper.py
@@ -7,9 +7,11 @@
 from . import slb
 from . import mie_grueneisen_debye as mgd
 from . import birch_murnaghan as bm
+from . import birch_murnaghan_4th as bm4
 from . import modified_tait as mt
 from . import hp 
 from . import cork
+from . import vinet
 from .equation_of_state import EquationOfState
 
 
@@ -21,6 +23,8 @@ def create(method):
     if isinstance(method, str):
         if method == "slb2":
             return slb.SLB2()
+        elif method == "vinet":
+            return vinet.Vinet()
         elif method == "mgd2":
             return mgd.MGD2()
         elif method == "mgd3":
@@ -31,6 +35,8 @@ def create(method):
             return bm.BM2()
         elif method == "bm3":
             return bm.BM3()
+        elif method == "bm4":
+            return bm4.BM4()
         elif method == "mt":
             return mt.MT()
         elif method == "hp_tmt":

diff --git a/burnman/eos/vinet.py b/burnman/eos/vinet.py
@@ -0,0 +1,127 @@
+# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for the Earth and Planetary Sciences
+# Copyright (C) 2012 - 2015 by the BurnMan team, released under the GNU GPL v2 or later.
+
+
+import scipy.optimize as opt
+from . import equation_of_state as eos
+import warnings
+from math import exp
+
+def bulk_modulus(volume, params):
+    """
+    compute the bulk modulus as per the third order
+    Vinet equation of state.  Returns bulk
+    modulus in the same units as the reference bulk
+    modulus.  Pressure must be in :math:`[Pa]`.
+    """
+
+    x = volume/params['V_0']
+    eta = (3./2.)*(params['Kprime_0']-1.)
+
+    K = (params['K_0']*pow(x,-2./3.))*(1+((eta*pow(x,1./3.)+1.)*(1.-pow(x,1./3.))))*exp(eta*(1.-pow(x,1./3.)))
+    return K
+
+def vinet(x, params):
+    """
+    equation for the third order Vinet equation of state, returns
+    pressure in the same units that are supplied for the reference bulk
+    modulus (params['K_0'])
+    """
+    eta = (3./2.)*(params['Kprime_0']-1.)
+    return 3.*params['K_0'] * (pow(x, -2./3.))*(1.-(pow(x,1./3.))) \
+    * exp(eta*(1.-pow(x,1./3.)))
+
+def volume(pressure, params):
+    """
+    Get the Vinet volume at a reference temperature for a given
+    pressure :math:`[Pa]`. Returns molar volume in :math:`[m^3]`
+    """
+
+    func = lambda x: vinet(x/params['V_0'], params) - pressure
+    V = opt.brentq(func, 0.1*params['V_0'], 1.5*params['V_0'])
+    return V
+
+class Vinet(eos.EquationOfState):
+    """
+    Base class for the isothermal Vinet equation of state.  This is third order in strain, and
+    has no temperature dependence.
+    """
+    def volume(self,pressure, temperature, params):
+        """
+        Returns volume :math:`[m^3]` as a function of pressure :math:`[Pa]`.
+        """
+        return volume(pressure,params)
+
+    def pressure(self, temperature, volume, params):
+        return vinet(volume/params['V_0'], params)
+
+    def isothermal_bulk_modulus(self,pressure,temperature, volume, params):
+        """
+        Returns isothermal bulk modulus :math:`K_T` :math:`[Pa]` as a function of pressure :math:`[Pa]`,
+        temperature :math:`[K]` and volume :math:`[m^3]`. 
+        """
+        return bulk_modulus(volume, params)
+
+    def adiabatic_bulk_modulus(self,pressure, temperature, volume, params):
+        """
+        Returns adiabatic bulk modulus :math:`K_s` of the mineral. :math:`[Pa]`.
+        """
+        return bulk_modulus(volume,params)
+
+    def shear_modulus(self,pressure, temperature, volume, params):
+        """
+        Returns shear modulus :math:`G` of the mineral. :math:`[Pa]`
+        Currently not included in the Vinet EOS, so omitted.
+        """
+        return 0.
+
+    def heat_capacity_v(self,pressure, temperature, volume, params):
+        """
+        Since this equation of state does not contain temperature effects, simply return a very large number. :math:`[J/K/mol]`
+        """
+        return 1.e99
+
+    def heat_capacity_p(self,pressure, temperature, volume, params):
+        """
+        Since this equation of state does not contain temperature effects, simply return a very large number. :math:`[J/K/mol]`
+        """
+        return 1.e99
+
+    def thermal_expansivity(self,pressure, temperature, volume, params):
+        """
+        Since this equation of state does not contain temperature effects, simply return zero. :math:`[1/K]`
+        """
+        return 0.
+
+    def grueneisen_parameter(self,pressure,temperature,volume,params):
+        """
+        Since this equation of state does not contain temperature effects, simply return zero. :math:`[unitless]`
+        """
+        return 0.
+
+    def validate_parameters(self, params):
+        """
+        Check for existence and validity of the parameters
+        """
+
+        #G is not included in the Vinet EOS so we shall set them to NaN's
+        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 = ['V_0', 'K_0', 'Kprime_0']
+        for k in expected_keys:
+            if k not in params:
+                raise KeyError('params object missing parameter : ' + k)
+
+        #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['V_0'] < 1.e-7 or params['V_0'] > 1.e-3:
+            warnings.warn( 'Unusual value for V_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'] < -5. or params['Kprime_0'] > 10.:
+            warnings.warn( 'Unusual value for Kprime_0', stacklevel=2 )