[cig-commits] [commit] add_thermodynamic_potentials: Updated chemical potential processing to accept solid solutions (2547bee)

[cig_noreply at geodynamics.org]

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

On branch  : add_thermodynamic_potentials
Link       : https://github.com/geodynamics/burnman/compare/f9a1e215ff875e89f69b5e717b85d9af5d44e93b...83c31a2149ddfb6f77c6a72fb5cf3eef08996a42

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

commit 2547bee82c10c9f43780a62d2afb28efb9cdb8fd
Author: Bob Myhill <myhill.bob at gmail.com>
Date:   Tue Nov 11 22:06:06 2014 +0100

    Updated chemical potential processing to accept solid solutions


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

2547bee82c10c9f43780a62d2afb28efb9cdb8fd
 burnman/chemicalpotentials.py | 37 +++++++++++++++++++++++++++----------
 burnman/solidsolution.py      |  1 -
 2 files changed, 27 insertions(+), 11 deletions(-)

diff --git a/burnman/chemicalpotentials.py b/burnman/chemicalpotentials.py
index fb57d96..89ab1fc 100644
--- a/burnman/chemicalpotentials.py
+++ b/burnman/chemicalpotentials.py
@@ -30,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
@@ -41,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):
diff --git a/burnman/solidsolution.py b/burnman/solidsolution.py
index 6297db6..80c864f 100644
--- a/burnman/solidsolution.py
+++ b/burnman/solidsolution.py
@@ -60,7 +60,6 @@ class SolidSolution(Mineral):
 
         self.V= sum([ self.base_material[i][0].V * self.molar_fraction[i] for i in range(self.n_endmembers) ]) + self.excess_volume
         
-        
         '''
         for prop in self.base_materials[0].params:
            try: