From 4918309a384bce8b7dd6e5d3c9144cf36459c3df Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Wed, 28 Oct 2015 16:18:24 +0100 Subject: [PATCH 1/2] Added activity functions to solid solution, updated example document --- burnman/solidsolution.py | 30 ++++++++++--- examples/example_solid_solution.py | 68 ++++++++++++++++++++++++------ 2 files changed, 81 insertions(+), 17 deletions(-) diff --git a/burnman/solidsolution.py b/burnman/solidsolution.py index 0e353f54f..7d736423a 100644 --- a/burnman/solidsolution.py +++ b/burnman/solidsolution.py @@ -10,8 +10,7 @@ from .mineral import Mineral, material_property from .solutionmodel import * from .averaging_schemes import reuss_average_function - - +from . import constants class SolidSolution(Mineral): """ @@ -122,9 +121,30 @@ def set_state(self, pressure, temperature): Mineral.set_state(self, pressure, temperature) for i in range(self.n_endmembers): self.endmembers[i][0].set_state(pressure, temperature) + + @material_property + def activities(self): + """ + Returns a list of endmember activities + """ + if self.temperature > 1.e-10: + return np.exp(self.excess_partial_gibbs/(constants.gas_constant*self.temperature)) + else: + return None + @material_property + def activity_coefficients(self): + """ + Returns a list of endmember activity coefficients (the non-ideal part of the activities + """ + if self.temperature > 1.e-10: + return np.exp(self.excess_partial_gibbs/(constants.gas_constant*self.temperature) + - IdealSolution._log_ideal_activities(self.solution_model, self.molar_fractions)) + else: + return None + @material_property def internal_energy(self): """ @@ -152,8 +172,7 @@ def partial_gibbs(self): """ return np.array([self.endmembers[i][0].gibbs for i in range(self.n_endmembers)]) + self.excess_partial_gibbs - - + @material_property def excess_gibbs(self): """ @@ -273,7 +292,7 @@ def adiabatic_bulk_modulus(self): return self.isothermal_bulk_modulus else: return self.isothermal_bulk_modulus*self.heat_capacity_p/self.heat_capacity_v - + @material_property def isothermal_compressibility(self): @@ -368,3 +387,4 @@ def heat_capacity_p(self): Aliased with self.C_p """ return sum([ self.endmembers[i][0].heat_capacity_p * self.molar_fractions[i] for i in range(self.n_endmembers) ]) + diff --git a/examples/example_solid_solution.py b/examples/example_solid_solution.py index 3f1ffe9e2..850818845 100644 --- a/examples/example_solid_solution.py +++ b/examples/example_solid_solution.py @@ -87,18 +87,38 @@ def __init__(self, molar_fractions=None): comp = np.linspace(0.001, 0.999, 100) g1_gibbs = np.empty_like(comp) g1_excess_gibbs = np.empty_like(comp) + g1_pyrope_activity = np.empty_like(comp) + g1_almandine_activity = np.empty_like(comp) + g1_pyrope_gamma = np.empty_like(comp) + g1_almandine_gamma = np.empty_like(comp) for i,c in enumerate(comp): molar_fractions=[1.0-c, c] g1.set_composition(molar_fractions) g1.set_state(P,T) g1_gibbs[i] = g1.gibbs g1_excess_gibbs[i] = g1.excess_gibbs + g1_pyrope_activity[i] = g1.activities[0] + g1_almandine_activity[i] = g1.activities[1] + g1_pyrope_gamma[i] = g1.activity_coefficients[0] + g1_almandine_gamma[i] = g1.activity_coefficients[1] plt.plot( [0., 1.], [py_gibbs/1000., alm_gibbs/1000.], 'b-', linewidth=1., label='Mechanical mixing') plt.plot( comp, g1_gibbs/1000., 'r-', linewidth=1., label='Ideal solid solution') - plt.title("Pyrope-almandine join") + plt.title("Ideal pyrope-almandine join") plt.ylabel("Gibbs free energy of solution (kJ/mol)") - plt.xlabel("Pyrope fraction") + plt.xlabel("Almandine fraction") + plt.legend(loc='lower right') + plt.show() + + plt.plot( comp, g1_pyrope_activity, 'r-', linewidth=1., label='Pyrope activity') + plt.plot( comp, g1_almandine_activity, 'b-', linewidth=1., label='Almandine activity') + + plt.plot( comp, g1_pyrope_gamma, 'r-.', linewidth=1., label='Pyrope activity coefficient') + plt.plot( comp, g1_almandine_gamma, 'b--', linewidth=1., label='Almandine activity coefficient') + plt.title("Ideal pyrope-almandine join") + plt.ylim(0.0, 1.01) + plt.ylabel("Activities") + plt.xlabel("Almandine fraction") plt.legend(loc='lower right') plt.show() @@ -128,22 +148,44 @@ def __init__(self, molar_fractions=None): burnman.SolidSolution.__init__(self, molar_fractions) + g2=mg_fe_garnet_2() g2_excess_gibbs = np.empty_like(comp) + g2_pyrope_activity = np.empty_like(comp) + g2_almandine_activity = np.empty_like(comp) + g2_pyrope_gamma = np.empty_like(comp) + g2_almandine_gamma = np.empty_like(comp) for i,c in enumerate(comp): - molar_fractions=[1.0-c, c] + molar_fractions=[1.-c, c] g2.set_composition(molar_fractions) g2.set_state(P,T) g2_excess_gibbs[i] = g2.excess_gibbs + g2_pyrope_activity[i] = g2.activities[0] + g2_almandine_activity[i] = g2.activities[1] + g2_pyrope_gamma[i] = g2.activity_coefficients[0] + g2_almandine_gamma[i] = g2.activity_coefficients[1] plt.plot( comp, g1_excess_gibbs, 'r-', linewidth=1., label='Ideal solution') plt.plot( comp, g2_excess_gibbs, 'g-', linewidth=1., label='Symmetric solution, 2.5 kJ/mol') plt.title("Pyrope-almandine join (model comparison)") plt.ylabel("Excess gibbs free energy of solution (J/mol)") - plt.xlabel("Pyrope fraction") + plt.xlabel("Almandine fraction") plt.legend(loc='upper left') plt.show() + + plt.plot( comp, g2_pyrope_activity, 'r-', linewidth=1., label='Pyrope activity') + plt.plot( comp, g2_almandine_activity, 'b-', linewidth=1., label='Almandine activity') + + plt.plot( comp, g2_pyrope_gamma, 'r-.', linewidth=1., label='Pyrope activity coefficient') + plt.plot( comp, g2_almandine_gamma, 'b--', linewidth=1., label='Almandine activity coefficient') + plt.title("Ideal pyrope-almandine join") + plt.ylabel("Activities") + plt.xlabel("Almandine fraction") + plt.legend(loc='lower right') + plt.show() + + ''' Adding more endmembers is very straightforward. Interaction terms must be added in the order [[12,13,14,...],[23,24,...],[34,...],...] @@ -187,7 +229,7 @@ def __init__(self, molar_fractions=None): plt.plot( comp, g3_excess_entropy, 'r-', linewidth=1., label='Excess entropy') plt.title("Pyrope-majorite join") plt.ylabel("Entropy (J/K/mol)") - plt.xlabel("Pyrope fraction") + plt.xlabel("Majorite fraction") plt.legend(loc='upper left') plt.show() @@ -231,7 +273,7 @@ def __init__(self, molar_fractions=None): plt.plot( comp, g4_excess_gibbs_1200, 'b-', linewidth=1., label='1200 K') plt.title("Pyrope-grossular join (asymmetric model)") plt.ylabel("Excess gibbs free energy of solution (J/mol)") - plt.xlabel("Pyrope fraction") + plt.xlabel("Grossular fraction") plt.legend(loc='lower left') plt.show() @@ -260,17 +302,19 @@ def __init__(self, molar_fractions=None): burnman.SolidSolution.__init__(self, molar_fractions) g5=mg_fe_ca_garnet_Ganguly() - g5_excess_enthalpy = np.empty_like(comp) - + g5_excess_gibbs = np.empty_like(comp) for i,c in enumerate(comp): - molar_fractions=[1.0-c, 0., c, 0.] + molar_fractions=[1.-c, 0., c, 0.] g5.set_composition(molar_fractions) g5.set_state(1.e5,298.15) - g5_excess_enthalpy[i] = g5.excess_enthalpy + g5_excess_gibbs[i] = g5.excess_gibbs + + plt.plot( comp, g5_excess_gibbs/3., 'r-', linewidth=1., label='Py-Gr excess enthalpy (J/cation-mole)') - plt.plot( comp, g5_excess_enthalpy/3., 'r-', linewidth=1., label='Py-Gr excess enthalpy (J/cation-mole)') plt.title("Asymmetric py-gr join (Ganguly et al., 1996; Figure 5)") plt.ylabel("Excess enthalpy of solution (J/cation-mol)") - plt.xlabel("Pyrope fraction") + plt.xlabel("XCa") plt.legend(loc='lower left') plt.show() + + From a585b52b29f91930e161f9c12ec7ae75c0e7c544 Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Wed, 24 Feb 2016 14:47:14 +0000 Subject: [PATCH 2/2] Rejigged activity formulation for efficiency, added tests --- burnman/solidsolution.py | 23 +++++++--------------- burnman/solutionmodel.py | 27 +++++++++++++++++++++++++- examples/example_solid_solution.py | 2 +- tests/test_solidsolution.py | 31 ++++++++++++++++++++++++++++-- 4 files changed, 63 insertions(+), 20 deletions(-) diff --git a/burnman/solidsolution.py b/burnman/solidsolution.py index 7d736423a..bca2bc4b3 100644 --- a/burnman/solidsolution.py +++ b/burnman/solidsolution.py @@ -126,23 +126,16 @@ def set_state(self, pressure, temperature): @material_property def activities(self): """ - Returns a list of endmember activities + Returns a list of endmember activities [unitless] """ - if self.temperature > 1.e-10: - return np.exp(self.excess_partial_gibbs/(constants.gas_constant*self.temperature)) - else: - return None + return self.solution_model.activities( self.pressure, self.temperature, self.molar_fractions) @material_property def activity_coefficients(self): """ - Returns a list of endmember activity coefficients (the non-ideal part of the activities + Returns a list of endmember activity coefficients (gamma = activity / ideal activity) [unitless] """ - if self.temperature > 1.e-10: - return np.exp(self.excess_partial_gibbs/(constants.gas_constant*self.temperature) - - IdealSolution._log_ideal_activities(self.solution_model, self.molar_fractions)) - else: - return None + return self.solution_model.activity_coefficients( self.pressure, self.temperature, self.molar_fractions) @material_property @@ -151,8 +144,7 @@ def internal_energy(self): Returns internal energy of the mineral [J] Aliased with self.energy """ - raise NotImplementedError("need to implement internal_energy() for solid solution!") - return None + return self.molar_helmholtz - self.pressure*self.molar_volume @material_property @@ -197,8 +189,7 @@ def molar_helmholtz(self): Returns Helmholtz free energy of the solid solution [J] Aliased with self.helmholtz """ - raise NotImplementedError("need to implement molar_helmholtz() for solid solution!") - return None + return self.molar_gibbs + self.temperature*self.molar_entropy @material_property @@ -248,7 +239,7 @@ def excess_entropy(self): @material_property def molar_entropy(self): """ - Returns enthalpy of the solid solution [J] + Returns entropy of the solid solution [J] Aliased with self.S """ return sum([ self.endmembers[i][0].S * self.molar_fractions[i] for i in range(self.n_endmembers) ]) + self.excess_entropy diff --git a/burnman/solutionmodel.py b/burnman/solutionmodel.py index 41ed172e9..3167ca9d2 100644 --- a/burnman/solutionmodel.py +++ b/burnman/solutionmodel.py @@ -229,6 +229,13 @@ def _ideal_activities (self, molar_fractions): activities[e]=normalisation_constant*activities[e] return activities + def activity_coefficients(self, pressure, temperature, molar_fractions): + return np.ones_like(molar_fractions) + + def activities (self, pressure, temperature, molar_fractions): + return self._ideal_activities(molar_fractions) + + class AsymmetricRegularSolution (IdealSolution): """ @@ -314,7 +321,15 @@ def excess_enthalpy(self, pressure, temperature, molar_fractions): H_excess=np.dot(self.alpha.T,molar_fractions)*np.dot(phi.T,np.dot(self.Wh,phi)) return H_excess + pressure*self.excess_volume ( pressure, temperature, molar_fractions ) - + def activity_coefficients(self, pressure, temperature, molar_fractions): + if temperature > 1.e-10: + return np.exp(self._non_ideal_excess_partial_gibbs(pressure, temperature, molar_fractions)/(constants.gas_constant*temperature)) + else: + raise Exception("Activity coefficients not defined at 0 K.") + + def activities(self, pressure, temperature, molar_fractions): + return IdealSolution.activities(self, pressure, temperature, molar_fractions) * self.activity_coefficients(pressure, temperature, molar_fractions) + class SymmetricRegularSolution (AsymmetricRegularSolution): """ Solution model implementing the symmetric regular solution model @@ -401,3 +416,13 @@ def excess_entropy(self, pressure, temperature, molar_fractions): def excess_enthalpy(self, pressure, temperature, molar_fractions): H_excess=np.dot(molar_fractions, self._non_ideal_function(self.Wh, molar_fractions)) return H_excess + pressure*self.excess_volume (pressure, temperature, molar_fractions) + + def activity_coefficients(self, pressure, temperature, molar_fractions): + if temperature > 1.e-10: + return np.exp(self._non_ideal_excess_partial_gibbs(pressure, temperature, molar_fractions)/(constants.gas_constant*temperature)) + else: + raise Exception("Activity coefficients not defined at 0 K.") + + def activities(self, pressure, temperature, molar_fractions): + return IdealSolution.activities(self, pressure, temperature, molar_fractions) * self.activity_coefficients(pressure, temperature, molar_fractions) + diff --git a/examples/example_solid_solution.py b/examples/example_solid_solution.py index 850818845..fad5a148b 100644 --- a/examples/example_solid_solution.py +++ b/examples/example_solid_solution.py @@ -179,7 +179,7 @@ def __init__(self, molar_fractions=None): plt.plot( comp, g2_pyrope_gamma, 'r-.', linewidth=1., label='Pyrope activity coefficient') plt.plot( comp, g2_almandine_gamma, 'b--', linewidth=1., label='Almandine activity coefficient') - plt.title("Ideal pyrope-almandine join") + plt.title("Non-ideal pyrope-almandine join") plt.ylabel("Activities") plt.xlabel("Almandine fraction") plt.legend(loc='lower right') diff --git a/tests/test_solidsolution.py b/tests/test_solidsolution.py index e4ec3340e..4d9737d28 100644 --- a/tests/test_solidsolution.py +++ b/tests/test_solidsolution.py @@ -70,6 +70,15 @@ def __init__(self, molar_fractions=None): burnman.SolidSolution.__init__(self, molar_fractions) +# Ideal solid solution +class olivine_ideal_ss(burnman.SolidSolution): + def __init__(self, molar_fractions=None): + self.name='Fo-Fo solid solution' + self.type='ideal' + self.endmembers = [[forsterite(), '[Mg]2SiO4'], [fayalite(), '[Fe]2SiO4']] + + burnman.SolidSolution.__init__(self, molar_fractions) + # Olivine solid solution class olivine_ss(burnman.SolidSolution): def __init__(self, molar_fractions=None): @@ -86,8 +95,8 @@ def __init__(self, molar_fractions=None): # Name self.name='orthopyroxene' self.type='symmetric' - self.endmembers = [[forsterite(), 'Mg[Mg]Si2O6'], [forsterite(), '[Mg1/2Al1/2]2AlSiO6']] - self.enthalpy_interaction=[[10.0e3]] + self.endmembers = [[forsterite(), '[Mg][Mg]Si2O6'], [forsterite(), '[Mg1/2Al1/2][Mg1/2Al1/2]AlSiO6']] + self.enthalpy_interaction=[[burnman.constants.gas_constant*1.0e3]] burnman.SolidSolution.__init__(self, molar_fractions) @@ -217,6 +226,24 @@ def test_subregular(self): self.assertArraysAlmostEqual(ss0.excess_partial_gibbs, ss1.excess_partial_gibbs) + def test_activities_ideal(self): + ol = olivine_ideal_ss() + ol.set_composition( np.array([0.5, 0.5]) ) + ol.set_state(1.e5,1000.) + self.assertArraysAlmostEqual(ol.activities, [0.25, 0.25]) + + def test_activity_coefficients_ideal(self): + ol = olivine_ideal_ss() + ol.set_composition( np.array([0.5, 0.5]) ) + ol.set_state(1.e5,1000.) + self.assertArraysAlmostEqual(ol.activity_coefficients, [1., 1.]) + + def test_activity_coefficients_non_ideal(self): + opx = orthopyroxene() + opx.set_composition( np.array([0.0, 1.0]) ) + opx.set_state(1.e5,1000.) + self.assertArraysAlmostEqual(opx.activity_coefficients, [np.exp(1.), 1.]) + if __name__ == '__main__': unittest.main()