Merge remote-tracking branch 'geodynamics/master' into Planets

burnman/eos/birch_murnaghan.py

@@ -5,6 +5,7 @@
 
 import scipy.optimize as opt
 from . import equation_of_state as eos
+from ..tools import bracket
 import warnings
 
 def bulk_modulus(volume, params):
@@ -56,8 +57,11 @@ def volume(pressure, params):
     """
 
     func = lambda x: birch_murnaghan(params['V_0']/x, params) - pressure
-    V = opt.brentq(func, 0.5*params['V_0'], 1.5*params['V_0'])
-    return V
+    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 volume_fourth_order(pressure,params):
 

burnman/eos/debye.py

@@ -111,6 +111,7 @@ def debye_fn_cheb(x):
         return ((val_infinity/x)/x)/x;
 
 
+@jit
 def thermal_energy(T, debye_T, n):
     """
     calculate the thermal energy of a substance.  Takes the temperature,
@@ -122,6 +123,7 @@ def thermal_energy(T, debye_T, n):
     E_th = 3.*n*constants.gas_constant*T * debye_fn_cheb(debye_T/T)
     return E_th
 
+@jit
 def heat_capacity_v(T,debye_T,n):
     """
     Heat capacity at constant volume.  In J/K/mol
@@ -132,6 +134,7 @@ def heat_capacity_v(T,debye_T,n):
     C_v = 3.0*n*constants.gas_constant* ( 4.0*debye_fn_cheb(x) - 3.0*x/(np.exp(x)-1.0) )
     return C_v
 
+@jit
 def helmholtz_free_energy(T, debye_T, n):
     """
     Helmholtz free energy of lattice vibrations in the Debye model.

burnman/eos/equation_of_state.py

@@ -324,8 +324,8 @@ def validate_parameters(self, params):
         and an equation of state is not required to implement it.  This function will
         not return anything, though it may raise warnings or errors.
 
-        Params
-        ------
+        Parameters
+        ----------
         params : dictionary
             Dictionary containing material parameters required by the equation of state.
         """

burnman/eos/mie_grueneisen_debye.py

@@ -11,6 +11,7 @@
 from . import birch_murnaghan as bm
 from . import debye
 from .. import constants
+from ..tools import bracket
 
 class MGDBase(eos.EquationOfState):
     """
@@ -38,8 +39,11 @@ def volume(self, pressure,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
-        V = opt.brentq(func, 0.5*params['V_0'], 1.5*params['V_0'])
-        return V
+        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 isothermal_bulk_modulus(self, pressure,temperature,volume, params):
         """
@@ -123,11 +127,14 @@ def pressure(self, temperature, 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
-        return G_th
+        if T > 1.e-10:
+            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
+            return G_th
+        else:
+            return 0.
 
     #compute the Debye temperature in K.  Takes the
     #parameter x, which is V_0/V (molar volumes).
@@ -155,11 +162,14 @@ def _thermal_pressure(self,T,V, params):
     #calculate the thermal correction for the mgd
     #bulk modulus (see matas et al, 2007)
     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 * \
-            ((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
+        if T > 1.e-10:
+            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 * \
+                   ((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
+        else:
+            return 0.
 
 
     def validate_parameters(self, params):

burnman/eos/slb.py

@@ -7,11 +7,45 @@
 import scipy.optimize as opt
 import warnings
 
+# Try to import the jit from numba.  If it is
+# not available, just go with the standard
+# python interpreter
+try:
+    from numba import jit
+except ImportError:
+    def jit(fn):
+        return fn
+
+
 from . import birch_murnaghan as bm
 from . import debye
 from . import equation_of_state as eos
+from ..tools import bracket
+
+
+@jit
+def _grueneisen_parameter_fast(V_0, volume, gruen_0, q_0):
+    """global function with plain parameters so jit will work"""
+    x = V_0 / volume
+    f = 1./2. * (pow(x, 2./3.) - 1.)
+    a1_ii = 6. * gruen_0 # EQ 47
+    a2_iikk = -12.*gruen_0 + 36.*gruen_0*gruen_0 - 18.*q_0*gruen_0 # EQ 47
+    nu_o_nu0_sq = 1.+ a1_ii*f + (1./2.)*a2_iikk * f*f # EQ 41
+    return 1./6./nu_o_nu0_sq * (2.*f+1.) * ( a1_ii + a2_iikk*f )
+
+
+@jit
+def _delta_pressure(x,pressure,temperature,V_0,T_0,Debye_0,n,a1_ii,a2_iikk,b_iikk,b_iikkmm):
 
+    f= 0.5*(pow(V_0/x,2./3.)-1.)
+    debye_temperature =  Debye_0 * np.sqrt(1. + a1_ii * f + 1./2. * a2_iikk*f*f)
+    E_th =  debye.thermal_energy(temperature, debye_temperature, n) #thermal energy at temperature T
+    E_th_ref = debye.thermal_energy(T_0, debye_temperature, n) #thermal energy at reference temperature
+    nu_o_nu0_sq = 1.+ a1_ii*f + (1./2.)*a2_iikk * f*f # EQ 41
+    gr = 1./6./nu_o_nu0_sq * (2.*f+1.) * ( a1_ii + a2_iikk*f )
 
+    return (1./3.)*(pow(1.+2.*f,5./2.))*((b_iikk*f)+(0.5*b_iikkmm*f*f)) \
+            + gr*(E_th - E_th_ref)/x - pressure #EQ 21
 
 class SLBBase(eos.EquationOfState):
     """
@@ -70,29 +104,7 @@ def _thermal_pressure(self,T,V, params):
         P_th = gr * debye.thermal_energy(T,Debye_T, params['n'])/V
         return P_th
 
-    def __volume(self,x,pressure,temperature,V_0,T_0,Debye_0,n,a1_ii,a2_iikk,b_iikk,b_iikkmm):
-
-        f= 0.5*(pow(V_0/x,2./3.)-1.)
-        debye_temperature =  Debye_0 * np.sqrt(1. + a1_ii * f + 1./2. * a2_iikk*f*f)
-        E_th =  debye.thermal_energy(temperature, debye_temperature, n) #thermal energy at temperature T
-        E_th_ref = debye.thermal_energy(T_0, debye_temperature, n) #thermal energy at reference temperature
-        nu_o_nu0_sq = 1.+ a1_ii*f + (1./2.)*a2_iikk * f*f # EQ 41
-        gr = 1./6./nu_o_nu0_sq * (2.*f+1.) * ( a1_ii + a2_iikk*f )
-
-        return (1./3.)*(pow(1.+2.*f,5./2.))*((b_iikk*f)+(0.5*b_iikkmm*f*f)) \
-                + gr*(E_th - E_th_ref)/x - pressure #EQ 21
 
-    def _volume(self,x,pressure,temperature,V_0,T_0,Debye_0,n,a1_ii,a2_iikk,b_iikk,b_iikkmm):
-
-        f= 0.5*(pow(V_0/x,2./3.)-1.)
-        debye_temperature =  Debye_0 * np.sqrt(1. + a1_ii * f + 1./2. * a2_iikk*f*f)
-        E_th =  debye.thermal_energy(temperature, debye_temperature, n) #thermal energy at temperature T
-        E_th_ref = debye.thermal_energy(T_0, debye_temperature, n) #thermal energy at reference temperature
-        nu_o_nu0_sq = 1.+ a1_ii*f + (1./2.)*a2_iikk * f*f # EQ 41
-        gr = 1./6./nu_o_nu0_sq * (2.*f+1.) * ( a1_ii + a2_iikk*f )
-
-        return (1./3.)*(pow(1.+2.*f,5./2.))*((b_iikk*f)+(0.5*b_iikkmm*f*f)) \
-                + gr*(E_th - E_th_ref)/x - pressure #EQ 21
 
     def volume(self, pressure, temperature, params):
         """
@@ -111,22 +123,12 @@ def volume(self, pressure, temperature, params):
 
         # we need to have a sign change in [a,b] to find a zero. Let us start with a
         # conservative guess:
-        a = 0.6*params['V_0']
-        b = 1.2*params['V_0']
-        args=pressure,temperature,V_0,T_0,Debye_0,n,a1_ii,a2_iikk,b_iikk,b_iikkmm
-        # if we have a sign change, we are done:
-        if self._volume(a,pressure,temperature,V_0,T_0,Debye_0,n,a1_ii,a2_iikk,b_iikk,b_iikkmm)*self._volume(b,pressure,temperature,V_0,T_0,Debye_0,n,a1_ii,a2_iikk,b_iikk,b_iikkmm)<0:
-            return opt.brentq(self._volume, a, b,args=args)
-        else:
-            tol = 0.0001
-            sol = opt.fmin(lambda x : self._volume(x,args)*self._volume(x,args), 1.0*params['V_0'], ftol=tol, full_output=1, disp=0)
-            if sol[1] > tol*2:
-                raise ValueError('Cannot find volume, likely outside of the range of validity for EOS')
-            else:
-                warnings.warn("May be outside the range of validity for EOS")
-                return sol[0][0]
-
-
+        args = (pressure,temperature,V_0,T_0,Debye_0,n,a1_ii,a2_iikk,b_iikk,b_iikkmm)
+        try:
+            sol = bracket(_delta_pressure, params['V_0'], 1.e-2*params['V_0'], args)
+        except ValueError:
+            raise Exception('Cannot find a volume, perhaps you are outside of the range of validity for the equation of state?')
+        return opt.brentq(_delta_pressure, sol[0], sol[1] ,args=args)
 
     def pressure( self, temperature, volume, params):
         """
@@ -145,18 +147,11 @@ def pressure( self, temperature, volume, params):
 
         return P
 
-
     def grueneisen_parameter(self, pressure, temperature, volume, params):
         """
         Returns grueneisen parameter :math:`[unitless]` 
         """
-        x = params['V_0'] / volume
-        f = 1./2. * (pow(x, 2./3.) - 1.)
-        gruen_0 = params['grueneisen_0']
-        a1_ii = 6. * gruen_0 # EQ 47
-        a2_iikk = -12.*gruen_0 + 36.*gruen_0*gruen_0 - 18.*params['q_0']*gruen_0 # EQ 47
-        nu_o_nu0_sq = 1.+ a1_ii*f + (1./2.)*a2_iikk * f*f # EQ 41
-        return 1./6./nu_o_nu0_sq * (2.*f+1.) * ( a1_ii + a2_iikk*f )
+        return _grueneisen_parameter_fast(params['V_0'], volume, params['grueneisen_0'], params['q_0'])
 
     def isothermal_bulk_modulus(self, pressure,temperature, volume, params):
         """

burnman/mineral_helpers.py

@@ -17,70 +17,6 @@
 from .composite import Composite
 
 
-
-class HelperSolidSolution(Mineral):
-    """
-    This material is deprecated!
-
-    Class for coming up with a new mineral based based on a solid
-    solution between two or more end member minerals.  It is not
-    completely clear how to do this, or how valid this approximation
-    is, but here we just do a weighted arithmetic average of the
-    thermoelastic properties of the end members according to their molar fractions
-    """
-    def __init__(self, endmembers, molar_fractions):
-        """
-        Takes a list of end member minerals, and a matching list of
-        molar fractions of those minerals for mixing them.  Simply
-        comes up with a new mineral by doing a weighted arithmetic
-        average of the end member minerals
-        """
-        Mineral.__init__(self)
-        self.endmembers = endmembers
-        self.molar_fractions = molar_fractions
-        assert(len(endmembers) == len(molar_fractions))
-        assert(sum(molar_fractions) > 0.9999)
-        assert(sum(molar_fractions) < 1.0001)
-
-        self.method = endmembers[0].method
-
-        #does not make sense to do a solid solution with different number of
-        #atoms per formula unit or different equations of state, at least not simply...
-        for m in endmembers:
-            m.set_method(self.method)
-            if('n' in endmembers[0].params):
-                assert(m.params['n'] == endmembers[0].params['n'])
-
-        self.params = {}
-
-
-    def debug_print(self, indent=""):
-        print("%sHelperSolidSolution(%s):" % (indent, self.to_string()))
-        indent += "  "
-        for (fraction, mat) in zip(self.molar_fractions, self.endmembers):
-            print("%s%g of" % (indent, fraction))
-            mat.debug_print(indent + "  ")
-
-    def set_method(self, method):
-        for mat in self.endmembers:
-            mat.set_method(method)
-        self.method = self.endmembers[0].method
-
-    def set_state(self, pressure, temperature):
-        for mat in self.endmembers:
-            mat.set_state(pressure, temperature)
-
-        itrange = range(0, len(self.endmembers))
-        self.params = {}
-        for prop in self.endmembers[0].params:
-           try:
-                self.params[prop] = sum([ self.endmembers[i].params[prop] * self.molar_fractions[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 endmembers.
-                self.params[prop] = self.endmembers[0].params[prop]
-
-        Mineral.set_state(self,pressure,temperature)
-
 class HelperSpinTransition(Composite):
     """
     Helper class that makes a mineral that switches between two materials
@@ -112,58 +48,3 @@ def set_state(self, pressure, temperature):
             Composite.set_fractions(self, [0.0, 1.0])
 
         Composite.set_state(self, pressure, temperature)
-
-
-
-
-
-
-
-
-class HelperFeDependent(Mineral):
-
-    """
-    This material is deprecated!
-
-    Helper to implement a rock that does iron exchange (two minerals with
-    changing iron content based on P,T).
-
-    Classes deriving from this helper need to implement
-    create_inner_material() and provide a function in __init__ that computes
-    the iron exchange.
-    """
-
-    def __init__(self, iron_number_with_pt, idx):
-        Material.__init__(self)
-        self.iron_number_with_pt = iron_number_with_pt
-        self.which_index = idx  # take input 0 or 1 from iron_number_with_pt()
-        self.method = None
-
-    def debug_print(self, indent=""):
-        print("%sHelperFeDependent:" % indent)
-        print("%s  %s" % (indent, self.to_string()))
-
-    def create_inner_material(self, iron_number):
-        raise NotImplementedError("need to implement create_inner_material() in derived class!")
-        return None
-
-    def set_method(self, method):
-        self.method = method
-
-    def iron_number(self):
-        return self.iron_number_with_pt(self.pressure, self.temperature)[self.which_index]
-
-    def set_state(self, pressure, temperature):
-        Material.set_state(self, pressure, temperature)
-        self.endmembers = self.create_inner_material(self.iron_number())
-        if self.method:
-            self.endmembers.set_method(self.method)
-        self.endmembers.set_state(pressure, temperature)
-
-    def unroll(self):
-        """ return (fractions, minerals) where both are arrays. May depend on current state """
-        return ([self.endmembers],[1.0])
-
-
-
-

burnman/minerals/Murakami_2013.py

@@ -54,20 +54,6 @@ def __init__(self):
 
         Mineral.__init__(self)
 
-class ferropericlase(helpers.HelperSolidSolution):
-    def __init__(self, fe_num):
-        endmembers = [periclase(), wuestite()]
-        molar_fractions = [1. - fe_num, 0.0 + fe_num] # keep the 0.0 +, otherwise it is an array sometimes
-        helpers.HelperSolidSolution.__init__(self, endmembers, molar_fractions)
-
-
-class mg_fe_perovskite(helpers.HelperSolidSolution):
-    def __init__(self, fe_num):
-        endmembers = [mg_perovskite(), fe_perovskite()]
-        molar_fractions = [1. - fe_num, 0.0 + fe_num] # keep the 0.0 +, otherwise it is an array sometimes
-        helpers.HelperSolidSolution.__init__(self, endmembers, molar_fractions)
-
-
 class mg_perovskite(Mineral):
     def __init__(self):
         self.params = {
@@ -106,4 +92,3 @@ def __init__(self):
 
 mg_bridgmanite = mg_perovskite
 fe_bridgmanite = fe_perovskite
-mg_fe_bridgmanite = mg_fe_perovskite