[cig-commits] [commit] add_thermodynamic_potentials: Move things around (d5ddad0)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Dec 9 09:57:28 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_thermodynamic_potentials
Link : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51
>---------------------------------------------------------------
commit d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51
Author: Ian Rose <ian.r.rose at gmail.com>
Date: Tue Oct 7 09:00:30 2014 -0700
Move things around
>---------------------------------------------------------------
d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51
burnman/equilibriumequations.py | 60 ----------------------------
burnman/gibbsminimization.py | 52 ++++++++++++++++++++++++
burnman/test_tc_clone.py => test_tc_clone.py | 12 ++----
3 files changed, 55 insertions(+), 69 deletions(-)
diff --git a/burnman/equilibriumequations.py b/burnman/equilibriumequations.py
deleted file mode 100644
index 1318777..0000000
--- a/burnman/equilibriumequations.py
+++ /dev/null
@@ -1,60 +0,0 @@
-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/gibbsminimization.py b/burnman/gibbsminimization.py
index c77dc09..895054c 100644
--- a/burnman/gibbsminimization.py
+++ b/burnman/gibbsminimization.py
@@ -175,3 +175,55 @@ def sparsify_basis ( basis ):
return new_basis
+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/test_tc_clone.py
similarity index 92%
rename from burnman/test_tc_clone.py
rename to test_tc_clone.py
index e9e905d..5fbd91f 100644
--- a/burnman/test_tc_clone.py
+++ b/test_tc_clone.py
@@ -1,15 +1,9 @@
-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
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.}
mineral_list = [mg_fe_olivine(), mg_fe_wadsleyite()]
@@ -43,8 +37,8 @@ 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 ( mineral_list )
-null=sparsify_basis(compute_nullspace(stoic[0]))
+stoic=burnman.assemble_stoichiometric_matrix ( mineral_list )
+null=burnman.sparsify_basis(burnman.compute_nullspace(stoic[0]))
# Create a compositional vector with elements in the same order as the
# stoichiometric matrix
@@ -136,7 +130,7 @@ for i in np.linspace(fvars[0][1],fvars[0][2],fvars[0][3]):
guesses[3]=1.0
elif fvars[0][0] == 2 and fvars[1][0] == 3:
guesses[4]=1.0
- soln=opt.fsolve(set_eqns,guesses,args=(comp_vector, mineral_list, stoic, null, fixed_vars), full_output=1, xtol=1e-10)
+ soln=opt.fsolve(burnman.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]
More information about the CIG-COMMITS
mailing list