[cig-commits] [commit] add_thermodynamic_potentials: More disorganized work (cc6b47d)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Dec 9 09:56:34 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_thermodynamic_potentials
Link : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51
>---------------------------------------------------------------
commit cc6b47d6957bec5eae9c3c1340bded52c57a182e
Author: ian-r-rose <ian.r.rose at gmail.com>
Date: Thu Sep 25 18:51:57 2014 -0700
More disorganized work
>---------------------------------------------------------------
cc6b47d6957bec5eae9c3c1340bded52c57a182e
burnman/composite.py | 1 +
burnman/equilibriumassemblage.py | 57 ++++++++--------------------
test_gibbs.py | 80 ++++++++++++++++++++++++++++++++--------
3 files changed, 81 insertions(+), 57 deletions(-)
diff --git a/burnman/composite.py b/burnman/composite.py
index 9c5c109..431a7c1 100644
--- a/burnman/composite.py
+++ b/burnman/composite.py
@@ -105,6 +105,7 @@ class Composite(Material):
self.temperature = temperature
for (fraction, phase) in self.children:
phase.set_state(pressure, temperature)
+ self.gibbs = np.sum( np.array([ ph.gibbs*frac for (frac, ph) in self.children]) )
def density(self):
"""
diff --git a/burnman/equilibriumassemblage.py b/burnman/equilibriumassemblage.py
index bfb8a84..64504fe 100644
--- a/burnman/equilibriumassemblage.py
+++ b/burnman/equilibriumassemblage.py
@@ -62,7 +62,6 @@ class EquilibriumAssemblage(burnman.Material):
self.__setup_subspaces()
self.__compute_baseline_assemblage()
-
def set_method(self, method):
for phase in self.phases:
phase.set_method(method)
@@ -70,15 +69,13 @@ class EquilibriumAssemblage(burnman.Material):
def set_state( self, pressure, temperature):
- self.__compute_baseline_assemblage()
n = len(self.endmember_formulae)
- minimize_gibbs = lambda x : self.__compute_gibbs( pressure, temperature, x )
- non_negative_constraint = lambda x : self.__species_vector( x )
- sol = opt.fmin( minimize_gibbs, self.reduced_species_vector, full_output=1, retall=1)
+ 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)
self.reduced_species_vector = sol[0]
- self.gibbs = sol[1]
self.species_vector = self.__species_vector(self.reduced_species_vector)
- self.print_assemblage()
+ self.gibbs = sol[1]*np.abs(ref_gibbs) + ref_gibbs
def __compute_gibbs( self, pressure, temperature, reduced_vector ):
@@ -91,21 +88,18 @@ class EquilibriumAssemblage(burnman.Material):
for phase in self.phases:
if isinstance (phase, burnman.SolidSolution):
n = len(phase.base_material)
- total_frac = np.sum( species_vector[i:(i+n)] )/np.sum(species_vector)
phase.set_method('slb3')
phase.set_composition( np.array( species_vector[i:(i+n)]/np.sum(species_vector[i:(i+n)])) )
phase.set_state( pressure, temperature )
- tmp_gibbs += phase.gibbs * total_frac
+ tmp_gibbs += phase.gibbs
i+=n
elif isinstance(phase, burnman.Mineral):
phase.set_method('slb3')
phase.set_state( pressure, temperature )
- tmp_gibbs += phase.gibbs * species_vector[i]/np.sum(species_vector)
+ tmp_gibbs += phase.gibbs * species_vector[i]
i+=1
else:
raise Exception('Unsupported mineral type, can only read burnman.Mineral or burnman.SolidSolution')
- if np.any(species_vector < -1.e-6):
- return 1.e30
return tmp_gibbs
@@ -128,8 +122,8 @@ class EquilibriumAssemblage(burnman.Material):
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.sparsify_basis(self.right_nullspace)
- for i in range( self.right_nullspace.shape[1]):
- print self.right_nullspace[:,i]/np.max(np.abs(self.right_nullspace[:,i]))
+ 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)
@@ -139,40 +133,19 @@ class EquilibriumAssemblage(burnman.Material):
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 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.
+ #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)
- 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)
-
- baseline_assemblage = np.dot(linalg.pinv2( self.stoichiometric_matrix ) , self.bulk_composition_vector)
-
- i = 0
- while i < len( baseline_assemblage) :
- if baseline_assemblage[i] < 0.0:
- print baseline_assemblage
- react_id = np.argmax( self.right_nullspace[i,:])
- reaction = self.right_nullspace[:,react_id]
- baseline_assemblage = baseline_assemblage - reaction * baseline_assemblage[i]/reaction[i]
- print baseline_assemblage
- i = 0
- elif baseline_assemblage[i] < eps:
- baseline_assemblage[i] = 0.0
- i = i+1
- else:
- i = i+1
+ 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." )
- assert( np.all(np.abs(np.dot(self.stoichiometric_matrix, baseline_assemblage)\
+ assert( np.all(np.abs(np.dot(self.stoichiometric_matrix, baseline_assemblage[0])\
- self.bulk_composition_vector) < eps) )
- self.baseline_assemblage = baseline_assemblage
+ 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)
diff --git a/test_gibbs.py b/test_gibbs.py
index 9de25d4..424b37d 100644
--- a/test_gibbs.py
+++ b/test_gibbs.py
@@ -4,36 +4,86 @@ import numpy as np
import matplotlib.pyplot as plt
-#minlist = [mg_fe_olivine(), mg_fe_wadsleyite()]
-#composition = {'Fe': 1./35., 'Mg': 9./35., 'O': 4./7., 'Si': 1./7}
composition = { 'Mg': 2./7., 'O': 4./7., 'Si': 1./7}
-minlist = [forsterite(), mg_wadsleyite(), mg_ringwoodite(), mg_perovskite(), periclase(), stishovite()]
-#minlist = [forsterite(), mg_wadsleyite(), mg_ringwoodite()]
+minlist = [forsterite(), mg_wadsleyite(), mg_ringwoodite(), mg_perovskite(), periclase(), stishovite(),mg_akimotoite(), mg_majorite() ]
fo = forsterite()
fo.set_method('slb3')
wa = mg_wadsleyite()
wa.set_method('slb3')
-
-
+ri = mg_ringwoodite()
+ri.set_method('slb3')
+pvpe = burnman.Composite( [0.5, 0.5], [mg_perovskite(), periclase()])
+pvpe.set_method('slb3')
+pest = burnman.Composite( [0.5, 0.5], [periclase(), stishovite()])
+pest.set_method('slb3')
ea = burnman.EquilibriumAssemblage(composition, minlist)
ea.set_method('slb3')
-n= 60.
+n= 100.
pressures = np.linspace(1.e5, 30.e9, n)
temperatures = np.linspace(800., 2800.0, n)
+G = np.empty_like(pressures)
+Gr = np.empty_like(pressures)
+Gf = np.empty_like(pressures)
+Gw = np.empty_like(pressures)
+Gp = np.empty_like(pressures)
+Gs = np.empty_like(pressures)
+Ga = np.empty_like(pressures)
phase = np.empty( (n,n) )
-eps = 1.e-5
+T = 1600.
+
+for i,p in enumerate(pressures):
+ ea.set_state(p, T)
+ G[i] = ea.gibbs*7.
+ fo.set_state(p, T)
+ Gf[i] = fo.gibbs
+ wa.set_state(p, T)
+ Gw[i] = wa.gibbs
+ ri.set_state(p, T)
+ Gr[i] = ri.gibbs
+ pvpe.set_state(p, T)
+ Gp[i] = pvpe.gibbs*2.
+Ga = (Gw+Gr+Gf+Gp)*1./4.
+plt.plot(pressures, G-Ga, 'k--', linewidth=4)
+plt.plot(pressures, Gr-Ga)
+plt.plot(pressures, Gw-Ga)
+plt.plot(pressures, Gf-Ga)
+plt.plot(pressures, Gp-Ga)
+plt.show()
+plt.clf()
-#for i,p in enumerate(pressures):
-# for j,t in enumerate(temperatures):
-# ea.set_state(p, t)
-# phase[j,i] = np.argmax(ea.species_vector)
-#plt.imshow( np.flipud(phase))
-#plt.colorbar()
-#plt.show()
+eps = 1.e-4
+for i,p in enumerate(pressures):
+ for j,t in enumerate(temperatures):
+ ea.set_state(p, t)
+ #forsterite field
+ if ea.species_vector[0]/np.sum(ea.species_vector) > eps:
+ phase[j,i] = 0.
+ #wadsleyite field
+ if ea.species_vector[1]/np.sum(ea.species_vector) > eps:
+ phase[j,i] = 1.
+ #ringwoodite field
+ if ea.species_vector[2]/np.sum(ea.species_vector) > eps:
+ phase[j,i] = 2.
+ #perovskite+periclase field
+ if ea.species_vector[3]/np.sum(ea.species_vector) > eps and ea.species_vector[4]/np.sum(ea.species_vector) > eps:
+ phase[j,i] = 3.
+ #periclase + stishovite field
+ if ea.species_vector[4]/np.sum(ea.species_vector) > eps and ea.species_vector[5]/np.sum(ea.species_vector) > eps:
+ phase[j,i] = 4.
+ #akimotoite + periclase field
+ if ea.species_vector[4]/np.sum(ea.species_vector) > eps and ea.species_vector[6]/np.sum(ea.species_vector) > eps:
+ phase[j,i] = 5.
+ #majorite + periclase field
+ if ea.species_vector[4]/np.sum(ea.species_vector) > eps and ea.species_vector[7]/np.sum(ea.species_vector) > eps:
+ phase[j,i] = 6.
+plt.imshow( phase, origin='lower', extent=[0., 30., 800., 2800.], aspect='auto')
+plt.xlabel("Pressure (GPa)")
+plt.ylabel("Temperature (K)")
+plt.show()
More information about the CIG-COMMITS
mailing list