[cig-commits] [commit] add_thermodynamic_potentials: Cleanup. (4f5cce6)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Dec 9 09:56:38 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_thermodynamic_potentials
Link : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51
>---------------------------------------------------------------
commit 4f5cce60c8214932785c93c27dbf4d869bc38ac6
Author: ian-r-rose <ian.r.rose at gmail.com>
Date: Fri Sep 26 19:54:56 2014 -0700
Cleanup.
>---------------------------------------------------------------
4f5cce60c8214932785c93c27dbf4d869bc38ac6
burnman/equilibriumassemblage.py | 52 +++++++++++++++-------------------------
burnman/gibbsminimization.py | 3 +--
2 files changed, 20 insertions(+), 35 deletions(-)
diff --git a/burnman/equilibriumassemblage.py b/burnman/equilibriumassemblage.py
index 65b2849..63a2762 100644
--- a/burnman/equilibriumassemblage.py
+++ b/burnman/equilibriumassemblage.py
@@ -17,7 +17,7 @@ class EquilibriumAssemblage(burnman.Material):
Class for taking a certain assemblage of elements,
then calculating the equilibrium assemblage of minerals
that minimizes the Gibbs free energy. This should
- have similar functionality to ``burnman.Composite'',
+ have similar functionality to :class:`burnman.Composite`,
but instead of having a fixed set of minerals, it
dynamically updates the minerals when you call set_state().
As such, it should be significantly slower to use.
@@ -31,7 +31,7 @@ class EquilibriumAssemblage(burnman.Material):
Parameters
----------
- elements : dictionary
+ composition : dictionary
Dictionary where the keys are strings corresponding to
elements (e.g. 'Mg') and the values are molar fractions
of that element.
@@ -56,12 +56,10 @@ 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()
self.__compute_baseline_assemblage()
+
def set_method(self, method):
for phase in self.phases:
phase.set_method(method)
@@ -72,14 +70,21 @@ class EquilibriumAssemblage(burnman.Material):
n = len(self.endmember_formulae)
ref_gibbs = self.__compute_gibbs( pressure, temperature, self.reduced_species_vector *0.)
minimize_gibbs = lambda x : (self.__compute_gibbs( pressure, temperature, x ) - ref_gibbs)/np.abs(ref_gibbs)
- sol = opt.fmin_slsqp( minimize_gibbs, self.reduced_species_vector, f_ieqcons = lambda x : self.__species_vector(x), full_output=1, disp=False)
+ sol = opt.fmin_slsqp( minimize_gibbs, self.reduced_species_vector, f_ieqcons = lambda x : self.__compute_species_vector(x), full_output=1, disp=False)
self.reduced_species_vector = sol[0]
- self.species_vector = self.__species_vector(self.reduced_species_vector)
+ self.species_vector = self.__compute_species_vector(self.reduced_species_vector)
self.gibbs = sol[1]*np.abs(ref_gibbs) + ref_gibbs
+
+ def print_assemblage(self):
+ tot = np.sum(self.species_vector)
+ for f,s in zip(self.endmember_formulae, self.species_vector):
+ print f, s/tot
+
+
def __compute_gibbs( self, pressure, temperature, reduced_vector ):
- species_vector = self.__species_vector(reduced_vector)
+ species_vector = self.__compute_species_vector(reduced_vector)
assert( len(species_vector) == len(self.endmember_formulae) )
tmp_gibbs = 0.0
@@ -105,39 +110,24 @@ class EquilibriumAssemblage(burnman.Material):
return tmp_gibbs
- def __species_vector ( self, reduced_vector ):
+
+ def __compute_species_vector ( self, reduced_vector ):
return self.baseline_assemblage + np.dot( self.right_nullspace, np.transpose(reduced_vector) )
- def __bulk_composition (self, species_vector):
+
+ def __compute_bulk_composition (self, species_vector):
return np.dot(self.stoichiometric_matrix, species_vector)
- def print_assemblage(self):
- tot = np.sum(self.species_vector)
- for f,s in zip(self.endmember_formulae, self.species_vector):
- print f, s/tot
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))
+ self.right_nullspace = gm.compute_nullspace( self.stoichiometric_matrix )
self.right_nullspace = gm.sparsify_basis(self.right_nullspace)
- for col in range(self.right_nullspace.shape[1]):
- self.right_nullspace[:,col] = self.right_nullspace[:,col]/np.sum(np.abs(self.right_nullspace[:,col]))
- left_null_mask = ( S <= eps)
- self.left_nullspace = np.compress(left_null_mask, U, axis=1)
def __compute_baseline_assemblage(self):
eps = 1.e-10
- #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 throw an error if so.
- null_power = np.dot(self.left_nullspace.T, self.bulk_composition_vector)
-
baseline_assemblage = opt.nnls( self.stoichiometric_matrix, self.bulk_composition_vector)
if baseline_assemblage[1] > eps :
raise Exception( "Composition cannot be represented by the given minerals." )
@@ -147,8 +137,4 @@ class EquilibriumAssemblage(burnman.Material):
self.baseline_assemblage = baseline_assemblage[0]
self.reduced_species_vector = np.zeros( self.right_nullspace.shape[1] )
- self.species_vector = self.__species_vector( self.reduced_species_vector)
-
-
-
-
+ self.species_vector = self.__compute_species_vector( self.reduced_species_vector)
diff --git a/burnman/gibbsminimization.py b/burnman/gibbsminimization.py
index f3e4320..b3ceba4 100644
--- a/burnman/gibbsminimization.py
+++ b/burnman/gibbsminimization.py
@@ -105,8 +105,7 @@ def compute_nullspace ( stoichiometric_matrix ):
# The columns of V that correspond to small singular values
# (or do not exist) are the left nullspace for the matrix.
# select them, appropriately transpose them, and return them
- S=np.append(S, [0. for i in range(len(Vh)-len(S))])
- null_mask = (S <= eps)
+ null_mask = ( np.append(S, np.zeros(len(Vh)-len(S))) <= eps)
null_space = np.compress(null_mask, Vh, axis=0)
return np.transpose(null_space)
More information about the CIG-COMMITS
mailing list