[cig-commits] [commit] add_thermodynamic_potentials: Disorganized work on getting dGdN stuff working (86bc2e5)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Dec 10 14:49:05 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_thermodynamic_potentials
Link : https://github.com/geodynamics/burnman/compare/8b9fb89bc6a4ecfef8ceaffb61e81e8da56c09ef...86bc2e51c012edd73c6b8365a35a8445d216a9f4
>---------------------------------------------------------------
commit 86bc2e51c012edd73c6b8365a35a8445d216a9f4
Author: Ian Rose <ian.r.rose at gmail.com>
Date: Wed Oct 22 11:44:25 2014 -0700
Disorganized work on getting dGdN stuff working
Conflicts:
burnman/solidsolution.py
>---------------------------------------------------------------
86bc2e51c012edd73c6b8365a35a8445d216a9f4
burnman/equilibriumassemblage.py | 83 ++++++++++++++++++++++++++-----------
burnman/solidsolution.py | 10 ++---
burnman/solutionmodel.py | 89 ++++++++++++++++++++++------------------
olivine_wadsleyite.py | 10 +++--
solidsolution_benchmarks.py | 4 +-
5 files changed, 122 insertions(+), 74 deletions(-)
diff --git a/burnman/equilibriumassemblage.py b/burnman/equilibriumassemblage.py
index 135e9b8..84aaba6 100644
--- a/burnman/equilibriumassemblage.py
+++ b/burnman/equilibriumassemblage.py
@@ -83,21 +83,52 @@ class EquilibriumAssemblage(burnman.Material):
temperature : float
The temperature in [K] at which to set the state.
"""
+ self._minimize_gibbs( pressure, temperature )
+# self._solve_equilibrium_equations (pressure, temperature )
+
+
+ def _solve_equilibrium_equations( self, pressure, temperature ):
+
+ ref_gibbs = self.__compute_gibbs( pressure, temperature, self.__compute_species_vector(self.reduced_species_vector * 0.))
+ print ref_gibbs
+ equilibrium = lambda x : (self._equilibrium_equations( pressure, temperature, self.__compute_species_vector( x ) ))
+ sol = opt.root(equilibrium, self.reduced_species_vector, method='hybr', options={'xtol':1.e-11})
+
+ self.reduced_species_vector = sol.x
+ self.species_vector = self.__compute_species_vector(self.reduced_species_vector)
+ self.gibbs = self.__compute_gibbs(pressure, temperature, self.species_vector)
+
+ def _equilibrium_equations ( self, pressure, temperature, species_vector):
+ partial_gibbs = self._compute_partial_gibbs( pressure, temperature, species_vector)
+ gibbs_inequalities = np.dot(partial_gibbs, self.right_nullspace)
+
+ if np.any(species_vector < -1.e-6):
+ gibbs_inequalitites = np.ones_like( gibbs_inequalities) * 10000000.
+
+ return gibbs_inequalities
+
+ def _minimize_gibbs( self, pressure, temperature):
n = len(self.endmember_formulae)
# Calculate a reference gibbs free energy at the baseline assemblage. This is a kind of ill-conditioned
# minimization problem, and it seems to work better if we subtract off a reference state and normalize.
- ref_gibbs = self.__compute_gibbs( pressure, temperature, self.reduced_species_vector * 0.)
+ ref_gibbs = self.__compute_gibbs( pressure, temperature, self.__compute_species_vector(self.reduced_species_vector * 0.))
#define the function to minimize and then do it.
- 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.__compute_species_vector(x), full_output=1, disp=False)
+ minimize_gibbs = lambda x : (self.__compute_gibbs( pressure, temperature, self.__compute_species_vector(x) ) - ref_gibbs)/np.abs(ref_gibbs)
+# minimize_gibbs = lambda x : (self.__compute_gibbs( pressure, temperature, self.__compute_species_vector(x) ))
+ equilibrium = lambda x : (self._equilibrium_equations( pressure, temperature, self.__compute_species_vector( x ) )/np.abs(ref_gibbs))
+# equilibrium=None
#Set the solution
- self.reduced_species_vector = sol[0]
+ constraints={'type':'ineq', 'fun':self.__compute_species_vector}
+ sol = opt.minimize( minimize_gibbs, self.reduced_species_vector, method='SLSQP', jac=equilibrium, constraints=constraints )
+ self.reduced_species_vector = sol.x
self.species_vector = self.__compute_species_vector(self.reduced_species_vector)
- self.gibbs = sol[1]*np.abs(ref_gibbs) + ref_gibbs
+ self.gibbs = self.__compute_gibbs(pressure, temperature, self.species_vector)
+
+
def print_assemblage(self):
@@ -108,25 +139,13 @@ class EquilibriumAssemblage(burnman.Material):
for f,s in zip(self.endmember_formulae, self.species_vector):
print f, s/tot
-
- def __compute_gibbs( self, pressure, temperature, reduced_vector ):
- """
- Given a pressure, temperature, and vector in the nullspace,
- calculate the gibbs free energy of the assemblage. This
- is basically the function to minimize when taking the
- assemblage to a P-T, subject to the bulk composition contraint
- (which is parameterized by vectors in the nullspace)
- """
-
- species_vector = self.__compute_species_vector(reduced_vector)
- assert( len(species_vector) == len(self.endmember_formulae) )
-
- tmp_gibbs = 0.0
- i = 0
+ def _compute_partial_gibbs( self, pressure, temperature, species_vector):
+ partial_gibbs = np.empty_like(species_vector)
#Loop over the various phases and compute the gibbs free energy
#of each one at P,T. We have to treat solid solutions and
#single phases somewhat differently.
+ i = 0
for phase in self.phases:
if isinstance (phase, burnman.SolidSolution):
@@ -136,17 +155,32 @@ class EquilibriumAssemblage(burnman.Material):
phase.set_composition( molar_fractions )
phase.set_state( pressure, temperature )
- tmp_gibbs += phase.gibbs * (frac if frac > 1.e-6 else 0.0)
+ partial_gibbs[i:(i+n)] = phase.partial_gibbs
i+=n
elif isinstance(phase, burnman.Mineral):
phase.set_state( pressure, temperature )
- tmp_gibbs += phase.gibbs * species_vector[i]
+ partial_gibbs[i] = phase.gibbs
i+=1
else:
raise Exception('Unsupported mineral type, can only read burnman.Mineral or burnman.SolidSolution')
- return tmp_gibbs
+ return partial_gibbs
+
+
+ def __compute_gibbs( self, pressure, temperature, species_vector ):
+ """
+ Given a pressure, temperature, and vector in the nullspace,
+ calculate the gibbs free energy of the assemblage. This
+ is basically the function to minimize when taking the
+ assemblage to a P-T, subject to the bulk composition contraint
+ (which is parameterized by vectors in the nullspace)
+ """
+
+ assert( len(species_vector) == len(self.endmember_formulae) )
+
+ gibbs = np.dot( species_vector, self._compute_partial_gibbs(pressure, temperature, species_vector) )
+ return gibbs
def __compute_species_vector ( self, reduced_vector ):
@@ -154,7 +188,8 @@ class EquilibriumAssemblage(burnman.Material):
Given a vector in the nullspace, return a full species vector
in the endmember space.
"""
- return self.baseline_assemblage + np.dot( self.right_nullspace, np.transpose(reduced_vector) )
+ species_vector = self.baseline_assemblage + np.dot( self.right_nullspace, np.transpose(reduced_vector) )
+ return species_vector
def __compute_bulk_composition (self, species_vector):
diff --git a/burnman/solidsolution.py b/burnman/solidsolution.py
index e03676c..3c488bc 100644
--- a/burnman/solidsolution.py
+++ b/burnman/solidsolution.py
@@ -52,11 +52,14 @@ class SolidSolution(Mineral):
for i in range(self.n_endmembers):
self.base_material[i][0].set_state(pressure, temperature)
+ self.excess_partial_gibbs = self.solution_model.excess_partial_gibbs_free_energies( pressure, temperature, self.molar_fraction)
self.excess_gibbs = self.solution_model.excess_gibbs_free_energy( pressure, temperature, self.molar_fraction)
+ self.partial_gibbs = np.array([self.base_material[i][0].gibbs for i in range(self.n_endmembers)]) + self.excess_partial_gibbs
+ self.gibbs= sum([ self.base_material[i][0].gibbs * self.molar_fraction[i] for i in range(self.n_endmembers) ]) + self.excess_gibbs
+
self.excess_volume = self.solution_model.excess_volume( pressure, temperature, self.molar_fraction)
self.V= sum([ self.base_material[i][0].V * self.molar_fraction[i] for i in range(self.n_endmembers) ]) + self.excess_volume
- self.gibbs= sum([ self.base_material[i][0].gibbs * self.molar_fraction[i] for i in range(self.n_endmembers) ]) + self.excess_gibbs
'''
@@ -74,7 +77,4 @@ class SolidSolution(Mineral):
return sum([ self.base_material[i][0].calcgibbs(pressure, temperature) * molar_fractions[i] for i in range(self.n_endmembers) ]) + self.solution_model.excess_gibbs_free_energy( pressure, temperature, molar_fractions)
def calcpartialgibbsexcesses(self, pressure, temperature, molar_fractions):
- Hint, Sint, Vint = self.solution_model.non_ideal_interactions(molar_fractions)
- partialgibbsexcesses=np.empty(len(molar_fractions))
- partialgibbsexcesses = np.array([0.+R*temperature*self.solution_model.ln_ideal_activities(molar_fractions)[i] + Hint[i] - temperature*Sint[i] + pressure*Vint[i] for i in range(self.n_endmembers)])
- return partialgibbsexcesses
+ return self.solution_model.excess_partial_gibbs_free_energies( pressure, temperature, molar_fractions)
diff --git a/burnman/solutionmodel.py b/burnman/solutionmodel.py
index ad97188..a812b49 100644
--- a/burnman/solutionmodel.py
+++ b/burnman/solutionmodel.py
@@ -25,8 +25,9 @@ class SolutionModel:
def excess_gibbs_free_energy( self, pressure, temperature, molar_fractions):
"""
Given a list of molar fractions of different phases,
- compute the excess Gibbs free energy of the solution,
- relative to an ideal model.
+ compute the excess Gibbs free energy of the solution.
+ The base class implementation assumes that the excess gibbs
+ free energy is zero.
Parameters
----------
@@ -44,14 +45,38 @@ class SolutionModel:
G_excess : float
The excess Gibbs free energy
"""
- return 0.0
+ return np.dot(np.array(molar_fractions), self.excess_partial_gibbs_free_energies( pressure, temperature, molar_fractions))
+
+ def excess_partial_gibbs_free_energies( self, pressure, temperature, molar_fractions):
+ """
+ Given a list of molar fractions of different phases,
+ compute the excess Gibbs free energy for each endmember of the solution.
+ The base class implementation assumes that the excess gibbs
+ free energy is zero.
+
+ Parameters
+ ----------
+ pressure : float
+ Pressure at which to evaluate the solution model. [Pa]
+
+ temperature : float
+ Temperature at which to evaluate the solution. [K]
+
+ molar_fractions : list of floats
+ List of molar fractions of the different endmembers in solution
+
+ Returns
+ -------
+ partial_G_excess : numpy array
+ The excess Gibbs free energy of each endmember
+ """
+ return np.empty_like( np.array(molar_fractions) )
def excess_volume( self, pressure, temperature, molar_fractions):
"""
Given a list of molar fractions of different phases,
compute the excess Gibbs free energy of the solution,
- relative to an ideal model. The base class implementation
- assumes that the excess volume is zero.
+ The base class implementation assumes that the excess volume is zero.
Parameters
----------
@@ -87,6 +112,12 @@ class IdealSolution ( SolutionModel ):
self.solution_formulae, self.n_sites, self.sites, self.n_occupancies, self.endmember_occupancies, self.site_multiplicities = \
ProcessSolidSolutionChemistry(self.formulas)
+ self._calculate_endmember_configurational_entropies()
+
+ def excess_partial_gibbs_free_energies( self, pressure, temperature, molar_fractions ):
+ return self._ideal_excess_partial_gibbs( temperature, molar_fractions )
+
+ def _calculate_endmember_configurational_entropies( self ):
self.endmember_configurational_entropies=np.zeros(shape=(self.n_endmembers))
for idx, endmember_occupancy in enumerate(self.endmember_occupancies):
for occ in range(self.n_occupancies):
@@ -95,34 +126,23 @@ class IdealSolution ( SolutionModel ):
self.endmember_configurational_entropies[idx] - \
R*self.site_multiplicities[occ]*endmember_occupancy[occ]*np.log(endmember_occupancy[occ])
+ def _endmember_configurational_entropy_contribution(self, molar_fractions):
+ return np.dot(molar_fractions, self.endmember_configurational_entropies)
-
- def excess_gibbs_free_energy( self, pressure, temperature, molar_fractions ):
- return self.ideal_gibbs_excess( temperature, molar_fractions )
-
- def configurational_entropy (self, molar_fractions):
+ def _configurational_entropy (self, molar_fractions):
site_occupancies=np.dot(molar_fractions, self.endmember_occupancies)
conf_entropy=0
for idx, occupancy in enumerate(site_occupancies):
if occupancy > 1e-10:
conf_entropy=conf_entropy-R*occupancy*self.site_multiplicities[idx]*np.log(occupancy)
- # Alternative entropy construction
- #activities=self.ideal_activities(molar_fractions)
- #conf_entropy=0.
- #for i, activity in enumerate(activities):
- # if activity > 1e-10:
- # conf_entropy=conf_entropy - R*molar_fractions[i]*np.log(activity) + molar_fractions[i]*self.endmember_configurational_entropies[i]
-
return conf_entropy
- def endmember_configurational_entropy_contribution(self, molar_fractions):
- return np.dot(molar_fractions, self.endmember_configurational_entropies)
- def ideal_gibbs_excess( self, temperature, molar_fractions ):
- return 0.0-temperature*(self.configurational_entropy(molar_fractions) - self.endmember_configurational_entropy_contribution(molar_fractions))
+ def _ideal_excess_partial_gibbs( self, temperature, molar_fractions ):
+ return R * temperature * self._log_ideal_activities(molar_fractions)
- def ln_ideal_activities ( self, molar_fractions ):
+ def _log_ideal_activities ( self, molar_fractions ):
site_occupancies=np.dot(molar_fractions, self.endmember_occupancies)
lna=np.empty(shape=(self.n_endmembers))
@@ -137,7 +157,7 @@ class IdealSolution ( SolutionModel ):
return lna
- def ideal_activities ( self, molar_fractions ):
+ def _ideal_activities ( self, molar_fractions ):
site_occupancies=np.dot(molar_fractions, self.endmember_occupancies)
activities=np.empty(shape=(self.n_endmembers))
@@ -186,10 +206,8 @@ class AsymmetricRegularSolution ( IdealSolution ):
#initialize ideal solution model
IdealSolution.__init__(self, endmembers )
- def configurational_entropy( self, molar_fractions ):
- return IdealSolution.configurational_entropy( self, molar_fractions )
- def non_ideal_interactions( self, molar_fractions ):
+ def _non_ideal_interactions( self, molar_fractions ):
# -sum(sum(qi.qj.Wij*)
# equation (2) of Holland and Powell 2003
phi=np.array([self.alpha[i]*molar_fractions[i] for i in range(self.n_endmembers)])
@@ -210,24 +228,17 @@ class AsymmetricRegularSolution ( IdealSolution ):
return Hint, Sint, Vint
+ def _non_ideal_excess_partial_gibbs( self, pressure, temperature, molar_fractions) :
- def excess_gibbs_free_energy( self, pressure, temperature, molar_fractions ):
+ Hint, Sint, Vint = self._non_ideal_interactions( molar_fractions )
+ return Hint - temperature*Sint + pressure*Vint
- ideal_gibbs = IdealSolution.ideal_gibbs_excess( self, temperature, molar_fractions )
-
-
- phi=np.array([self.alpha[i]*molar_fractions[i] for i in range(self.n_endmembers)])
- phi=np.divide(phi, np.sum(phi))
-
- H_excess=np.dot(self.alpha.T,molar_fractions)*np.dot(phi.T,np.dot(self.Wh,phi))
- S_excess=np.dot(self.alpha.T,molar_fractions)*np.dot(phi.T,np.dot(self.Ws,phi))
- V_excess=np.dot(self.alpha.T,molar_fractions)*np.dot(phi.T,np.dot(self.Wv,phi))
-
- non_ideal_gibbs = H_excess - temperature*S_excess + pressure*V_excess
+ def excess_partial_gibbs_free_energies( self, pressure, temperature, molar_fractions ):
+ ideal_gibbs = IdealSolution._ideal_excess_partial_gibbs (self, temperature, molar_fractions )
+ non_ideal_gibbs = self._non_ideal_excess_partial_gibbs( pressure, temperature, molar_fractions)
return ideal_gibbs + non_ideal_gibbs
-
def excess_volume ( self, pressure, temperature, molar_fractions ):
phi=np.array([self.alpha[i]*molar_fractions[i] for i in range(self.n_endmembers)])
diff --git a/olivine_wadsleyite.py b/olivine_wadsleyite.py
index 4415ec1..345fd35 100644
--- a/olivine_wadsleyite.py
+++ b/olivine_wadsleyite.py
@@ -4,7 +4,8 @@ import numpy as np
import matplotlib.pyplot as plt
-composition = { 'Fe': 1./35., 'Mg': 9./35., 'O': 4./7., 'Si': 1./7.}
+#composition = { 'Fe': 1./35., 'Mg': 9./35., 'O': 4./7., 'Si': 1./7.}
+composition = { 'Mg': 1.8, 'Fe': 0.2, 'O': 4., 'Si': 1.}
minlist = [mg_fe_olivine(), mg_fe_wadsleyite()]
ol = mg_fe_olivine()
@@ -17,7 +18,7 @@ ea = burnman.EquilibriumAssemblage(composition, minlist)
ea.set_method('slb3')
co = burnman.Composite( [0.5, 0.5], [ol, wa] )
-pressures = np.linspace( 1.e9, 30.e9, 100 )
+pressures = np.linspace( 12.e9, 14.e9, 50 )
temperature = 1500.
Gol = np.empty_like(pressures)
Gwa = np.empty_like(pressures)
@@ -29,14 +30,15 @@ for i,p in enumerate(pressures):
ol.set_state(p,temperature)
wa.set_state(p,temperature)
ea.set_state(p,temperature)
+ ea.print_assemblage()
Gol[i] = ol.gibbs
Gwa[i] = wa.gibbs
- Gea[i] = ea.gibbs*7.
+ Gea[i] = ea.gibbs
Gav = (Gol+Gwa)*1./2.
plt.plot(pressures, Gol-Gav)
plt.plot(pressures, Gwa-Gav)
-plt.plot(pressures, Gea-Gav, 'k--', linewidth=4)
+plt.plot(pressures, Gea-Gav, 'k-', linewidth=2)
plt.show()
diff --git a/solidsolution_benchmarks.py b/solidsolution_benchmarks.py
index 92fef5d..b1a16e1 100644
--- a/solidsolution_benchmarks.py
+++ b/solidsolution_benchmarks.py
@@ -55,7 +55,7 @@ for i,c in enumerate(comp):
molar_fractions=[1.0-c, c]
sp.set_composition( np.array(molar_fractions) )
sp.set_state( 1e5, 298.15 )
- sp_entropies[i] = sp.solution_model.configurational_entropy( molar_fractions )
+ sp_entropies[i] = sp.solution_model._configurational_entropy( molar_fractions )
sp_entropies_NK1967[i] = -8.3145*(c*np.log(c) + (1.-c)*np.log(1.-c) + c*np.log(c/2.) + (2.-c)*np.log(1.-c/2.)) # eq. 7 in Navrotsky and Kleppa, 1967.
#fig1 = mpimg.imread('configurational_entropy.png') # Uncomment these two lines if you want to overlay the plot on a screengrab from SLB2011
@@ -177,7 +177,7 @@ for idx, model in enumerate(opx_models):
molar_fractions=[1.0-c, c]
model.set_composition( np.array(molar_fractions) )
model.set_state( 0.0, 0.0 )
- opx_entropies[idx][i] = model.solution_model.configurational_entropy( molar_fractions )
+ opx_entropies[idx][i] = model.solution_model._configurational_entropy( molar_fractions )
More information about the CIG-COMMITS
mailing list