[cig-commits] [commit] add_thermodynamic_potentials: Added chemical potential and fugacity functions and examples (5ef1a29)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/burnman

On branch  : add_thermodynamic_potentials
Link       : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51

>---------------------------------------------------------------

commit 5ef1a29b106ff6d68671f20dfe7d2c93241d3720
Author: Bob Myhill <myhill.bob at gmail.com>
Date:   Wed Sep 3 17:26:12 2014 +0200

    Added chemical potential and fugacity functions and examples


>---------------------------------------------------------------

5ef1a29b106ff6d68671f20dfe7d2c93241d3720
 burnman/chemicalpotentials.py      | 77 +++++++++++++++++++++++++++++++++-----
 burnman/test_chemicalpotentials.py | 34 ++++++++++++++---
 2 files changed, 95 insertions(+), 16 deletions(-)

diff --git a/burnman/chemicalpotentials.py b/burnman/chemicalpotentials.py
index 548681d..f2b5c23 100644
--- a/burnman/chemicalpotentials.py
+++ b/burnman/chemicalpotentials.py
@@ -6,28 +6,85 @@ import burnman
 from burnman import minerals
 from processchemistry import *
 import numpy as np
+from scipy.linalg import lu
 
+R=8.3145
 # This module computes chemical potentials (partial molar gibbs free energies) for an assemblage based on the Gibbs free energies and compositions of the individual phases.
 
 # It can also calculate fugacities based on the gibbs free energies of the endmembers corresponding to chemical components.
 
 # It can also calculate fugacities relative to other bulk compositions
 
-def chemicalpotentials(minerals, component_formulae):
-    mineral_gibbs=[]
-    mineral_formulae=[]
-    for mineral in minerals:
-        mineral_gibbs.append(mineral.gibbs)
-        mineral_formulae.append(mineral.params['formula'])
+def chemicalpotentials(assemblage, component_formulae):
+    """
+    The compositional space of the components does not have to be a 
+    superset of the compositional space of the assemblage. Nor do they have to
+    compose an orthogonal basis. 
 
+    The components must each be described by a linear mineral combination
+
+    The mineral compositions must be linearly independent
+
+        Parameters
+        ----------
+        assemblage : list of classes
+            List of material classes
+            set_method and set_state should already have been used
+
+        component_formulae : list of dictionaries
+            List of chemical component formula dictionaries
+            No restriction on length
+
+        Returns
+        -------
+        component_potentials : array of floats
+            Array of chemical potentials of components
+
+    """
+    mineral_gibbs=[mineral.gibbs for mineral in assemblage]
+    mineral_formulae=[mineral.params['formula'] for mineral in assemblage]
     mineral_compositions, elements=compositional_array(mineral_formulae)
 
-    for i in range(len(component_formulae)):
-        component_formulae[i]=dictionarize_formula(component_formulae[i])
+    pl, u = lu(mineral_compositions, permute_l=True)
+    assert(min(np.dot(np.square(u), np.ones(shape=len(elements)))) > 1e-6), \
+        'Mineral compositions do not form an independent set of basis vectors'
 
     component_compositions=ordered_compositional_array(component_formulae, elements)
-    component_proportions=np.linalg.lstsq(component_compositions.T, mineral_compositions.T)[0].T
-    component_potentials=np.linalg.solve(component_proportions,mineral_gibbs)
+
+    p=np.linalg.lstsq(mineral_compositions.T,component_compositions.T)
+    for idx, error in enumerate(p[1]):
+        assert (error < 1e-10), \
+            'Component %d not defined by prescribed assemblage' % (idx+1)
+
+    mineral_proportions=np.around(p[0],10).T
+
+    component_potentials=np.dot(mineral_proportions, mineral_gibbs)
 
     return component_potentials
 
+def fugacity(component_formula, standard_material, assemblage):
+    """
+        Parameters
+        ----------
+        component_formula : dictionary
+            Chemical formula dictionary
+
+        standard_material: class
+            Material class
+            set_method and set_state should already have been used
+
+        assemblage: list of classes
+            List of material classes
+            set_method and set_state should already have been used
+
+        Returns
+        -------
+        fugacity : float
+            Value of the fugacity of the component with respect to
+            the standard material
+
+    """
+    chemical_potential=chemicalpotentials(assemblage, [component_formula])[0]
+
+    fugacity=np.exp((chemical_potential - standard_material.gibbs)/(R*assemblage[0].temperature))
+    return fugacity
diff --git a/burnman/test_chemicalpotentials.py b/burnman/test_chemicalpotentials.py
index 1664d28..eb0b3ad 100644
--- a/burnman/test_chemicalpotentials.py
+++ b/burnman/test_chemicalpotentials.py
@@ -7,24 +7,46 @@ from burnman import minerals
 from chemicalpotentials import *
 import numpy as np
 
-
+'''
+Here we initialise the minerals we'll be using
+'''
 P=1.e8
 T=1000.
 
 py=minerals.HP_2011_ds62.py()
 fo=minerals.HP_2011_ds62.fo()
 en=minerals.HP_2011_ds62.en()
-
-minerals=[py, fo, en]
+per=minerals.HP_2011_ds62.per()
+minerals=[py, fo, en, per]
 
 for mineral in minerals:
     mineral.set_method('mtait')
     mineral.set_state(P, T)
 
+
+'''
+Here we find chemical potentials of MgO, Al2O3 and SiO2 for
+an assemblage containing pyrope, forsterite and enstatite 
+at 0.1 GPa, 1000 K
+'''
+
+assemblage=[py, fo, en]
 component_formulae=['MgO', 'Al2O3', 'SiO2']
+component_formulae_dict=[dictionarize_formula(f) for f in component_formulae]
+chem_potentials=chemicalpotentials(assemblage, component_formulae_dict)
+
+for idx, component_formula in enumerate(component_formulae):
+    print "mu("+component_formula+"):", chem_potentials[idx]
+print ''
 
-chem_potentials=chemicalpotentials(minerals, component_formulae)
 
-for idx, mineral in enumerate(minerals):
-    print mineral.params['name'], chem_potentials[idx]
+'''
+Here we find the 'periclase fugacity' of the assemblage
+Fugacity is often defined relative to some reference pressure
+Here we use room pressure, 100 kPa
+'''
 
+Pr=1.e5
+per.set_state(Pr, T)
+component=dictionarize_formula('MgO')
+print 'log(fugacity_per):', np.log10(fugacity(component, per, assemblage))