[cig-commits] [commit] add_thermodynamic_potentials: Adding non-functional solidsolution.py (5b03886)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Dec 9 09:53:31 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_thermodynamic_potentials
Link : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51
>---------------------------------------------------------------
commit 5b038869cdc2e926c84a047f1f7e27663d7ff099
Author: Bob Myhill <myhill.bob at gmail.com>
Date: Tue Jul 29 00:24:54 2014 -0700
Adding non-functional solidsolution.py
>---------------------------------------------------------------
5b038869cdc2e926c84a047f1f7e27663d7ff099
.../{mie_grueneisen_debye.py => modified_tait.py} | 9 +-
burnman/{mineral.py => solidsolution.py} | 116 ++++++++++++++-------
2 files changed, 82 insertions(+), 43 deletions(-)
diff --git a/burnman/mie_grueneisen_debye.py b/burnman/modified_tait.py
similarity index 95%
copy from burnman/mie_grueneisen_debye.py
copy to burnman/modified_tait.py
index 9b3a1e5..cc6786b 100644
--- a/burnman/mie_grueneisen_debye.py
+++ b/burnman/modified_tait.py
@@ -1,5 +1,5 @@
# BurnMan - a lower mantle toolkit
-# Copyright (C) 2012, 2013, Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
+# Copyright (C) 2012, 2013, Myhill, R., Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
# Released under GPL v2 or later.
import numpy as np
@@ -12,10 +12,9 @@ import burnman.debye as debye
class MGDBase(eos.EquationOfState):
"""
- Base class for a generic finite-strain-mie-grueneisen-debye
- equation of state. References for this can be found in many
- places, such as Shim, Duffym and Kenichi (2002) and Jackson and Rigedn
- (1996). Here we mostly follow the appendices of Matas et al (2007)
+ Base class for a generic modified Tait equation of state.
+ References for this can be found in Huang and Chow (1974)
+ and Holland and Powell (2011; followed here).
"""
#def __init__(self):
diff --git a/burnman/mineral.py b/burnman/solidsolution.py
similarity index 55%
copy from burnman/mineral.py
copy to burnman/solidsolution.py
index e641665..ed8bf33 100644
--- a/burnman/mineral.py
+++ b/burnman/solidsolution.py
@@ -6,54 +6,97 @@ import warnings
import numpy as np
-from burnman.material import Material
-import burnman.equation_of_state as eos
-import burnman.birch_murnaghan as bm
-import burnman.slb as slb
-import burnman.mie_grueneisen_debye as mgd
+# Kronecker delta
+def kd(x, y):
+ if x == y:
+ return 1
+ else:
+ return 0
-class Mineral(Material):
+import burnman.symmetric_model as sym
+import burnman.asymmetric_model as asym
+
+class SolidSolution(Mineral):
"""
- This is the base class for all minerals. States of the mineral
- can only be queried after setting the pressure and temperature
- using set_state(). The method for computing properties of
- the material is set using set_method(), which should be done
+ This is the base class for all solid solutions.
+ States of the solid solution can only be queried after setting
+ the pressure and temperature using set_state().
+ The method for computing properties of
+ the solution is set using set_method(), which should be done
once after creating the material.
- This class is available as ``burnman.Mineral``.
+ This class is available as ``burnman.SolidSolution``.
If deriving from this class, set the properties in self.params
to the desired values. For more complicated materials you
can overwrite set_state(), change the params and then call
set_state() from this class.
- All the material parameters are expected to be in plain SI units. This
- means that the elastic moduli should be in Pascals and NOT Gigapascals,
- and the Debye temperature should be in K not C. Additionally, the
- reference volume should be in m^3/(mol molecule) and not in unit cell
- volume and 'n' should be the number of atoms per molecule. Frequently in
- the literature the reference volume is given in Angstrom^3 per unit cell.
- To convert this to m^3/(mol molecule) you should multiply by 10^(-30) *
- N_a / Z, where N_a is Avogadro's number and Z is the number of atoms per
- unit cell. You can look up Z in many places, including www.mindat.org
+ All the solid solution parameters are expected to be in SI units. This
+ means that the interaction parameters should be in J/mol, with the T
+ and P derivatives in J/K/mol and m^3/(mol molecule).
"""
- def __init__(self):
- self.params = { 'name':'generic',
- 'equation_of_state': 'slb3', #Equation of state used to fit the parameters
- 'V_0': 0., #Molar volume [m^3/(mole molecules)] at room pressure/temperature
- 'K_0': 0., #Reference bulk modulus [Pa] at room pressure/temperature
- 'Kprime_0': 0., #pressure derivative of bulk modulus
- 'G_0': 0., #reference shear modulus at room pressure/temperature
- 'Gprime_0': 0., #pressure derivative of shear modulus
- 'molar_mass': 0., #molar mass in units of [kg/mol]
- 'n': 0., #number of atoms per molecule
- 'Debye_0': 0., #Debye temperature for material. See Stixrude & Lithgow-Bertelloni, 2005 for values
- 'grueneisen_0': 0., #Gruneisen parameter for material. See Stixrude & Lithgow-Bertelloni, 2005 for values
- 'q_0': 0., #q value used in caluclations. See Stixrude & Lithgow-Bertelloni, 2005 for values
- 'eta_s_0': 0.0} #eta value used in calculations. See Stixrude & Lithgow-Bertelloni, 2005 for values
- self.method = None
+ def __init__(self, base_material, molar_fraction, site_occupancy, interaction_parameter, van_laar_parameter):
+ self.base_material = base_material
+ self.molar_fraction = molar_fraction
+ assert(len(base_material) == len(molar_fraction))
+ assert(len(base_material) == len(van_laar_parameter))
+ assert(sum(molar_fraction) > 0.9999)
+ assert(sum(molar_fraction) < 1.0001)
+
+ # make sure the number of atoms in each endmember is the same
+ # (when we start adding minerals with vacancies
+ # this will need to be changed).
+ for m in base_material:
+ if(base_material[0].params.has_key('n')):
+ assert(m.params['n'] == base_material[0].params['n'])
+
+ # Set the state of all the endmembers
+ def set_state(self, pressure, temperature):
+ for mat in self.base_materials:
+ mat.method = self.method
+ mat.set_state(pressure, temperature)
+
+ itrange = range(0, len(self.base_material))
+
+ self.params = {}
+
+ # NEEDS CHANGING FROM HERE!!!
+ # Description of RTlngamma terms given by the asymmetric formalism
+ # is included in Holland and Powell (2003)
+
+ # the volume is calculated in the same way as Gex (p.493),
+ # replacing Wij with Wvij
+ Vex=-2.*sum([ sum([ self.molar_fraction[i]*self.van_laar_parameter[i]*self.molar_fraction[j]*self.van_laar_parameter[j]*self.interaction_parameter[i][2] / (self.van_laar_parameter[i] + self.van_laar_parameter[j]) for j in range(i, len(self.base_material)) ]) for i in range(0, len(self.base_material)-1) ])/sum([ self.van_laar_parameter[l]*self.molar_fraction[l] for l in itrange ])
+ V= sum([ self.base_materials[i].molar_volume() * self.molar_fraction[i] for i in itrange ]) + Vex
+
+ # the volume is calculated in the same way as Gex (p.493),
+ # replacing Wij with Wvij
+ # NEED TO CORRECT INDICES FOR INT PARAM
+ # W[i][j] == W[(itrange-1)*i+(j-1)]
+ Vex=-2.*sum([ sum([ self.molar_fraction[i]*self.van_laar_parameter[i]*self.molar_fraction[j]*self.van_laar_parameter[j]*self.interaction_parameter[(itrange-1)*i+(j-1)][2] / (self.van_laar_parameter[i] + self.van_laar_parameter[j]) for j in range(i, len(self.base_material)) ]) for i in range(0, len(self.base_material)-1) ])/sum([ self.van_laar_parameter[l]*self.molar_fraction[l] for l in itrange ])
+ V= sum([ self.base_materials[i].molar_volume() * self.molar_fraction[i] for i in itrange ]) + Vex
+
+ # Same for the Gibbs free energy
+ # NEED TO CORRECT INDICES FOR INT PARAM
+ Gex=-2.*sum([ sum([ self.molar_fraction[i]*self.van_laar_parameter[i]*self.molar_fraction[j]*self.van_laar_parameter[j]*(pressure*self.interaction_parameter[(itrange-1)*i+(j-1)][0] + temperature*self.interaction_parameter[(itrange-1)*i+(j-1)][1] + pressure*self.interaction_parameter[(itrange-1)*i+(j-1)][2]) / (self.van_laar_parameter[i] + self.van_laar_parameter[j]) for j in range(i, len(self.base_material)) ]) for i in range(0, len(self.base_material)-1) ])/sum([ self.van_laar_parameter[l]*self.molar_fraction[l] for l in itrange ])
+ G= sum([ self.base_materials[i].molar_gibbs() * self.molar_fraction[i] for i in itrange ]) + Gex
+
+
+
+
+
+ for prop in self.base_materials[0].params:
+ try:
+ self.params[prop] = sum([ self.base_materials[i].params[prop] * self.molar_fraction[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 base_material.
+ self.params[prop] = self.base_materials[0].params[prop]
+ Mineral.set_state(self, pressure, temperature)
+
+ # BELOW IS FROM THE MINERAL CLASS.
def set_method(self, method):
"""
Set the equation of state to be used for this mineral.
@@ -88,9 +131,6 @@ class Mineral(Material):
"""
return "'" + self.__class__.__module__.replace(".minlib_",".") + "." + self.__class__.__name__ + "'"
- def debug_print(self, indent=""):
- print "%s%s" % (indent, self.to_string())
-
def unroll(self):
return ([1.0],[self])
More information about the CIG-COMMITS
mailing list