[cig-commits] [commit] add_thermodynamic_potentials: More work on gibbs minimization. Not remotely working (269fc2c)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Dec 9 09:56:32 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_thermodynamic_potentials
Link : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51
>---------------------------------------------------------------
commit 269fc2c5c78eaa552a787fcc3c17c80ee5af595b
Author: ian-r-rose <ian.r.rose at gmail.com>
Date: Wed Sep 17 15:22:36 2014 -0700
More work on gibbs minimization. Not remotely working
>---------------------------------------------------------------
269fc2c5c78eaa552a787fcc3c17c80ee5af595b
burnman/__init__.py | 1 +
burnman/equilibriumassemblage.py | 101 +++++++++++++++++++++++++++++++++++++++
burnman/gibbsminimization.py | 3 +-
test_gibbs.py | 10 ++--
4 files changed, 109 insertions(+), 6 deletions(-)
diff --git a/burnman/__init__.py b/burnman/__init__.py
index ef94f19..7177c9f 100644
--- a/burnman/__init__.py
+++ b/burnman/__init__.py
@@ -34,6 +34,7 @@ from mineral import Mineral
from material import Material
from composite import Composite
from solidsolution import SolidSolution
+from equilibriumassemblage import EquilibriumAssemblage
from mineral_helpers import *
#high level functions
diff --git a/burnman/equilibriumassemblage.py b/burnman/equilibriumassemblage.py
new file mode 100644
index 0000000..e71c3ea
--- /dev/null
+++ b/burnman/equilibriumassemblage.py
@@ -0,0 +1,101 @@
+# BurnMan - a lower mantle toolkit
+# Copyright (C) 2012-2014, Myhill, R., Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
+# Released under GPL v2 or later.
+
+import numpy as np
+import scipy.optimize as opt
+import warnings
+
+import burnman
+import burnman.gibbsminimization as gm
+
+
+class EquilibriumAssemblage(burnman.Material):
+ """
+ Class for taking a certain assemblage of elements,
+ then calculating the equilibrium assemblage of minerals
+ that minimizes the Gibbs free energy. This should
+ have similar functionality to ``burnman.Composite'',
+ but instead of having a fixed set of minerals, it
+ dynamically updates the minerals when you call set_state().
+ As such, it should be significantly slower to use.
+ """
+
+ def __init__(self, composition, phases):
+ """
+ Initialize the equilibrium assemblage with the elements which
+ comprise it as well as a list of phases into which those
+ elements may go.
+
+ Parameters
+ ----------
+ elements : dictionary
+ Dictionary where the keys are strings corresponding to
+ elements (e.g. 'Mg') and the values are molar fractions
+ of that element.
+ phases : list of :class:`burnman.Mineral` or :class:`burnman.SolidSolution`
+ List of phases over which the equilibrium assemblage will try
+ to minimize the gibbs free energy. This class only understands
+ how to use instances of :class:`burnman.Mineral` or
+ :class:`burnman.SolidSolution`.
+ """
+ self.composition = composition
+ self.phases = phases
+
+ stoich, stoich_el, formulae = gm.assemble_stoichiometric_matrix(phases)
+ #The elements in the composition object should be a subset of the
+ #set of elements in the various phases
+ assert( set(composition.keys()).issubset(set(stoich_el) ))
+
+ self.endmember_formulae = formulae
+ self.stoichiometric_matrix = stoich
+ self.elements = stoich_el
+
+ self.bulk_composition_vector = np.array([ (composition[e] if e in composition.keys() else 0.0) \
+ for e in self.elements] )
+
+ def set_method(self, method):
+ for phase in self.phases:
+ phase.set_method(method)
+
+
+ def set_state( self, pressure, temperature):
+
+ n = len(self.endmember_formulae)
+ minimize_gibbs = lambda x : self.__compute_gibbs( pressure, temperature, x )
+ sol = opt.fmin_slsqp( minimize_gibbs, np.ones(n)/n, f_eqcons = self.__composition_constraint, bounds=[(0.0, 1.0),]*n, full_output=0)
+ self.species_vector = sol
+ print zip(self.endmember_formulae, self.species_vector)
+
+ def __compute_gibbs( self, pressure, temperature, species_vector ):
+
+ assert( len(species_vector) == len(self.endmember_formulae) )
+
+ tmp_gibbs = 0.0
+ i = 0
+
+ for phase in self.phases:
+ if isinstance (phase, burnman.SolidSolution):
+ n = len(phase.base_material)
+ total_frac = np.sum( species_vector[i:(i+n)] )
+ phase.set_method('slb3')
+ phase.set_composition( np.array( species_vector[i:(i+n)]/total_frac) )
+ phase.set_state( pressure, temperature )
+ tmp_gibbs += phase.gibbs * total_frac
+ i+=n
+ elif isinstance(phase, burnman.Mineral):
+ phase.set_method('slb3')
+ phase.set_state( pressure, temperature )
+ tmp_gibbs += phase.gibbs * species_vector[i]
+ i+=1
+ else:
+ raise Exception('Unsupported mineral type, can only read burnman.Mineral or burnman.SolidSolution')
+
+ return tmp_gibbs
+
+
+ def __composition_constraint (self, species_vector ):
+ assert( len(species_vector) == len(self.endmember_formulae) )
+ return np.dot( self.stoichiometric_matrix, species_vector) - self.bulk_composition_vector
+
+
diff --git a/burnman/gibbsminimization.py b/burnman/gibbsminimization.py
index 4c5db08..b63a463 100644
--- a/burnman/gibbsminimization.py
+++ b/burnman/gibbsminimization.py
@@ -8,7 +8,7 @@ import burnman
from burnman.processchemistry import *
-def assemble_stoichiometric_matrix ( minerals ):
+def assemble_stoichiometric_matrix ( minerals):
"""
This takes a list of minerals and assembles a matrix where
the rows are elements and the columns are species (or endmembers).
@@ -52,6 +52,7 @@ def assemble_stoichiometric_matrix ( minerals ):
return stoichiometric_matrix, elements, formulae
+
def compute_nullspace ( stoichiometric_matrix ):
"""
Given a stoichiometric matrix, compute a basis for the nullspace.
diff --git a/test_gibbs.py b/test_gibbs.py
index 327f5c2..e231407 100644
--- a/test_gibbs.py
+++ b/test_gibbs.py
@@ -4,9 +4,9 @@ import numpy as np
import matplotlib.pyplot as plt
-minlist = [mg_fe_olivine(), wuestite(), periclase(), stishovite()]
+minlist = [mg_fe_olivine(), wuestite(), periclase(), stishovite(), jadeite()]
+composition = {'Mg': 0.2, 'Si': 0.2, 'Fe': 0.1, 'O': 0.5}
-M, elements, formulae = burnman.assemble_stoichiometric_matrix( minlist )
-print M
-N = burnman.compute_nullspace(M)
-print N
+ea = burnman.EquilibriumAssemblage(composition, minlist)
+ea.set_method('slb3')
+ea.set_state(100.e9, 3000.)
More information about the CIG-COMMITS
mailing list