[cig-commits] [commit] add_thermodynamic_potentials: SolidSolution test now includes site occupancies and ideal activities (c6b4ada)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Dec 9 09:54:17 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_thermodynamic_potentials
Link : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51
>---------------------------------------------------------------
commit c6b4ada97fe2e2ef1a57276b07c82bed7d70067a
Author: Bob Myhill <myhill.bob at gmail.com>
Date: Thu Aug 28 01:11:43 2014 +0200
SolidSolution test now includes site occupancies and ideal activities
>---------------------------------------------------------------
c6b4ada97fe2e2ef1a57276b07c82bed7d70067a
burnman/test_process_solidsolution.py | 114 +++++++++++++++++++++++++---------
1 file changed, 83 insertions(+), 31 deletions(-)
diff --git a/burnman/test_process_solidsolution.py b/burnman/test_process_solidsolution.py
index bb0e265..978136b 100644
--- a/burnman/test_process_solidsolution.py
+++ b/burnman/test_process_solidsolution.py
@@ -2,56 +2,108 @@
# Copyright (C) 2012, 2013, Myhill, R., Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
# Released under GPL v2 or later.
-# This is a function which returns the number of atoms and molar mass of a compound given its unit formula as an argument.
+# This is a test script which is the precursor to implementing a solid solution BaseClass.
import re
+import numpy as np
+from fractions import Fraction
- # Endmembers
+# 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]]
+# Interaction parameters
+excess_enthalpy=[[2.5e3, 29.1e3, 15e3][10e3,18e3][48e3]]
+excess_entropy=[[0., 0., 0.][0., 0.][0.]]
+excess_volume=[[0., 0., 0.][0., 0.][0.]]
+interaction_parameter=[excess_enthalpy,excess_entropy,excess_volume]
+
+
+
+endmember_proportions = np.array([ 0.5, 0.2, 0.1, 0.2 ])
+#endmember_proportions = np.array([ 0.0, 0.0, 0.0, 1.0 ])
+#endmember_proportions = np.array([ 1.0, 0.0, 0.0, 0.0 ])
+
+# "sites" is a 2D list of sites and the elements which reside on them
+# "site_occupancies" is a 3D list describing the elemental site occupancies of each endmember
+# "array_occupancies" is a 2D np.array of site_occupancies, concatenating the 2nd and 3rd dimension
nsites=base_material[0][1].count('[')
print 'Number of sites:', nsites
+print ''
-for k in range(len(base_material)):
- print base_material[k][0]
- sites=re.split(r'\[', base_material[k][1])[1:]
- for i in range(nsites):
- site_occupancy=re.split(r'\]', sites[i])[0]
- site_multiplicity=re.split('[A-Z][^A-Z]*',re.split(r'\]', sites[i])[1])[0]
- print 'Multiplicity on site', i+1, ':', site_multiplicity
+sites=[[] for i in range(nsites)]
+site_occupancies=[]
+site_multiplicity=np.empty(shape=(nsites))
+n_elements=0
+for endmember in range(len(base_material)):
+ site_occupancies.append([[0]*len(sites[site]) for site in range(nsites)])
+ s=re.split(r'\[', base_material[endmember][1])[1:]
+ for site in range(nsites):
+ site_occupancy=re.split(r'\]', s[site])[0]
+ mult=re.split('[A-Z][^A-Z]*',re.split(r'\]', s[site])[1])[0]
+ if mult == '':
+ site_multiplicity[site]=1.0
+ else:
+ site_multiplicity[site]=mult
elements=re.findall('[A-Z][^A-Z]*',site_occupancy)
- for j in range(len(elements)):
- element_on_site=re.split('[0-9][^A-Z]*',elements[j])[0]
- proportion_element_on_site=re.findall('[0-9][^A-Z]*',elements[j])
+ 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=1.0
+ proportion_element_on_site=Fraction(1.0)
else:
- proportion_element_on_site=proportion_element_on_site[0]
- print element_on_site, proportion_element_on_site
- print ''
+ proportion_element_on_site=Fraction(proportion_element_on_site[0])
-print ''
+ if element_on_site not in sites[site]:
+ n_elements=n_elements+1
+ sites[site].append(element_on_site)
+ element_index=sites[site].index(element_on_site)
+ for parsed_mbr in range(len(site_occupancies)):
+ site_occupancies[parsed_mbr][site].append(0)
+ else:
+ element_index=sites[site].index(element_on_site)
+ site_occupancies[endmember][site][element_index]=proportion_element_on_site
-sites=[[] for i in range(nsites)]
-for i in range(len(base_material)):
- print base_material[i][1]
- for j in range(nsites):
- sites[j].append(j)
+
+array_occupancies=np.empty(shape=(len(base_material),n_elements))
+array_multiplicities=np.empty(shape=(n_elements))
+for endmember in range(len(base_material)):
+ n_element=0
+ for site in range(nsites):
+ for element in range(len(site_occupancies[endmember][site])):
+ array_occupancies[endmember][n_element]=site_occupancies[endmember][site][element]
+ array_multiplicities[n_element]=site_multiplicity[site]
+ n_element=n_element+1
+
+
+# Matrix operation to return site occupancies, given proportions of endmembers
+occupancies=np.dot(endmember_proportions, array_occupancies)
+print 'Site occupancies'
print sites
+print occupancies
+print ''
+
+
+
+# Ideal activities
+ideal_activity=np.empty(shape=(len(base_material)))
+for endmember in range(len(base_material)):
+ ideal_activity[endmember]=1.0
+ normalisation_constant=1.0
+ for element in range(n_elements):
+ if array_occupancies[endmember][element] != 0:
+ #print base_material[endmember][0], element, occupancies[element], array_occupancies[endmember][element]
+ ideal_activity[endmember]=ideal_activity[endmember]*pow(occupancies[element],array_multiplicities[element])
+ normalisation_constant=normalisation_constant/pow(array_occupancies[endmember][element],array_multiplicities[element])
+ ideal_activity[endmember]=normalisation_constant*ideal_activity[endmember]
-#2D vector of endmembers and alphas
-#2D vector of sites and multiplicities [[A,3],[B,2]]
-#2D vector of site occupancies e.g. [[x(Mg,A),x(Fe,A)],[x(Al,B),x(Mg,B),x(Si,B)]]
-#3D vector of site occupancies [[[1.,0.],[1.,0.,0.]],[[0.,1.],[1.,0.,0.]]]
-#2D vector of interaction parameters and temperature and pressure derivatives
+print 'Ideal activities'
+print ideal_activity
-#Operations required:
-#- matrix operation to return site occupancies, given proportions of endmembers
-#- operation to obtain ideal activities
-#- operation to obtain non ideal activities
+# Non ideal activities
+# STILL TO DO
More information about the CIG-COMMITS
mailing list