[cig-commits] [commit] add_thermodynamic_potentials: More work on setting up the gibbs minimization problem (655901b)

[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 655901ba5b2b9e782e6d8926346836a506291847
Author: ian-r-rose <ian.r.rose at gmail.com>
Date:   Fri Sep 19 10:00:26 2014 -0700

    More work on setting up the gibbs minimization problem


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

655901ba5b2b9e782e6d8926346836a506291847
 burnman/equilibriumassemblage.py | 42 ++++++++++++++++++++++++++++++++++++++++
 burnman/gibbsminimization.py     |  6 +++---
 test_gibbs.py                    |  5 ++---
 3 files changed, 47 insertions(+), 6 deletions(-)

diff --git a/burnman/equilibriumassemblage.py b/burnman/equilibriumassemblage.py
index e71c3ea..779a82f 100644
--- a/burnman/equilibriumassemblage.py
+++ b/burnman/equilibriumassemblage.py
@@ -4,6 +4,7 @@
 
 import numpy as np
 import scipy.optimize as opt
+import scipy.linalg as linalg
 import warnings
 
 import burnman
@@ -54,6 +55,12 @@ class EquilibriumAssemblage(burnman.Material):
         self.bulk_composition_vector = np.array([ (composition[e] if e in composition.keys() else 0.0) \
                                                    for e in self.elements] )
  
+        self.nullspace = gm.compute_nullspace( self.stoichiometric_matrix )
+        self.pseudoinverse = np.dot(linalg.pinv2( self.stoichiometric_matrix ) , self.bulk_composition_vector)
+
+        self.__setup_subspaces()
+        
+
     def set_method(self, method):
         for phase in self.phases:
             phase.set_method(method)
@@ -91,6 +98,7 @@ class EquilibriumAssemblage(burnman.Material):
             else:
                 raise Exception('Unsupported mineral type, can only read burnman.Mineral or burnman.SolidSolution')
 
+        print tmp_gibbs
         return tmp_gibbs
                 
         
@@ -99,3 +107,37 @@ class EquilibriumAssemblage(burnman.Material):
         return np.dot( self.stoichiometric_matrix, species_vector) - self.bulk_composition_vector
 
 
+    def __setup_subspaces (self):
+ 
+        eps = 1.e-10
+        U, S, Vh = linalg.svd( self.stoichiometric_matrix)
+
+        right_null_mask = ( np.append(S, np.zeros(len(Vh)-len(S))) <= eps)
+        self.right_nullspace = np.transpose(np.compress(right_null_mask, Vh, axis=0))
+
+        left_null_mask = ( S <= eps)
+        self.left_nullspace = np.compress(left_null_mask, U, axis=1)
+ 
+        #It is possible to give a bag of elements that cannot be represented by the list of
+        #minerals.  This corresponds to the bulk_composition_vector having power in the left nullspace
+        #of the stoichiometric matrix.  Here we check for this, and if it is the case, project
+        #it out of the left nullspace.  For most mantle assemblages, this is probably due to the
+        #amount of oxygen given being inconsistent with the possible minerals.
+        null_power = np.dot(self.left_nullspace.T, self.bulk_composition_vector)
+        if np.any( null_power > eps ):
+            print "Composition cannot be represented by the given minerals. We are projecting the composition onto the closest we can do, but you should probably recheck the composition vector"
+            for col in range(self.left_nullspace.shape[1]):
+                self.bulk_composition_vector -= self.left_nullspace[:,col]*null_power[col]
+            self.bulk_composition_vector = self.bulk_composition_vector/sum(self.bulk_composition_vector)
+            print "New vector: ", zip(self.elements, self.bulk_composition_vector)
+        
+
+        self.baseline_assemblage = np.dot(linalg.pinv2( self.stoichiometric_matrix ) , self.bulk_composition_vector)
+        assert( np.all(np.abs(np.dot(self.stoichiometric_matrix, self.baseline_assemblage)\
+                                                     - self.bulk_composition_vector) < eps) )
+        print zip(self.endmember_formulae, self.baseline_assemblage)#/np.sum(self.baseline_assemblage))
+    
+  
+
+
+
diff --git a/burnman/gibbsminimization.py b/burnman/gibbsminimization.py
index b63a463..fe04bda 100644
--- a/burnman/gibbsminimization.py
+++ b/burnman/gibbsminimization.py
@@ -65,9 +65,9 @@ def compute_nullspace ( stoichiometric_matrix ):
     """
 
     eps = 1.e-10
-    U, S, V = linalg.svd( stoichiometric_matrix)
-    S=np.append(S, [0. for i in range(len(V)-len(S))])
+    U, S, Vh = linalg.svd( stoichiometric_matrix)
+    S=np.append(S, [0. for i in range(len(Vh)-len(S))])
     null_mask = (S <= eps)
-    null_space = np.compress(null_mask, V, axis=0)
+    null_space = np.compress(null_mask, Vh, axis=0)
     return np.transpose(null_space)
     
diff --git a/test_gibbs.py b/test_gibbs.py
index e231407..a2bd4f7 100644
--- a/test_gibbs.py
+++ b/test_gibbs.py
@@ -4,9 +4,8 @@ import numpy as np
 import matplotlib.pyplot as plt
 
 
-minlist = [mg_fe_olivine(), wuestite(), periclase(), stishovite(), jadeite()]
-composition = {'Mg': 0.2, 'Si': 0.2, 'Fe': 0.1, 'O': 0.5}
+minlist = [mg_fe_olivine(), ferropericlase(), stishovite()]
+composition = {'Fe': 1./7., 'Mg': 1./7., 'O': 4./7., 'Si': 1./7}
 
 ea = burnman.EquilibriumAssemblage(composition, minlist)
 ea.set_method('slb3')
-ea.set_state(100.e9, 3000.)