[cig-commits] [commit] add_gibbs_energy: Added more attributes to the solidsolution class (7431641)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Dec 11 17:11:06 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_gibbs_energy
Link : https://github.com/geodynamics/burnman/compare/0000000000000000000000000000000000000000...2148b324d3e8aa7b527f831eb397590942563008
>---------------------------------------------------------------
commit 7431641be9e85360416ffebff1a8b3ddec5e81bf
Author: Bob Myhill <myhill.bob at gmail.com>
Date: Sun Aug 31 03:13:33 2014 +0200
Added more attributes to the solidsolution class
>---------------------------------------------------------------
7431641be9e85360416ffebff1a8b3ddec5e81bf
burnman/solidsolution.py | 16 +++++-----
burnman/test_solidsolution.py | 71 +++++++++++++++++++++++++++++++++++++++++++
2 files changed, 79 insertions(+), 8 deletions(-)
diff --git a/burnman/solidsolution.py b/burnman/solidsolution.py
index 0001d1c..375057d 100644
--- a/burnman/solidsolution.py
+++ b/burnman/solidsolution.py
@@ -68,12 +68,12 @@ class SolidSolution(Mineral):
# NOT IMPLEMENTED YET
# Number of unique site occupancies (e.g.. Mg on X etc.)
- sites=[[] for i in range(self.n_sites)]
+ self.sites=[[] for i in range(self.n_sites)]
list_occupancies=[]
list_multiplicity=np.empty(shape=(self.n_sites))
self.n_occupancies=0
for endmember in range(self.n_endmembers):
- list_occupancies.append([[0]*len(sites[site]) for site in range(self.n_sites)])
+ list_occupancies.append([[0]*len(self.sites[site]) for site in range(self.n_sites)])
s=re.split(r'\[', base_material[endmember][1])[1:]
for site in range(self.n_sites):
site_occupancy=re.split(r'\]', s[site])[0]
@@ -91,14 +91,14 @@ class SolidSolution(Mineral):
else:
proportion_element_on_site=Fraction(proportion_element_on_site[0])
- if element_on_site not in sites[site]:
+ if element_on_site not in self.sites[site]:
self.n_occupancies=self.n_occupancies+1
- sites[site].append(element_on_site)
- element_index=sites[site].index(element_on_site)
+ self.sites[site].append(element_on_site)
+ element_index=self.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)
+ element_index=self.sites[site].index(element_on_site)
list_occupancies[endmember][site][element_index]=proportion_element_on_site
# Site occupancies and multiplicities
@@ -134,14 +134,14 @@ class SolidSolution(Mineral):
self.molar_fraction=molar_fraction
# Ideal activities
- occupancies=np.dot(self.molar_fraction, self.site_occupancies)
+ self.occupancies=np.dot(self.molar_fraction, self.site_occupancies)
self.ideal_activity=np.empty(shape=(self.n_endmembers))
for endmember in range(self.n_endmembers):
self.ideal_activity[endmember]=1.0
normalisation_constant=1.0
for element in range(self.n_occupancies):
if self.site_occupancies[endmember][element] != 0:
- self.ideal_activity[endmember]=self.ideal_activity[endmember]*pow(occupancies[element],self.site_multiplicities[element])
+ self.ideal_activity[endmember]=self.ideal_activity[endmember]*pow(self.occupancies[element],self.site_multiplicities[element])
normalisation_constant=normalisation_constant/pow(self.site_occupancies[endmember][element],self.site_multiplicities[element])
self.ideal_activity[endmember]=normalisation_constant*self.ideal_activity[endmember]
diff --git a/burnman/test_solidsolution.py b/burnman/test_solidsolution.py
new file mode 100644
index 0000000..09e9056
--- /dev/null
+++ b/burnman/test_solidsolution.py
@@ -0,0 +1,71 @@
+# 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.
+
+
+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 import minerals
+
+
+garnet=minerals.HP_2011_ds62.garnet()
+composition=np.array([ 0.5, 0.2, 0.1, 0.2 ])
+garnet.set_composition(composition)
+
+print 'Molar fraction'
+print garnet.molar_fraction
+print ''
+print 'Site occupancies'
+print garnet.sites
+print garnet.occupancies
+print ''
+print 'Ideal activities'
+print garnet.ideal_activity
+print ''
+print 'Volume excess'
+print garnet.V_excess
+print ''
+
+# Excess volumes for the pyrope-grossular join
+n=100
+pyrope_proportion= np.empty(shape=(n+1))
+garnet_excess_volume= np.empty(shape=(n+1))
+for i in range(n+1):
+ pyrope_proportion[i]=float(i)/n
+ composition=([ pyrope_proportion[i], 0.0, 1.-pyrope_proportion[i], 0.0 ])
+ garnet.set_composition(composition)
+ garnet_excess_volume[i]=garnet.V_excess
+
+pressure=1.e9 # Pa
+temperature=573.15 # K
+composition=[0.9, 0.0, 0.1, 0.0]
+garnet.set_state(pressure, temperature, composition)
+
+# Excess gibbs for the pyrope-grossular join
+n=100
+pyrope_proportion= np.empty(shape=(n+1))
+garnet_excess_gibbs= np.empty(shape=(n+1))
+for i in range(n+1):
+ pyrope_proportion[i]=float(i)/n
+ composition=([ pyrope_proportion[i], 0.0, 1.-pyrope_proportion[i], 0.0 ])
+ garnet.set_state(pressure, temperature, composition)
+ garnet_excess_gibbs[i]=garnet.gibbs_excess
+
+
+import matplotlib.pyplot as plt
+plt.subplot(1,2,1)
+plt.plot(pyrope_proportion,garnet_excess_volume,color='r',linestyle='-',marker='o',markerfacecolor='r',markersize=0)
+plt.xlim(min(pyrope_proportion),max(pyrope_proportion))
+plt.xlabel("p(pyrope)")
+plt.title("V excess (m^3/mol) \nfor pyrope-grossular garnets")
+
+plt.subplot(1,2,2)
+plt.plot(pyrope_proportion,garnet_excess_gibbs,color='r',linestyle='-',marker='o',markerfacecolor='r',markersize=0)
+plt.xlim(min(pyrope_proportion),max(pyrope_proportion))
+plt.xlabel("p(pyrope)")
+plt.title("Excess Gibbs (J/mol) for pyrope-grossular garnets\n"+str(pressure/1.e9)+" GPa, "+str(temperature)+" K")
+plt.show()
+
More information about the CIG-COMMITS
mailing list