[cig-commits] [commit] inversion, master, validate_MT_params: Add pressure to the various EOSs that need it (494f1dd)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Fri Dec 12 18:25:04 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 494f1dd70a832763ce1ceee271483625ec3c0b07
Author: ian-r-rose <ian.r.rose at gmail.com>
Date: Mon Sep 1 17:15:34 2014 -0700
Add pressure to the various EOSs that need it
>---------------------------------------------------------------
494f1dd70a832763ce1ceee271483625ec3c0b07
burnman/birch_murnaghan.py | 3 +
burnman/equation_of_state.py | 18 +++
burnman/material.py | 19 ----
burnman/mie_grueneisen_debye.py | 6 +
burnman/processchemistry.py | 207 ++++++++++++++++------------------
burnman/slb.py | 15 +++
burnman/solidsolution.py | 16 +--
burnman/test_process_solidsolution.py | 177 -----------------------------
burnman/test_solidsolution_parsing.py | 66 +++++++++++
9 files changed, 214 insertions(+), 313 deletions(-)
diff --git a/burnman/birch_murnaghan.py b/burnman/birch_murnaghan.py
index e30f8e6..d0fc4f0 100644
--- a/burnman/birch_murnaghan.py
+++ b/burnman/birch_murnaghan.py
@@ -83,6 +83,9 @@ class BirchMurnaghanBase(eos.EquationOfState):
"""
return volume(pressure,params)
+ def pressure(self, temperature, volume, params):
+ return birch_murnaghan(params['V_0']/volume, params)
+
def isothermal_bulk_modulus(self,pressure,temperature, volume, params):
"""
Returns adiabatic bulk modulus [Pa] as a function of pressure [Pa],
diff --git a/burnman/equation_of_state.py b/burnman/equation_of_state.py
index 81c25c2..8b4cfd9 100644
--- a/burnman/equation_of_state.py
+++ b/burnman/equation_of_state.py
@@ -41,6 +41,24 @@ class EquationOfState(object):
"""
raise NotImplementedError("")
+ def pressure(self, temperature, volume, params):
+ """
+ Parameters
+ ----------
+ volume : float
+ Molar volume at which to evaluate the equation of state. [m^3]
+ temperature : float
+ Temperature at which to evaluate the equation of state. [K]
+ params : dictionary
+ Dictionary containing material parameters required by the equation of state.
+
+ Returns
+ -------
+ pressure : float
+ Pressure of the mineral, including cold and thermal parts. [m^3]
+ """
+ raise NotImplementedError("")
+
def density(self, pressure, temperature, params):
"""
Calculate the density of the mineral.
diff --git a/burnman/material.py b/burnman/material.py
index d75bd10..f49d9a9 100644
--- a/burnman/material.py
+++ b/burnman/material.py
@@ -116,22 +116,3 @@ class Material(object):
"""
raise NotImplementedError("need to implement density() in derived class!")
return None
-
- def eos_pressure(self, temperature, volume):
- """
- Returns the pressure of this material given temperature and volume as input parameters. The return value of this function does not depend on the current
- state (temperature, pressure).
-
- Notes
- -----
- Needs to be implemented in derived classes.
-
- Returns
- -------
- eos_pressure : float
- The pressure experienced by this material in Pa
-
- """
- raise NotImplementedError("need to implement eos_pressure() in derived class!")
- return None
-
diff --git a/burnman/mie_grueneisen_debye.py b/burnman/mie_grueneisen_debye.py
index fd606b2..b5b3603 100644
--- a/burnman/mie_grueneisen_debye.py
+++ b/burnman/mie_grueneisen_debye.py
@@ -40,6 +40,12 @@ class MGDBase(eos.EquationOfState):
V = opt.brentq(func, 0.5*params['V_0'], 1.5*params['V_0'])
return V
+ def pressure(self, temperature, volume, params):
+ P = bm.birch_murnaghan(params['V_0']/volume, params) + \
+ self.__thermal_pressure(temperature, volume, params) - \
+ self.__thermal_pressure(300., volume, params)
+ return P
+
def isothermal_bulk_modulus(self, pressure,temperature,volume, params):
"""
Returns isothermal bulk modulus [Pa] as a function of pressure [Pa],
diff --git a/burnman/processchemistry.py b/burnman/processchemistry.py
index b9208b5..fd2aaa6 100644
--- a/burnman/processchemistry.py
+++ b/burnman/processchemistry.py
@@ -10,116 +10,41 @@ import re
import numpy as np
from fractions import Fraction
-def ProcessChemistry(formula):
- filename="data/input_masses/atomic_masses.dat"
- file = open(filename, "r")
- el_name = []
- el_mass = []
-
- i=0
- for line in file:
+def read_masses(filename):
+ lookup=dict()
+ for line in open(filename, "r"):
data="%".join(line.split("%")[:1]).split()
if data != []:
- el_name.append(data[0])
- el_mass.append(float(data[1]))
- file.close()
+ lookup[data[0]]=float(data[1])
+ return lookup
-# Loop over elements
- n=0.
- molar_mass=0.
+def dictionarize_formula(formula):
+ f=dict()
for element in re.findall('[A-Z][^A-Z]*', formula):
list=filter(None, re.split(r'(\d+)', element))
- # Look up number of atoms of element
+ # Look up number of atoms of element
if len(list) == 1:
nel=1.
else:
nel=float(list[1])
- # Increment size of compound
- n=n+nel
+ f[list[0]]=nel
+ return f
- # Find atomic mass of element
- molar_mass=molar_mass+nel*el_mass[el_name.index(list[0])]
+def formula_mass(formula, atomic_masses):
+ mass=0.
+ for element in formula:
+ mass=mass+formula[element]*atomic_masses[element]
+ return mass
- return n, molar_mass
-
-
-def ProcessSolidSolutionChemistry(formula):
- n_sites=formula[0].count('[')
- n_endmembers=len(formula)
- # Check the number of sites is the same for each endmember
- for i in range(n_endmembers):
- assert(formula[i].count('[') == n_sites)
-
- # 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))
- 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'\[', formula[endmember])[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 and multiplicities
- endmember_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])):
- endmember_occupancies[endmember][n_element]=list_occupancies[endmember][site][element]
- site_multiplicities[n_element]=list_multiplicity[site]
- n_element=n_element+1
-
- return n_sites, sites, n_occupancies, endmember_occupancies, site_multiplicities
-
-def CompositionEquality(endmember_formula, solution_formula):
- # Parse endmember formula
- parsed_endmember_formula=[]
- elist=[]
- for enamenumber in re.findall('[A-Z][^A-Z]*', endmember_formula):
- element=filter(None, re.split(r'(\d+)', enamenumber))
- if element[0] not in elist:
- elist.append(element[0])
- parsed_endmember_formula.append([element[0], float(element[1])])
- else:
- element_index=elist.index(element[0])
- parsed_endmember_formula[element_index][1]=parsed_endmember_formula[element_index][1]+float(element[1])
-
- # Parse solution formula
- parsed_solution_formula=[]
-
- s=re.split(r'\[', solution_formula)[1:]
+def dictionarize_site_formula(formula):
+ solution_formulae=dict()
sites=[[] for i in range(len(s))]
list_occupancies=[]
list_multiplicity=np.empty(shape=(len(s)))
n_occupancies=0
- elist=[]
+ s=re.split(r'\[', formula)[1:]
+
for site in range(len(s)):
site_occupancy=re.split(r'\]', s[site])[0]
mult=re.split('[A-Z][^A-Z]*',re.split(r'\]', s[site])[1])[0]
@@ -130,7 +55,7 @@ def CompositionEquality(endmember_formula, solution_formula):
else:
list_multiplicity[site]=mult
- # Loop over elements on site
+ # Loop over elements on a site
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]
@@ -140,23 +65,87 @@ def CompositionEquality(endmember_formula, solution_formula):
else:
proportion_element_on_site=Fraction(proportion_element_on_site[0])
n_element=float(mult)*proportion_element_on_site
- if element_on_site not in elist:
- elist.append(element_on_site)
- parsed_solution_formula.append([element_on_site, n_element])
- else:
- element_index=elist.index(element_on_site)
- parsed_solution_formula[element_index][1]=parsed_solution_formula[element_index][1] + n_element
+ f[element_on_site]=f.get(element_on_site, 0.0) + n_element
+ # Loop over elements not on a site
for enamenumber in re.findall('[A-Z][^A-Z]*', not_in_site):
element=filter(None, re.split(r'(\d+)', enamenumber))
- if element[0] not in elist:
- elist.append(element[0])
- parsed_solution_formula.append([element[0], float(element[1])])
+ f[element[0]]=f.get(element[0], 0.0) + float(element[1])
+
+ return f
+
+
+def ProcessSolidSolutionChemistry(formula):
+ n_sites=formula[0].count('[')
+ n_endmembers=len(formula)
+
+ # Check the number of sites is the same for all endmembers
+ for i in range(n_endmembers):
+ assert(formula[i].count('[') == n_sites)
+
+ solution_formulae=[]
+ sites=[[] for i in range(n_sites)]
+ list_occupancies=[]
+ list_multiplicity=np.empty(shape=(n_sites))
+ n_occupancies=0
+
+ s=re.split(r'\[', formula[i])[1:]
+
+ # Number of unique site occupancies (e.g.. Mg on X etc.)
+ for endmember in range(n_endmembers):
+ solution_formula=dict()
+ list_occupancies.append([[0]*len(sites[site]) for site in range(n_sites)])
+ s=re.split(r'\[', formula[endmember])[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]
+ not_in_site=filter(None, re.split(r'\]', s[site]))[1]
+ not_in_site=not_in_site.replace(mult, '', 1)
+ if mult == '':
+ list_multiplicity[site]=1.0
else:
- element_index=elist.index(element[0])
- parsed_solution_formula[element_index][1]=parsed_solution_formula[element_index][1]+float(element[1])
+ list_multiplicity[site]=mult
+
+ # Loop over elements on a site
+ 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])
+ solution_formula[element_on_site]=solution_formula.get(element_on_site, 0.0) + float(mult)*proportion_element_on_site
+
+ 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
+
+ # Loop over elements not on a site
+ for enamenumber in re.findall('[A-Z][^A-Z]*', not_in_site):
+ element=filter(None, re.split(r'(\d+)', enamenumber))
+ solution_formula[element[0]]=solution_formula.get(element[0], 0.0) + float(element[1])
+
+ solution_formulae.append(solution_formula)
+
+ # Site occupancies and multiplicities
+ endmember_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])):
+ endmember_occupancies[endmember][n_element]=list_occupancies[endmember][site][element]
+ site_multiplicities[n_element]=list_multiplicity[site]
+ n_element=n_element+1
+
- parsed_endmember_formula=sorted(parsed_endmember_formula,key=lambda l:l[0])
- parsed_solution_formula=sorted(parsed_solution_formula,key=lambda l:l[0])
+ return solution_formulae, n_sites, sites, n_occupancies, endmember_occupancies, site_multiplicities
- return parsed_endmember_formula == parsed_solution_formula
diff --git a/burnman/slb.py b/burnman/slb.py
index 8bfada2..358d182 100644
--- a/burnman/slb.py
+++ b/burnman/slb.py
@@ -100,7 +100,22 @@ class SLBBase(eos.EquationOfState):
warnings.warn("May be outside the range of validity for EOS")
return sol[0]
+ def pressure( self, temperature, volume, params):
+ """
+ Returns the pressure of the mineral at a given temperature and volume [Pa]
+ """
+ debye_T = self.__debye_temperature(params['V_0']/volume, params)
+ gr = self.grueneisen_parameter(0.0, temperature, volume, params) #does not depend on pressure
+ E_th = debye.thermal_energy(temperature, debye_T, params['n'])
+ E_th_ref = debye.thermal_energy(300., debye_T, params['n']) #thermal energy at reference temperature
+
+ b_iikk= 9.*params['K_0'] # EQ 28
+ b_iikkmm= 27.*params['K_0']*(params['Kprime_0']-4.) # EQ 29
+ f = 0.5*(pow(params['V_0']/volume,2./3.)-1.) # EQ 24
+ P = (1./3.)*(pow(1.+2.*f,5./2.))*((b_iikk*f) \
+ +(0.5*b_iikkmm*pow(f,2.))) + gr*(E_th - E_th_ref)/volume #EQ 21
+ return P
def grueneisen_parameter(self, pressure, temperature, volume, params):
"""
diff --git a/burnman/solidsolution.py b/burnman/solidsolution.py
index 92096ca..8828394 100644
--- a/burnman/solidsolution.py
+++ b/burnman/solidsolution.py
@@ -5,7 +5,6 @@
import numpy as np
from burnman.mineral import Mineral
from burnman.processchemistry import ProcessSolidSolutionChemistry
-from burnman.processchemistry import CompositionEquality
import warnings
R = 8.3145 # J/K/mol
@@ -51,17 +50,18 @@ class SolidSolution(Mineral):
self.excess_entropy=interaction_parameter[1]
self.excess_volume=interaction_parameter[2]
+ # Process solid solution chemistry
+ self.solution_formulae, self.n_sites, self.sites, self.n_occupancies, self.endmember_occupancies, self.site_multiplicities = ProcessSolidSolutionChemistry([base_material[i][1] for i in range(self.n_endmembers)])
+
# Check the chemical composition of each endmember is consistent with the dataset
for i in range(self.n_endmembers):
- endmember_formula=base_material[i][0].params['formula']
- solution_formula=base_material[i][1]
- if CompositionEquality(endmember_formula, solution_formula) != True:
- print 'Formula of endmember', base_material[i][0].params['name'], 'does not agree with formula in the', self.name, 'SolidSolution model'
+ endmember_formula=self.base_material[i][0].params['formula']
+ solution_formula=self.solution_formulae[i]
+ if endmember_formula != solution_formula:
+ print 'Formula of endmember', self.base_material[i][0].params['name'], 'does not agree with formula in the', self.name, 'SolidSolution model'
+ print endmember_formula, self.solution_formulae[i]
exit()
- # Process solid solution chemistry
- self.n_sites, self.sites, self.n_occupancies, self.endmember_occupancies, self.site_multiplicities = ProcessSolidSolutionChemistry([base_material[i][1] for i in range(self.n_endmembers)])
-
# Create array of van Laar parameters
self.alpha=np.array([base_material[i][2] for i in range(self.n_endmembers)])
diff --git a/burnman/test_process_solidsolution.py b/burnman/test_process_solidsolution.py
deleted file mode 100644
index 92fbbde..0000000
--- a/burnman/test_process_solidsolution.py
+++ /dev/null
@@ -1,177 +0,0 @@
-# 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()
diff --git a/burnman/test_solidsolution_parsing.py b/burnman/test_solidsolution_parsing.py
new file mode 100644
index 0000000..ba7b863
--- /dev/null
+++ b/burnman/test_solidsolution_parsing.py
@@ -0,0 +1,66 @@
+# 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
+from processchemistry import ProcessSolidSolutionChemistry
+
+# 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 ])
+
+
+
+solution_formulae, n_sites, sites, n_occupancies, endmember_occupancies, site_multiplicities = ProcessSolidSolutionChemistry([base_material[i][1] for i in range(len(base_material))])
+
+
+print solution_formulae
+print n_sites
+print n_occupancies
+print sites
+print endmember_occupancies
+print site_multiplicities
+
More information about the CIG-COMMITS
mailing list