Merge pull request #710 from bobmyhill/equilibrate_relaxed_composite

Equilibrate relaxed composite

Author: bobmyhill

burnman/classes/material.py

@@ -55,7 +55,8 @@ class Material(object):
     """
 
     def __init__(self):
-        self._number_of_moles = 1.0
+        if not hasattr(self, "number_of_moles"):
+            self.number_of_moles = 1.0
         if not hasattr(self, "name"):
             # if a derived class decides to set .name before calling this
             # constructor (I am looking at you, SLB_2011.py!), do not
@@ -212,6 +213,7 @@ def number_of_moles(self, value):
         :param value: Number of moles.
         :type value: float
         """
+        self.reset()
         if value < 0:
             raise ValueError("The number of moles cannot be negative")
         self._number_of_moles = value

burnman/classes/relaxedcomposite.py

@@ -108,8 +108,17 @@ def set_state(self, pressure, temperature, bulk_composition=None):
         )
         assert sol.success, "Equilibration failed in RelaxedComposite.set_state()"
 
+        self.copy_state_from_unrelaxed()
+
+    def copy_state_from_unrelaxed(self):
+        """
+        Copies the state (P, T, molar_fractions) from the unrelaxed
+        composite to the relaxed composite. This can be useful
+        if the unrelaxed composite has been modified externally
+        and the relaxed composite needs to be updated to match.
+        """
         # Set the relaxed composite to the equilibrated state
-        Composite.set_state(self, pressure, temperature)
+        Composite.set_state(self, self.unrelaxed.pressure, self.unrelaxed.temperature)
         Composite.set_fractions(self, self.unrelaxed.molar_fractions)
         self.number_of_moles = self.unrelaxed.number_of_moles
 

burnman/tools/equilibration.py

@@ -159,8 +159,6 @@ def set_compositions_and_state_from_parameters(assemblage, parameters):
     :param parameters: The current parameter values.
     :type parameters: numpy.array
     """
-
-    assemblage.set_state(parameters[0] * P_scaling, parameters[1] * T_scaling)
     i = 2
     phase_amounts = np.zeros(len(assemblage.phases))
     for phase_idx, ph in enumerate(assemblage.phases):
@@ -179,6 +177,15 @@ def set_compositions_and_state_from_parameters(assemblage, parameters):
     phase_amounts = np.abs(phase_amounts)
     assemblage.number_of_moles = sum(phase_amounts)
     assemblage.set_fractions(phase_amounts / assemblage.number_of_moles)
+
+    # We set_state last because setting the state may depend on the
+    # current composition of the assemblage. This would occur, for example,
+    # in the case of relaxed composites, where the relaxed state depends
+    # on the current bulk composition.
+    # Note that in the case of relaxed composites, the equilibrate function
+    # actually equilibrates the unrelaxed assemblage, and then copies the
+    # state to the relaxed assemblage afterwards.
+    assemblage.set_state(parameters[0] * P_scaling, parameters[1] * T_scaling)
     return None
 
 
@@ -922,7 +929,7 @@ def equilibrate(
     :type composition: dict
 
     :param assemblage: The assemblage to be equilibrated.
-    :type assemblage: :class:`burnman.Composite`
+    :type assemblage: :class:`burnman.Composite` or :class:`burnman.RelaxedComposite`
 
     :param equality_constraints: The list of equality constraints. See above
         for valid formats.
@@ -978,6 +985,11 @@ def equilibrate(
     outer_logger_level = logger.level
     logger.setLevel(logging.INFO if verbose else outer_logger_level)
 
+    relaxed_assemblage = None
+    if hasattr(assemblage, "unrelaxed"):
+        relaxed_assemblage = assemblage
+        assemblage = assemblage.unrelaxed
+
     for ph in assemblage.phases:
         if isinstance(ph, Solution) and not hasattr(ph, "molar_fractions"):
             raise Exception(
@@ -1155,7 +1167,11 @@ def lambda_bounds_fn(dx, x):
                 )
             sol = sol3
 
-        if sol.success and len(assemblage.reaction_affinities) > 0.0:
+        if relaxed_assemblage is not None:
+            relaxed_assemblage.copy_state_from_unrelaxed()
+            assemblage = relaxed_assemblage
+
+        if sol.success and len(assemblage.reaction_affinities) > 0:
             maxres = np.max(np.abs(assemblage.reaction_affinities)) + 1.0e-5
             assemblage.equilibrium_tolerance = maxres
 

tests/test_composite.py

@@ -6,6 +6,43 @@
 from burnman import minerals, Composite, RelaxedComposite
 from burnman.minerals import HP_2011_ds62, mp50MnNCKFMASHTO
 from burnman.tools.eos import check_eos_consistency
+from burnman.tools.equilibration import equilibrate
+
+properties_extensive = [
+    "internal_energy",
+    "helmholtz",
+    "gibbs",
+    "enthalpy",
+    "entropy",
+    "volume",
+    "heat_capacity_p",
+    "heat_capacity_v",
+    "mass",
+    "H",
+    "S",
+    "V",
+    "C_p",
+    "C_v",
+]
+
+properties_intensive = [
+    "molar_internal_energy",
+    "molar_helmholtz",
+    "molar_gibbs",
+    "molar_enthalpy",
+    "molar_entropy",
+    "molar_volume",
+    "molar_heat_capacity_p",
+    "molar_heat_capacity_v",
+    "isothermal_bulk_modulus_reuss",
+    "isentropic_bulk_modulus_reuss",
+    "isothermal_compressibility_reuss",
+    "isentropic_compressibility_reuss",
+    "grueneisen_parameter",
+    "thermal_expansivity",
+    "molar_mass",
+    "density",
+]
 
 
 def setup_assemblage():
@@ -523,6 +560,37 @@ def test_single_phase_composite_aliases(self):
         for i in range(len(properties)):
             self.assertFloatEqual(rock_prps[i], alias_prps[i])
 
+    def test_number_of_moles_composite(self):
+        min1 = minerals.SLB_2011.mg_perovskite()
+        min2 = minerals.SLB_2011.periclase()
+        molar_fractions = [0.5, 0.5]
+        c1 = burnman.Composite([min1, min2], molar_fractions, fraction_type="molar")
+        c2 = burnman.Composite([min1, min2], molar_fractions, fraction_type="molar")
+        c1.set_state(40.0e9, 2000.0)
+        c2.set_state(40.0e9, 2000.0)
+
+        for prop in properties_extensive:
+            val1 = getattr(c1, prop)
+            val2 = getattr(c2, prop)
+            self.assertFloatEqual(val1, val2)
+
+        for prop in properties_intensive:
+            val1 = getattr(c1, prop)
+            val2 = getattr(c2, prop)
+            self.assertFloatEqual(val1, val2)
+
+        c2.number_of_moles = 2.0
+
+        for prop in properties_extensive:
+            val1 = getattr(c1, prop)
+            val2 = getattr(c2, prop)
+            self.assertFloatEqual(val1, val2 / 2.0)
+
+        for prop in properties_intensive:
+            val1 = getattr(c1, prop)
+            val2 = getattr(c2, prop)
+            self.assertFloatEqual(val1, val2)
+
     def test_relaxed_composite(self):
         assemblage = setup_assemblage()
         old_formula = assemblage.formula.copy()
@@ -595,47 +663,11 @@ def test_property_extensivity_composite(self):
         c1.set_state(40.0e9, 2000.0)
         c2.set_state(40.0e9, 2000.0)
 
-        properties_extensive = [
-            "internal_energy",
-            "helmholtz",
-            "gibbs",
-            "enthalpy",
-            "entropy",
-            "volume",
-            "heat_capacity_p",
-            "heat_capacity_v",
-            "mass",
-            "V",
-            "S",
-            "H",
-            "C_v",
-            "C_p",
-        ]
-
         for prop in properties_extensive:
             val1 = getattr(c1, prop)
             val2 = getattr(c2, prop)
             self.assertFloatEqual(val1, val2 / 2.0)
 
-        properties_intensive = [
-            "molar_internal_energy",
-            "molar_helmholtz",
-            "molar_gibbs",
-            "molar_enthalpy",
-            "molar_entropy",
-            "molar_volume",
-            "molar_heat_capacity_p",
-            "molar_heat_capacity_v",
-            "isothermal_bulk_modulus_reuss",
-            "isentropic_bulk_modulus_reuss",
-            "isothermal_compressibility_reuss",
-            "isentropic_compressibility_reuss",
-            "grueneisen_parameter",
-            "thermal_expansivity",
-            "molar_mass",
-            "density",
-        ]
-
         for prop in properties_intensive:
             val1 = getattr(c1, prop)
             val2 = getattr(c2, prop)
@@ -663,18 +695,28 @@ def test_formula_extensivity_relaxed_composite(self):
         formula = c1_unrelaxed.formula.copy()
         formula2 = c2_unrelaxed.formula.copy()
 
-        self.assertFloatEqual(formula["Mg"] * 2.0, formula2["Mg"])
+        for element in ["Na", "Fe", "Mg", "Al", "Si", "O", "Cr"]:
+            self.assertFloatEqual(formula[element], formula2[element] / 2.0)
 
         c1 = RelaxedComposite(c1_unrelaxed, c1_unrelaxed.reaction_basis)
         c2 = RelaxedComposite(c2_unrelaxed, c2_unrelaxed.reaction_basis)
 
+        for element in ["Na", "Fe", "Mg", "Al", "Si", "O", "Cr"]:
+            self.assertFloatEqual(c1.formula[element], c2.formula[element] / 2.0)
+            self.assertFloatEqual(c1.unrelaxed.formula[element], c1.formula[element])
+            self.assertFloatEqual(c2.unrelaxed.formula[element], c2.formula[element])
+
         c1.set_state(40.0e9, 2000.0)
         c2.set_state(40.0e9, 2000.0)
 
+        self.assertArraysAlmostEqual(c1.molar_fractions, c2.molar_fractions)
+
+        for element in ["Na", "Fe", "Mg", "Al", "Si", "O", "Cr"]:
+            self.assertFloatEqual(c1.formula[element], c2.formula[element] / 2.0)
+            self.assertFloatEqual(c1.unrelaxed.formula[element], c1.formula[element])
+            self.assertFloatEqual(c2.unrelaxed.formula[element], c2.formula[element])
+
         self.assertFloatEqual(c1.number_of_moles, c2.number_of_moles / 2.0)
-        self.assertFloatEqual(c1.formula["Mg"], c2.formula["Mg"] / 2.0)
-        self.assertFloatEqual(c1.formula["Si"], c2.formula["Si"] / 2.0)
-        self.assertFloatEqual(c1.formula["O"], c2.formula["O"] / 2.0)
         self.assertFloatEqual(c1.molar_mass, c2.molar_mass)
         self.assertFloatEqual(c1.mass, c2.mass / 2.0)
 
@@ -699,52 +741,33 @@ def test_property_extensivity_relaxed_composite(self):
         c1.set_state(40.0e9, 2000.0)
         c2.set_state(40.0e9, 2000.0)
 
-        properties_extensive = [
-            "internal_energy",
-            "helmholtz",
-            "gibbs",
-            "enthalpy",
-            "entropy",
-            "volume",
-            "heat_capacity_p",
-            "heat_capacity_v",
-            "mass",
-            "H",
-            "S",
-            "V",
-            "C_p",
-            "C_v",
-        ]
-
         for prop in properties_extensive:
             val1 = getattr(c1, prop)
             val2 = getattr(c2, prop)
             self.assertFloatEqual(val1, val2 / 2.0)
 
-        properties_intensive = [
-            "molar_internal_energy",
-            "molar_helmholtz",
-            "molar_gibbs",
-            "molar_enthalpy",
-            "molar_entropy",
-            "molar_volume",
-            "molar_heat_capacity_p",
-            "molar_heat_capacity_v",
-            "isothermal_bulk_modulus_reuss",
-            "isentropic_bulk_modulus_reuss",
-            "isothermal_compressibility_reuss",
-            "isentropic_compressibility_reuss",
-            "grueneisen_parameter",
-            "thermal_expansivity",
-            "molar_mass",
-            "density",
-        ]
-
         for prop in properties_intensive:
             val1 = getattr(c1, prop)
             val2 = getattr(c2, prop)
             self.assertFloatEqual(val1, val2)
 
+    def test_equilibrate_relaxed_composite(self):
+
+        assemblage = setup_assemblage()
+        relaxed_assemblage = RelaxedComposite(assemblage, assemblage.reaction_basis)
+        equality_constraints = [["P", 0.5e9], ["T", 900.0]]
+        equilibrate(
+            relaxed_assemblage.formula, relaxed_assemblage, equality_constraints
+        )
+        f1 = relaxed_assemblage.molar_fractions.copy()
+
+        assemblage = setup_assemblage()
+        relaxed_assemblage = RelaxedComposite(assemblage, assemblage.reaction_basis)
+        relaxed_assemblage.set_state(0.5e9, 900.0)
+        f2 = relaxed_assemblage.molar_fractions.copy()
+
+        self.assertArraysAlmostEqual(f1, f2)
+
 
 if __name__ == "__main__":
     unittest.main()