[cig-commits] [commit] add_thermodynamic_potentials: Added alpha version of tc-like equilibrium assemblage calculations (4451e82)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Dec 9 09:57:24 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_thermodynamic_potentials
Link : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51
>---------------------------------------------------------------
commit 4451e82be12799c32cd0edbf2af0030efd59d03b
Author: Bob Myhill <myhill.bob at gmail.com>
Date: Sun Oct 5 23:16:40 2014 +0200
Added alpha version of tc-like equilibrium assemblage calculations
>---------------------------------------------------------------
4451e82be12799c32cd0edbf2af0030efd59d03b
burnman/equilibriumequations.py | 60 ++++++++++++++++++
burnman/test_tc_clone.py | 133 +++++++++++++++++++++++++++++++++-------
2 files changed, 170 insertions(+), 23 deletions(-)
diff --git a/burnman/equilibriumequations.py b/burnman/equilibriumequations.py
new file mode 100644
index 0000000..1318777
--- /dev/null
+++ b/burnman/equilibriumequations.py
@@ -0,0 +1,60 @@
+import os, sys, numpy as np, matplotlib.pyplot as plt
+if not os.path.exists('burnman') and os.path.exists('../burnman'):
+ sys.path.insert(1,os.path.abspath('..'))
+
+import burnman
+import numpy as np
+
+
+def set_eqns(arg, bulk_composition, mineral_list, stoic, null, fixed_vars):
+ P=float(arg[0])
+ T=float(arg[1]) # arg[0] and arg[1] are P, T
+
+ mbridx=0
+ minidx=0
+ # mineral proportions and solid solution compositions
+ propscomps=np.array(arg[2:])
+ partialgibbs=np.empty(len(propscomps))
+ endmember_proportions=np.empty(len(propscomps))
+
+ # The following for loop
+ # finds the partial Gibbs free energies of the endmembers
+ for mineral in mineral_list:
+ mineral_proportion=propscomps[mbridx]
+ if isinstance(mineral, burnman.SolidSolution):
+ mbrprops=np.zeros(len(mineral.base_material))
+ mbrprops[0:len(mineral.base_material)-1]=propscomps[mbridx+1:mbridx+len(mineral.base_material)]
+ mbrprops[len(mineral.base_material)-1]=1.0-sum(mbrprops)
+ mineral.set_composition(mbrprops)
+ partial_gibbs_excesses=mineral.calcpartialgibbsexcesses(P, T, mbrprops)
+ for idx, endmember in enumerate(mineral.base_material):
+ endmember_proportions[minidx]=mineral_proportion*mbrprops[idx]
+ minidx=minidx+1
+ partialgibbs[mbridx]=endmember[0].calcgibbs(P, T) + partial_gibbs_excesses[idx]
+ mbridx=mbridx+1
+ else:
+ endmember_proportions[minidx]=mineral_proportion
+ minidx=minidx+1
+ partialgibbs[mbridx]=mineral.calcgibbs(P, T)
+ mbridx=mbridx+1
+
+ #### ADD EQUATIONS TO SOLVE ####
+ eqns=[]
+
+ # Add independent endmember reactions
+ gibbs_inequalities=np.dot(partialgibbs,null)
+ for ieq in gibbs_inequalities:
+ eqns.append(ieq)
+
+ # Add fixed variable constraints
+ eqns.append(arg[fixed_vars[0][0]]-fixed_vars[0][1])
+ eqns.append(arg[fixed_vars[1][0]]-fixed_vars[1][1])
+
+ # Add bulk compositional constraints
+ mineral_bulk_composition=np.dot(endmember_proportions,stoic[0].T)
+
+# for i in range(len(mineral_bulk_composition)):
+ for i in range(2):
+ eqns.append(mineral_bulk_composition[i] - bulk_composition[i])
+
+ return eqns
diff --git a/burnman/test_tc_clone.py b/burnman/test_tc_clone.py
index 56e0d9d..e9e905d 100644
--- a/burnman/test_tc_clone.py
+++ b/burnman/test_tc_clone.py
@@ -7,19 +7,47 @@ from burnman.minerals.SLB_2011 import *
from burnman.gibbsminimization import *
import numpy as np
from scipy.linalg import lstsq
+import scipy.optimize as opt
import matplotlib.pyplot as plt
+from burnman.equilibriumequations import *
+composition = { 'Mg': 1.8, 'Fe': 0.2, 'O': 4., 'Si': 1.}
-composition = { 'Mg': 1.8/7., 'Fe': 0.2/7., 'O': 4./7., 'Si': 1./7}
+mineral_list = [mg_fe_olivine(), mg_fe_wadsleyite()]
+for mineral in mineral_list:
+ mineral.set_method('slb3')
+ if isinstance(mineral, burnman.SolidSolution):
+ molar_fractions=np.zeros(len(mineral.base_material))
+ molar_fractions[0]=1.0
+ mineral.set_composition(molar_fractions)
+ mineral.set_state(0.,0.)
-minlist = [mg_fe_olivine(), mg_fe_wadsleyite()]
+'''
+molar_fractions=np.array([0.9,0.1])
+ol=mg_fe_olivine()
+fo=forsterite()
+fa=fayalite()
+
+ol.set_composition(molar_fractions)
+ol.set_state(P, T)
+fo.set_method('slb3')
+fa.set_method('slb3')
+fo.set_state(P,T)
+fa.set_state(P,T)
+
+# Test excess function
+print ol.gibbs
+print molar_fractions[0]*fo.gibbs + molar_fractions[1]*fa.gibbs + ol.excess_gibbs
+
+# Test partial gibbs functions
+print ol.excess_gibbs
+print 0.9*ol.calcpartialgibbsexcesses(P,T,molar_fractions)[0] + 0.1*ol.calcpartialgibbsexcesses(P,T,molar_fractions)[1]
+'''
-stoic=assemble_stoichiometric_matrix ( minlist )
+stoic=assemble_stoichiometric_matrix ( mineral_list )
null=sparsify_basis(compute_nullspace(stoic[0]))
-print null
-
# Create a compositional vector with elements in the same order as the
-# stoiciometric matrix
+# stoichiometric matrix
comp_vector=np.empty( (len(stoic[1])) )
for idx, element in enumerate(stoic[1]):
comp_vector[idx]=composition[element]
@@ -31,6 +59,7 @@ resid = np.dot(stoic[0],x) - comp_vector
tol=1e-10
+
if (any(abs(resid[i]) > tol for i in range(len(resid)))):
print "Oh no! Bulk composition outside the compositional range of the chosen minerals ( Maximum residual:", max(resid), ")"
@@ -45,27 +74,35 @@ else:
print 'Composition good'
-
-
# Create list of unknowns and their starting guesses
-gvars=['P (kbar)', 'T (K)']
-guesses=[250, 1000, 1.0] # First *mineral* proportion starts at 1.0
-
-gvars.extend(['p(' + mineral.name + ')' for mineral in minlist])
-
-
-
-for mineral in minlist:
+gvars=['P (Pa)', 'T (K)']
+P=1e9
+T=2000.
+guesses=[P, T]
+
+for mineral in mineral_list:
+ gvars.extend(['p(' + mineral.name + ')'])
+ guesses.extend([1.0])
if isinstance(mineral, burnman.SolidSolution):
for i in range(len(mineral.base_material)-1):
gvars.extend(['x(' + str(mineral.base_material[0][0].params['name']) + ')'])
- #guesses.extend(map(float,zip(*gueslist)[1]))
+ if i==0:
+ guesses.extend([1.0])
+ else:
+ guesses.extend([0.0])
print gvars
+print guesses
# Fix T, P, proportions or mixing parameters
fvars=[]
+fvars.append([1, 1000., 2000., 11])
+fvars.append([4, 0., 0., 1.]) # For some reason the solver fails when I set "2, 0., 0., 1", i.e. p(olivine)=0. This seems to be order dependent ... the solver works fine when I change the order of olivine and wadsleyite... :/
+
+'''
+# Fix T, P, proportions or mixing parameters
+fvars=[]
print '\nSelect two variables to fix'
for i in range(len(gvars)):
print i, gvars[i]
@@ -80,30 +117,80 @@ for i in range(2):
else:
print 'Incorrect number of lines'
sys.exit()
+'''
print gvars[fvars[0][0]], fvars[0][1], '<= x <=', fvars[0][2], 'nsteps=', fvars[0][3]
print gvars[fvars[1][0]], fvars[1][1], '<= x <=', fvars[1][2], 'nsteps=', fvars[1][3]
-'''
-# Set up problem and solve it with fsolve
+# Set up problem and solve it
print "\n", ' '.join(x.rjust(10) for x in gvars)
fixed_vars=[]
-for i in linspace(fvars[0][1],fvars[0][2],fvars[0][3]):
+for i in np.linspace(fvars[0][1],fvars[0][2],fvars[0][3]):
guesses[fvars[0][0]]=i
if fvars[0][0] == 2:
guesses[3]=1.0
- for j in linspace(fvars[1][1],fvars[1][2],fvars[1][3]):
+ for j in np.linspace(fvars[1][1],fvars[1][2],fvars[1][3]):
guesses[fvars[1][0]]=j
fixed_vars=[[fvars[0][0], i],[fvars[1][0], j]]
if fvars[1][0] == 2:
guesses[3]=1.0
elif fvars[0][0] == 2 and fvars[1][0] == 3:
guesses[4]=1.0
- soln=fsolve(set_eqns,guesses,args=(comp, minlist, mbrlist, mbrcomp, mbr, model, null, lookup, fixed_vars), full_output=1, xtol=1e-10)
+ soln=opt.fsolve(set_eqns,guesses,args=(comp_vector, mineral_list, stoic, null, fixed_vars), full_output=1, xtol=1e-10)
if soln[2]==1:
print " ".join(str(("%10.5f" % x)) for x in soln[0])
guesses=soln[0]
else:
print " ".join(str(("%10.5f" % x)) for x in soln[0]), soln[3]
-'''
+
+
+# Plot up the last iteration
+
+P=soln[0][0]
+T=soln[0][1]
+
+ol=mg_fe_olivine()
+wad=mg_fe_wadsleyite()
+
+comp = np.linspace(0, 0.2, 100)
+ol_gibbs = np.empty_like(comp)
+wad_gibbs = np.empty_like(comp)
+
+
+for i,c in enumerate(comp):
+ ol.set_composition( np.array([1.0-c, c]) )
+ wad.set_composition( np.array([1.0-c, c]) )
+ ol.set_state( P, T )
+ wad.set_state( P, T )
+ ol_gibbs[i] = ol.gibbs
+ wad_gibbs[i] = wad.gibbs
+
+
+# Detrend
+a=(ol_gibbs[len(comp)-1] - ol_gibbs[0])/(comp[len(comp)-1] - comp[0])
+b=ol_gibbs[0] - a*comp[0]
+
+for i,c in enumerate(comp):
+ ol_gibbs[i] = ol_gibbs[i] - (a*comp[i] + b)
+ wad_gibbs[i] = wad_gibbs[i] - (a*comp[i] + b)
+
+
+comp_tangent=np.array([1.-soln[0][3], 1.-soln[0][5]])
+ol.set_composition( np.array([1.-comp_tangent[0], comp_tangent[0]]) )
+wad.set_composition( np.array([1.-comp_tangent[1], comp_tangent[1]]) )
+ol.set_state( P, T )
+wad.set_state( P, T )
+
+gibbs_tangent=np.array([ol.gibbs - (a*comp_tangent[0] + b), wad.gibbs - (a*comp_tangent[1] + b)])
+
+
+plt.plot( comp, ol_gibbs, 'b', linewidth=1., label='ol')
+plt.plot( comp, wad_gibbs, 'g', linewidth=1., label='wad')
+plt.plot( comp_tangent, gibbs_tangent, 'r', linewidth=1., label='ol-wad equilibrium')
+plt.title('olivine-wadsleyite equilibrium at '+str(round(P/1.e9,2))+' GPa and '+str(round(T,0))+' K')
+plt.xlim(0.0,0.2)
+plt.ylabel("Detrended Gibbs free energy (kJ/mol)")
+plt.xlabel("x(ol), x(wad)")
+plt.legend(loc='upper right')
+plt.show()
More information about the CIG-COMMITS
mailing list