[cig-commits] [commit] add_thermodynamic_potentials: Added chemicalpotentials stub (9ac0935)

[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 9ac0935c40a10c03e56753c81ba1958aa035214c
Author: Bob Myhill <myhill.bob at gmail.com>
Date:   Wed Sep 3 03:12:21 2014 +0200

    Added chemicalpotentials stub


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

9ac0935c40a10c03e56753c81ba1958aa035214c
 burnman/chemicalpotentials.py      | 33 +++++++++++++++++++++++++++++++++
 burnman/test_chemicalpotentials.py | 30 ++++++++++++++++++++++++++++++
 2 files changed, 63 insertions(+)

diff --git a/burnman/chemicalpotentials.py b/burnman/chemicalpotentials.py
new file mode 100644
index 0000000..548681d
--- /dev/null
+++ b/burnman/chemicalpotentials.py
@@ -0,0 +1,33 @@
+# 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 burnman
+from burnman import minerals
+from processchemistry import *
+import numpy as np
+
+# 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'])
+        
+    mineral_compositions, elements=compositional_array(mineral_formulae)
+
+    for i in range(len(component_formulae)):
+        component_formulae[i]=dictionarize_formula(component_formulae[i])
+
+    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)
+
+    return component_potentials
+
diff --git a/burnman/test_chemicalpotentials.py b/burnman/test_chemicalpotentials.py
new file mode 100644
index 0000000..1664d28
--- /dev/null
+++ b/burnman/test_chemicalpotentials.py
@@ -0,0 +1,30 @@
+import os, sys, numpy as np, matplotlib.pyplot as plt
+if not os.path.exists('burnman') and os.path.exists('../burnman'):
+    sys.path.insert(1,os.path.abspath('..'))
+
+import burnman
+from burnman import minerals
+from chemicalpotentials import *
+import numpy as np
+
+
+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]
+
+for mineral in minerals:
+    mineral.set_method('mtait')
+    mineral.set_state(P, T)
+
+component_formulae=['MgO', 'Al2O3', 'SiO2']
+
+chem_potentials=chemicalpotentials(minerals, component_formulae)
+
+for idx, mineral in enumerate(minerals):
+    print mineral.params['name'], chem_potentials[idx]
+