[cig-commits] [commit] add_gibbs_energy: Added temporary test for previously buggy gibbs excesses (4ac0c70)

[cig_noreply at geodynamics.org]

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

On branch  : add_gibbs_energy
Link       : https://github.com/geodynamics/burnman/compare/0000000000000000000000000000000000000000...2148b324d3e8aa7b527f831eb397590942563008

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

commit 4ac0c706fcd53e6ea18cc0199918813cc0a9af10
Author: Bob Myhill <myhill.bob at gmail.com>
Date:   Sun Oct 5 23:13:57 2014 +0200

    Added temporary test for previously buggy gibbs excesses


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

4ac0c706fcd53e6ea18cc0199918813cc0a9af10
 solidsolution_benchmarks.py | 23 +++++++++++++++++++++--
 1 file changed, 21 insertions(+), 2 deletions(-)

diff --git a/solidsolution_benchmarks.py b/solidsolution_benchmarks.py
index 67eef9c..92fef5d 100644
--- a/solidsolution_benchmarks.py
+++ b/solidsolution_benchmarks.py
@@ -148,7 +148,7 @@ class orthopyroxene_long_dashed(burnman.SolidSolution):
         base_material = [[enstatite(), 'Mg[Mg]Si2O6'],[mg_tschermaks_molecule(), '[Mg1/2Al1/2]2AlSiO6'] ]
 
         # Interaction parameters
-        enthalpy_interaction=[[0.0]]
+        enthalpy_interaction=[[10.0e3]]
 
         burnman.SolidSolution.__init__(self, base_material, \
                           burnman.solutionmodel.SymmetricRegularSolution(base_material, enthalpy_interaction) )
@@ -280,4 +280,23 @@ plt.ylabel("Excess enthalpy of solution (kJ/mol)")
 plt.xlabel("cats fraction")
 plt.show()
 
-# Excess volume of solution
+# Check endmember excess Gibbs goes to zero... 
+
+opx = orthopyroxene_long_dashed()
+opx.set_method('slb3')
+
+comp = np.linspace(0, 1.0, 100)
+gibbs = np.empty_like(comp)
+
+for i,c in enumerate(comp):
+   opx.set_composition( np.array([1.0-c, c]) )
+   opx.set_state( 0.0, 2000.0 )
+   gibbs[i] = opx.excess_gibbs
+
+
+plt.plot( comp, gibbs/1000., 'b--', linewidth=3.)
+plt.xlim(0.0,1.0)
+plt.ylim(-2.,8.0)
+plt.ylabel("Excess enthalpy of solution (kJ/mol)")
+plt.xlabel("mgts fraction")
+plt.show()