[cig-commits] [commit] inversion, master, validate_MT_params: First alpha of SolidSolution base class (3fd4238)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Fri Dec 12 18:24:54 PST 2014
Repository : https://github.com/geodynamics/burnman
On branches: inversion,master,validate_MT_params
Link : https://github.com/geodynamics/burnman/compare/80c2a295c42dfdb38f83f6c1334bf7d8f97a8463...409647ff05dfad6a686198cac1481bd46b5e2e62
>---------------------------------------------------------------
commit 3fd4238f38679631371933f46990954c0df94c1b
Author: Bob Myhill <myhill.bob at gmail.com>
Date: Sat Aug 30 20:50:05 2014 +0200
First alpha of SolidSolution base class
>---------------------------------------------------------------
3fd4238f38679631371933f46990954c0df94c1b
burnman/solidsolution.py | 156 ++++++++++++++++++++++--------
burnman/test_process_solidsolution.py | 177 ++++++++++++++++++++++++++++++++++
2 files changed, 293 insertions(+), 40 deletions(-)
diff --git a/burnman/solidsolution.py b/burnman/solidsolution.py
index ebb75d2..44198d6 100644
--- a/burnman/solidsolution.py
+++ b/burnman/solidsolution.py
@@ -1,10 +1,12 @@
# BurnMan - a lower mantle toolkit
-# Copyright (C) 2012, 2013, Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
+# Copyright (C) 2012-2014, Myhill, R., Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
# Released under GPL v2 or later.
import warnings
+import re
import numpy as np
+from fractions import Fraction
kd = lambda x,y : 1 if x==y else 0
@@ -29,58 +31,130 @@ class SolidSolution(Mineral):
and P derivatives in J/K/mol and m^3/(mol molecule).
"""
- # init sets up matrices to speed up calculations when P, T, X is defined.
+ # init sets up matrices to speed up calculations for when P, T, X is defined.
def __init__(self, base_material, interaction_parameter):
+ # Initialise the solid soltion inputs
self.base_material = base_material
- self.molar_fraction = molar_fraction
- assert(len(base_material) == len(molar_fraction))
- assert(len(base_material) == len(van_laar_parameter))
+ self.interaction_parameter = interaction_parameter
+
+ # Number of endmembers in the solid solution
+ self.n_endmembers=len(base_material)
+
+ # Check that base_material and interaction parameter each have three parts
+ assert(len(map(list, zip(*base_material))) == 3)
+ assert(len(interaction_parameter) == 3)
+
+ # Check that the interaction parameters have the correct number of variables
+ for i in range(3):
+ assert(len(interaction_parameter[i]) == n_endmembers-1)
+ for j in range(n_endmembers-1):
+ assert(len(interaction_parameter[i][j]) == (n_endmembers-1)-j)
+
+ # Split interaction parameter into H, S, V terms
+ self.excess_enthalpy=interaction_parameter[0]
+ self.excess_entropy=interaction_parameter[1]
+ self.excess_volume=interaction_parameter[2]
+
+ # Number of sites
+ self.n_sites=base_material[0][1].count('[')
+
+ # Check the number of sites is the same for each endmember
+ for i in range(n_endmembers):
+ assert(base_material[i][1].count('[') == n_sites)
+
+ # Check the chemical composition of each endmember is consistent with the dataset
+ # NOT IMPLEMENTED YET
+
+ # Number of unique site occupancies (e.g.. Mg on X etc.)
+ sites=[[] for i in range(n_sites)]
+ list_occupancies=[]
+ list_multiplicity=np.empty(shape=(n_sites))
+ self.n_occupancies=0
+ for endmember in range(n_endmembers):
+ list_occupancies.append([[0]*len(sites[site]) for site in range(n_sites)])
+ s=re.split(r'\[', base_material[endmember][1])[1:]
+ for site in range(n_sites):
+ site_occupancy=re.split(r'\]', s[site])[0]
+ mult=re.split('[A-Z][^A-Z]*',re.split(r'\]', s[site])[1])[0]
+ if mult == '':
+ list_multiplicity[site]=1.0
+ else:
+ list_multiplicity[site]=mult
+ elements=re.findall('[A-Z][^A-Z]*',site_occupancy)
+ for i in range(len(elements)):
+ element_on_site=re.split('[0-9][^A-Z]*',elements[i])[0]
+ proportion_element_on_site=re.findall('[0-9][^A-Z]*',elements[i])
+ if len(proportion_element_on_site) == 0:
+ proportion_element_on_site=Fraction(1.0)
+ else:
+ proportion_element_on_site=Fraction(proportion_element_on_site[0])
+
+ if element_on_site not in sites[site]:
+ self.n_occupancies=n_occupancies+1
+ sites[site].append(element_on_site)
+ element_index=sites[site].index(element_on_site)
+ for parsed_mbr in range(len(list_occupancies)):
+ list_occupancies[parsed_mbr][site].append(0)
+ else:
+ element_index=sites[site].index(element_on_site)
+ list_occupancies[endmember][site][element_index]=proportion_element_on_site
+
+ # Site occupancies and multiplicities
+ self.site_occupancies=np.empty(shape=(n_endmembers,n_occupancies))
+ self.site_multiplicities=np.empty(shape=(n_occupancies))
+ for endmember in range(n_endmembers):
+ n_element=0
+ for site in range(n_sites):
+ for element in range(len(list_occupancies[endmember][site])):
+ self.site_occupancies[endmember][n_element]=list_occupancies[endmember][site][element]
+ self.site_multiplicities[n_element]=list_multiplicity[site]
+ n_element=n_element+1
+
+ self.alpha=np.array([base_material[i][2] for i in range(n_endmembers)])
+
+ self.Wh=np.zeros(shape=(n_endmembers,n_endmembers))
+ self.Ws=np.zeros(shape=(n_endmembers,n_endmembers))
+ self.Wv=np.zeros(shape=(n_endmembers,n_endmembers))
+ for i in range(n_endmembers):
+ for j in range(i+1, n_endmembers):
+ self.Wh[i][j]=2.*excess_enthalpy[i][j-i-1]/(alpha[i]+alpha[j])
+ self.Ws[i][j]=2.*excess_entropy[i][j-i-1]/(alpha[i]+alpha[j])
+ self.Wv[i][j]=2.*excess_volume[i][j-i-1]/(alpha[i]+alpha[j])
+
+ # END INITIALISATION
+
+
+ def set_composition(self, molar_fraction):
+ assert(len(self.base_material) == len(molar_fraction))
assert(sum(molar_fraction) > 0.9999)
assert(sum(molar_fraction) < 1.0001)
- # make sure the number of atoms in each endmember is the same
- # (when we start adding minerals with vacancies
- # this will need to be changed).
- for m in base_material:
- if(base_material[0].params.has_key('n')):
- assert(m.params['n'] == base_material[0].params['n'])
+ self.molar_fraction=molar_fraction
- # Set the state of all the endmembers
- def set_state(self, pressure, temperature):
- for mat in self.base_materials:
- mat.method = self.method
- mat.set_state(pressure, temperature)
+ phi=np.array([self.alpha[i]*molar_fraction[i] for i in range(self.n_endmembers)])
+ phi=np.divide(phi, np.sum(phi))
- itrange = range(0, len(self.base_material))
+ self.H_excess=np.dot(self.alpha.T,molar_fraction)*np.dot(phi.T,np.dot(self.Wh,phi))
+ self.S_deficit=np.dot(self.alpha.T,molar_fraction)*np.dot(phi.T,np.dot(self.Ws,phi))
+ self.V_excess=np.dot(self.alpha.T,molar_fraction)*np.dot(phi.T,np.dot(self.Wv,phi))
+ self.S_excess=0.0-S_deficit
- self.params = {}
-
- # NEEDS CHANGING FROM HERE!!!
- # Description of RTlngamma terms given by the asymmetric formalism
- # is included in Holland and Powell (2003)
-
- # the volume is calculated in the same way as Gex (p.493),
- # replacing Wij with Wvij
- Vex=-2.*sum([ sum([ self.molar_fraction[i]*self.van_laar_parameter[i]*self.molar_fraction[j]*self.van_laar_parameter[j]*self.interaction_parameter[i][2] / (self.van_laar_parameter[i] + self.van_laar_parameter[j]) for j in range(i, len(self.base_material)) ]) for i in range(0, len(self.base_material)-1) ])/sum([ self.van_laar_parameter[l]*self.molar_fraction[l] for l in itrange ])
- V= sum([ self.base_materials[i].molar_volume() * self.molar_fraction[i] for i in itrange ]) + Vex
-
- # the volume is calculated in the same way as Gex (p.493),
- # replacing Wij with Wvij
- # NEED TO CORRECT INDICES FOR INT PARAM
- # W[i][j] == W[(itrange-1)*i+(j-1)]
- Vex=-2.*sum([ sum([ self.molar_fraction[i]*self.van_laar_parameter[i]*self.molar_fraction[j]*self.van_laar_parameter[j]*self.interaction_parameter[(itrange-1)*i+(j-1)][2] / (self.van_laar_parameter[i] + self.van_laar_parameter[j]) for j in range(i, len(self.base_material)) ]) for i in range(0, len(self.base_material)-1) ])/sum([ self.van_laar_parameter[l]*self.molar_fraction[l] for l in itrange ])
- V= sum([ self.base_materials[i].molar_volume() * self.molar_fraction[i] for i in itrange ]) + Vex
-
-
- # Same for the Gibbs free energy
- # NEED TO CORRECT INDICES FOR INT PARAM
- Gex=-2.*sum([ sum([ self.molar_fraction[i]*self.van_laar_parameter[i]*self.molar_fraction[j]*self.van_laar_parameter[j]*(pressure*self.interaction_parameter[(itrange-1)*i+(j-1)][0] + temperature*self.interaction_parameter[(itrange-1)*i+(j-1)][1] + pressure*self.interaction_parameter[(itrange-1)*i+(j-1)][2]) / (self.van_laar_parameter[i] + self.van_laar_parameter[j]) for j in range(i, len(self.base_material)) ]) for i in range(0, len(self.base_material)-1) ])/sum([ self.van_laar_parameter[l]*self.molar_fraction[l] for l in itrange ])
- G= sum([ self.base_materials[i].molar_gibbs() * self.molar_fraction[i] for i in itrange ]) + Gex
+ def set_state(self, pressure, temperature, molar_fraction):
+ # Set the state of all the endmembers
+ for i in range(self.n_endmembers):
+ self.base_material[i][0].method = self.method
+ self.base_material[i][0].set_state(pressure, temperature)
+ # Find excess properties
+ self.set_composition(molar_fraction)
+ self.G_excess=self.H_excess - temperature*self.S_excess + pressure*self.V_excess
+ self.molar_volume= sum([ self.base_material[i][0].molar_volume() * self.molar_fraction[i] for i in n_endmembers ]) + self.V_excess
+ self.molar_gibbs= sum([ self.base_materials[i][0].molar_gibbs() * self.molar_fraction[i] for i in n_endmembers ]) + self.G_excess
+ '''
for prop in self.base_materials[0].params:
try:
self.params[prop] = sum([ self.base_materials[i].params[prop] * self.molar_fraction[i] for i in itrange ])
@@ -88,4 +162,6 @@ class SolidSolution(Mineral):
#if there is a type error, it is probably a string. Just go with the value of the first base_material.
self.params[prop] = self.base_materials[0].params[prop]
Mineral.set_state(self, pressure, temperature)
+ '''
+
diff --git a/burnman/test_process_solidsolution.py b/burnman/test_process_solidsolution.py
new file mode 100644
index 0000000..92fbbde
--- /dev/null
+++ b/burnman/test_process_solidsolution.py
@@ -0,0 +1,177 @@
+# BurnMan - a lower mantle toolkit
+# Copyright (C) 2012-2014, Myhill, R., Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
+# Released under GPL v2 or later.
+
+# This is a test script which is the precursor to implementing a solid solution BaseClass.
+
+
+'''
+# Inputs
+Solid solution model
+Endmember proportions
+P,T
+
+# Things we can initialise without endmember proportions
+n_sites
+n_occupancies
+n_endmembers
+alpha
+Wh, Ws, Wv
+site_occupancies
+site_multiplicities
+
+# Endmember proportions needed
+phi
+ideal activities
+occupancies
+
+# P, T needed
+Wtotal / nonideal contribution to gibbs
+'''
+
+import re
+import numpy as np
+from fractions import Fraction
+
+# Endmembers
+base_material = [['pyrope()', '[Mg]3[Al]2Si3O12', 1.0], \
+ ['almandine()', '[Fe]3[Al]2Si3O12', 1.0], \
+ ['grossular()', '[Ca]3[Al]2Si3O12', 2.7], \
+ ['majorite()', '[Mg]3[Mg1/2Si1/2]2Si3O12', 1.0]]
+
+n_endmembers=len(base_material)
+
+# Interaction parameters
+excess_enthalpy=[[2.5e3, 30.1e3, 15e3],[10e3,18e3],[48e3]]
+excess_entropy=[[0., 0., 0.],[0., 0.],[0.]]
+excess_volume=[[0., 0.164e-5, 0.],[0., 0.],[0.]]
+interaction_parameter=[excess_enthalpy,excess_entropy,excess_volume]
+
+
+# INPUT PROPORTIONS
+molar_fraction = np.array([ 0.5, 0.2, 0.1, 0.2 ])
+
+
+# "sites" is a 2D list of sites and the elements which reside on them
+# "list_occupancies" is a 3D list describing the elemental site occupancies of each endmember
+# "site_occupancies" is a 2D np.array of list_occupancies, concatenating the 2nd and 3rd dimension
+n_sites=base_material[0][1].count('[')
+print 'Number of sites:', n_sites
+print ''
+
+sites=[[] for i in range(n_sites)]
+list_occupancies=[]
+list_multiplicity=np.empty(shape=(n_sites))
+n_occupancies=0
+for endmember in range(n_endmembers):
+ list_occupancies.append([[0]*len(sites[site]) for site in range(n_sites)])
+ s=re.split(r'\[', base_material[endmember][1])[1:]
+ for site in range(n_sites):
+ site_occupancy=re.split(r'\]', s[site])[0]
+ mult=re.split('[A-Z][^A-Z]*',re.split(r'\]', s[site])[1])[0]
+ if mult == '':
+ list_multiplicity[site]=1.0
+ else:
+ list_multiplicity[site]=mult
+ elements=re.findall('[A-Z][^A-Z]*',site_occupancy)
+ for i in range(len(elements)):
+ element_on_site=re.split('[0-9][^A-Z]*',elements[i])[0]
+ proportion_element_on_site=re.findall('[0-9][^A-Z]*',elements[i])
+ if len(proportion_element_on_site) == 0:
+ proportion_element_on_site=Fraction(1.0)
+ else:
+ proportion_element_on_site=Fraction(proportion_element_on_site[0])
+
+ if element_on_site not in sites[site]:
+ n_occupancies=n_occupancies+1
+ sites[site].append(element_on_site)
+ element_index=sites[site].index(element_on_site)
+ for parsed_mbr in range(len(list_occupancies)):
+ list_occupancies[parsed_mbr][site].append(0)
+ else:
+ element_index=sites[site].index(element_on_site)
+ list_occupancies[endmember][site][element_index]=proportion_element_on_site
+
+
+
+site_occupancies=np.empty(shape=(n_endmembers,n_occupancies))
+site_multiplicities=np.empty(shape=(n_occupancies))
+for endmember in range(n_endmembers):
+ n_element=0
+ for site in range(n_sites):
+ for element in range(len(list_occupancies[endmember][site])):
+ site_occupancies[endmember][n_element]=list_occupancies[endmember][site][element]
+ site_multiplicities[n_element]=list_multiplicity[site]
+ n_element=n_element+1
+
+
+# Matrix operation to return site occupancies, given proportions of endmembers
+occupancies=np.dot(molar_fraction, site_occupancies)
+print 'Site occupancies'
+print sites
+print occupancies
+print ''
+
+
+
+# Ideal activities
+ideal_activity=np.empty(shape=(n_endmembers))
+for endmember in range(n_endmembers):
+ ideal_activity[endmember]=1.0
+ normalisation_constant=1.0
+ for element in range(n_occupancies):
+ if site_occupancies[endmember][element] != 0:
+ #print base_material[endmember][0], element, occupancies[element], site_occupancies[endmember][element]
+ ideal_activity[endmember]=ideal_activity[endmember]*pow(occupancies[element],site_multiplicities[element])
+ normalisation_constant=normalisation_constant/pow(site_occupancies[endmember][element],site_multiplicities[element])
+ ideal_activity[endmember]=normalisation_constant*ideal_activity[endmember]
+
+
+print 'Ideal activities'
+print ideal_activity
+
+# Non ideal activities
+# alpha^T*p*(phi^T*W*phi)
+
+alpha=np.array([base_material[i][2] for i in range(n_endmembers)])
+phi=np.array([alpha[i]*molar_fraction[i] for i in range(n_endmembers)])
+phi=np.divide(phi, np.sum(phi))
+
+Wh=np.zeros(shape=(n_endmembers,n_endmembers))
+Ws=np.zeros(shape=(n_endmembers,n_endmembers))
+Wv=np.zeros(shape=(n_endmembers,n_endmembers))
+for i in range(n_endmembers):
+ for j in range(i+1, n_endmembers):
+ Wh[i][j]=2.*excess_enthalpy[i][j-i-1]/(alpha[i]+alpha[j])
+ Ws[i][j]=2.*excess_entropy[i][j-i-1]/(alpha[i]+alpha[j])
+ Wv[i][j]=2.*excess_volume[i][j-i-1]/(alpha[i]+alpha[j])
+
+print ''
+print 'Excess enthalpy, entropy, volume (non-configurational)'
+print np.dot(alpha.T,molar_fraction)*np.dot(phi.T,np.dot(Wh,phi)), 'J/mol'
+
+print np.dot(alpha.T,molar_fraction)*np.dot(phi.T,np.dot(Ws,phi)), 'J/K/mol'
+
+print np.dot(alpha.T,molar_fraction)*np.dot(phi.T,np.dot(Wv,phi)), 'm^3/mol'
+
+
+
+# Plot excess volumes for the pyrope-grossular join
+ppy= np.linspace(0, 1, 101)
+vex= np.empty(shape=(101))
+for p in range(len(ppy)):
+ a=ppy[p]
+ molar_fraction = np.array([ a, 0.0, 1.-a, 0.0 ])
+ phi=np.array([alpha[i]*molar_fraction[i] for i in range(n_endmembers)])
+ phi=np.divide(phi, np.sum(phi))
+
+ vex[p]=np.dot(alpha.T,molar_fraction)*np.dot(phi.T,np.dot(Wv,phi))
+
+
+import matplotlib.pyplot as plt
+plt.plot(ppy,vex,color='r',linestyle='-',marker='o',markerfacecolor='r',markersize=0)
+plt.xlim(min(ppy),max(ppy))
+plt.xlabel("p(pyrope)")
+plt.title("V excess (m^3/mol) for pyrope-grossular garnets")
+
+plt.show()
More information about the CIG-COMMITS
mailing list