Added the Vinet EOS

Fixed if to elif

burnman/eos/helper.py

@@ -10,6 +10,7 @@
 from . import modified_tait as mt
 from . import hp 
 from . import cork
+from . import vinet
 from .equation_of_state import EquationOfState
 
 
@@ -21,6 +22,8 @@ def create(method):
     if isinstance(method, str):
         if method == "slb2":
             return slb.SLB2()
+        elif method == "V":
+            return v.Vinet()
         elif method == "mgd2":
             return mgd.MGD2()
         elif method == "mgd3":
@@ -31,6 +34,8 @@ def create(method):
             return bm.BM2()
         elif method == "bm3":
             return bm.BM3()
+        elif method == "bm4":
+            return bm.BM4()
         elif method == "mt":
             return mt.MT()
         elif method == "hp_tmt":

burnman/eos/vinet.py

@@ -0,0 +1,164 @@
+# 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
+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 density(pressure, params):
+    """
+    Get the Vinet density at a reference temperature for a given
+    pressure :math:`[Pa]`. Returns density in :math:`[kg/m^3]'
+    """
+    return params['molar_mass']/volume(pressure,params)
+
+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
+
+def shear_modulus_second_order(volume, params):
+    """
+    Get the Vinet shear modulus at a reference temperature, for a
+    given volume.  Returns shear modulus in :math:`[Pa]` (the same units as in
+    params['G_0']).  This uses a second order finite strain expansion
+    """
+
+    x = params['V_0']/volume
+    G=params['G_0'] * pow(x,5./3.)*(1.-0.5*(pow(x,2./3.)-1.)*(5.-3.*params['Gprime_0']*params['K_0']/params['G_0']))
+    return G
+
+def shear_modulus_third_order(volume, params):
+    """
+    Get the Vinet shear modulus at a reference temperature, for a
+    given volume.  Returns shear modulus in :math:`[Pa]` (the same units as in
+    params['G_0']).  This uses a third order finite strain expansion
+    """
+
+    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)
+    return G
+
+
+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
+        """
+
+        #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', 'G_0', 'Gprime_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 )
+        if params['G_0'] < 0.0 or params['G_0'] > 1.e13:
+            warnings.warn( 'Unusual value for G_0', stacklevel=2 )
+        if params['Gprime_0'] < -5. or params['Gprime_0'] > 10.:
+            warnings.warn( 'Unusual value for Gprime_0', stacklevel=2 )