[cig-commits] [commit] add_thermodynamic_potentials: Add test script. Not working (57e145b)

[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 57e145b3167487ea9e17073ae7dcef3c68156d91
Author: ian-r-rose <ian.r.rose at gmail.com>
Date:   Fri Oct 3 10:01:15 2014 -0700

    Add test script. Not working


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

57e145b3167487ea9e17073ae7dcef3c68156d91
 olivine_wadsleyite.py | 42 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 42 insertions(+)

diff --git a/olivine_wadsleyite.py b/olivine_wadsleyite.py
new file mode 100644
index 0000000..4415ec1
--- /dev/null
+++ b/olivine_wadsleyite.py
@@ -0,0 +1,42 @@
+import burnman
+from burnman.minerals.SLB_2011 import *
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+composition = { 'Fe': 1./35., 'Mg': 9./35., 'O': 4./7., 'Si': 1./7.}
+minlist = [mg_fe_olivine(), mg_fe_wadsleyite()]
+
+ol = mg_fe_olivine()
+ol.set_method('slb3')
+ol.set_composition( [0.9, 0.1])
+wa = mg_fe_wadsleyite()
+wa.set_method('slb3')
+wa.set_composition( [0.9, 0.1])
+ea = burnman.EquilibriumAssemblage(composition, minlist)
+ea.set_method('slb3')
+co = burnman.Composite( [0.5, 0.5], [ol, wa] )
+
+pressures = np.linspace( 1.e9, 30.e9, 100 )
+temperature = 1500.
+Gol = np.empty_like(pressures)
+Gwa = np.empty_like(pressures)
+Gea = np.empty_like(pressures)
+Gav = np.empty_like(pressures)
+
+
+for i,p in enumerate(pressures):
+    ol.set_state(p,temperature)
+    wa.set_state(p,temperature)
+    ea.set_state(p,temperature)
+
+    Gol[i] = ol.gibbs
+    Gwa[i] = wa.gibbs
+    Gea[i] = ea.gibbs*7.
+
+
+Gav = (Gol+Gwa)*1./2.
+plt.plot(pressures, Gol-Gav)
+plt.plot(pressures, Gwa-Gav)
+plt.plot(pressures, Gea-Gav, 'k--', linewidth=4)
+plt.show()