[cig-commits] [commit] master: Add order-disorder solid solution test, fixed 0K solid solution behaviour (8e75822)

[cig_noreply at geodynamics.org]

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

On branch  : master
Link       : https://github.com/geodynamics/burnman/compare/01d4aa027e2cb33852b7efdd80e67d64e5ec8f9c...881b878da7bd7228b4e714d9d10c4a75c876efb6

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

commit 8e75822469df55461ff056a0afae436faba077c5
Author: Bob Myhill <myhill.bob at gmail.com>
Date:   Fri Dec 12 19:43:01 2014 -0800

    Add order-disorder solid solution test, fixed 0K solid solution behaviour


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

8e75822469df55461ff056a0afae436faba077c5
 burnman/solidsolution.py                    | 10 +++++++---
 misc/benchmarks/solidsolution_benchmarks.py |  4 ++--
 tests/test_solidsolution.py                 | 20 ++++++++++++++++++++
 3 files changed, 29 insertions(+), 5 deletions(-)

diff --git a/burnman/solidsolution.py b/burnman/solidsolution.py
index 48c4951..7dcc85e 100644
--- a/burnman/solidsolution.py
+++ b/burnman/solidsolution.py
@@ -89,10 +89,14 @@ class SolidSolution(Mineral):
 
         # Derived properties
         self.C_v = self.C_p - self.V*temperature*self.alpha*self.alpha*self.K_T
-        self.K_S = self.K_T*self.C_p/self.C_v
-        self.gr = self.alpha*self.K_T*self.V/self.C_v
-
 
+        # C_v and C_p -> 0 as T -> 0
+        if temperature<1e-10:
+            self.K_S = self.K_T
+            self.gr = float('nan')
+        else:
+            self.K_S = self.K_T*self.C_p/self.C_v
+            self.gr = self.alpha*self.K_T*self.V/self.C_v     
 
     def calcgibbs(self, pressure, temperature, molar_fractions): 
         return sum([ self.base_material[i][0].calcgibbs(pressure, temperature) * molar_fractions[i] for i in range(self.n_endmembers) ]) + self.solution_model.excess_gibbs_free_energy( pressure, temperature, molar_fractions)
diff --git a/misc/benchmarks/solidsolution_benchmarks.py b/misc/benchmarks/solidsolution_benchmarks.py
index 9a6cb5f..21c2986 100644
--- a/misc/benchmarks/solidsolution_benchmarks.py
+++ b/misc/benchmarks/solidsolution_benchmarks.py
@@ -136,7 +136,7 @@ for idx, model in enumerate(opx_models):
     for i,c in enumerate(comp):
         molar_fractions=[1.0-c, c]
         model.set_composition( np.array(molar_fractions) )
-        model.set_state( 1.e5, 300. )
+        model.set_state( 0., 0. )
         opx_entropies[idx][i] = model.solution_model._configurational_entropy( molar_fractions )
 
 
@@ -180,7 +180,7 @@ gibbs = np.empty_like(comp)
 
 for i,c in enumerate(comp):
    cpx.set_composition( np.array([1.0-c, c]) )
-   cpx.set_state( 1.e5, 300. )
+   cpx.set_state( 0., 0. )
    gibbs[i] = cpx.excess_gibbs
 
 
diff --git a/tests/test_solidsolution.py b/tests/test_solidsolution.py
index be5ab2a..a123db9 100644
--- a/tests/test_solidsolution.py
+++ b/tests/test_solidsolution.py
@@ -91,6 +91,20 @@ class olivine_ss(burnman.SolidSolution):
         burnman.SolidSolution.__init__(self, base_material, \
                           burnman.solutionmodel.SymmetricRegularSolution(base_material, enthalpy_interaction) )
 
+# Orthopyroxene solid solution
+class orthopyroxene(burnman.SolidSolution):
+    def __init__(self):
+        # Name
+        self.name='orthopyroxene'
+
+        base_material = [[forsterite(), 'Mg[Mg]Si2O6'], [forsterite(), '[Mg1/2Al1/2]2AlSiO6']]
+
+        # Interaction parameters
+        enthalpy_interaction=[[10.0e3]]
+
+        burnman.SolidSolution.__init__(self, base_material, \
+                          burnman.solutionmodel.SymmetricRegularSolution(base_material, enthalpy_interaction) )
+
 
 
 
@@ -155,5 +169,11 @@ class test_solidsolution(BurnManTest):
         Wh=ol_ss.solution_model.Wh[0][1]
         self.assertArraysAlmostEqual([Wh/4.0], [H_excess])
 
+    def test_order_disorder(self):
+        opx = orthopyroxene()
+        opx.set_composition( np.array([0.0, 1.0]) )
+        self.assertArraysAlmostEqual([opx.excess_gibbs], [0.])
+
+
 if __name__ == '__main__':
     unittest.main()