[cig-commits] [commit] add_thermodynamic_potentials: Expanded use of chemicalpotentials to solid solutions (f9b9521)

[cig_noreply at geodynamics.org]

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

On branch  : add_thermodynamic_potentials
Link       : https://github.com/geodynamics/burnman/compare/d615573d0c849392955acc74130c13b6109784e2...f9b95219360a9197b16a07248236bbdabe05485c

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

commit f9b95219360a9197b16a07248236bbdabe05485c
Author: Bob Myhill <myhill.bob at gmail.com>
Date:   Tue Dec 9 10:16:39 2014 -0800

    Expanded use of chemicalpotentials to solid solutions


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

f9b95219360a9197b16a07248236bbdabe05485c
 burnman/chemicalpotentials.py | 41 ++++++++++++++++++++++++++++-------------
 1 file changed, 28 insertions(+), 13 deletions(-)

diff --git a/burnman/chemicalpotentials.py b/burnman/chemicalpotentials.py
index b4d089e..8b054dc 100644
--- a/burnman/chemicalpotentials.py
+++ b/burnman/chemicalpotentials.py
@@ -7,9 +7,7 @@ from burnman import minerals
 from processchemistry import *
 import numpy as np
 from scipy.linalg import lu
-from burnman.constants import gas_constant
-
-R=gas_constant
+from burnman.constants import R
 
 # 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.
 
@@ -32,6 +30,7 @@ def chemicalpotentials(assemblage, component_formulae):
         assemblage : list of classes
             List of material classes
             set_method and set_state should already have been used
+            the composition of the solid solutions should also have been set
 
         component_formulae : list of dictionaries
             List of chemical component formula dictionaries
@@ -43,25 +42,41 @@ def chemicalpotentials(assemblage, component_formulae):
             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)
-
-    pl, u = lu(mineral_compositions, permute_l=True)
+    # Split solid solutions into their respective endmembers
+    # Find the chemical potentials of all the endmembers
+    endmember_list=[]
+    endmember_potentials=[]
+    for mineral in assemblage:
+        if isinstance(mineral, burnman.SolidSolution):
+            for member in mineral.base_material:
+                endmember_list.append(member[0])
+            for potential in mineral.partial_gibbs:
+                endmember_potentials.append(potential)
+        else:
+            endmember_list.append(mineral)
+            endmember_potentials.append(mineral.gibbs)
+   
+    # Make an array of all the endmember formulae
+    endmember_formulae=[endmember.params['formula'] for endmember in endmember_list]
+    endmember_compositions, elements=compositional_array(endmember_formulae)
+
+    pl, u = lu(endmember_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'
+        'Endmember compositions do not form an independent set of basis vectors'
 
+    # Make an array of component formulae with elements in the same order as the endmember array
     component_compositions=ordered_compositional_array(component_formulae, elements)
 
-    p=np.linalg.lstsq(mineral_compositions.T,component_compositions.T)
+    p=np.linalg.lstsq(endmember_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)
+    # Create an array of endmember proportions which sum to each component composition
+    endmember_proportions=np.around(p[0],10).T
 
+    # Calculate the chemical potential of each component
+    component_potentials=np.dot(endmember_proportions, endmember_potentials)
     return component_potentials
 
 def fugacity(component_formula, standard_material, assemblage):