Merge pull request #637 from bobmyhill/set_state_with_volume_mineral

Speed up set_state_with_volume for P(V, T) minerals

Author: bobmyhill

burnman/classes/material.py

@@ -186,10 +186,20 @@ def set_state_with_volume(
     ):
         """
         This function acts similarly to set_state, but takes volume and
-        temperature as input to find the pressure. In order to ensure
-        self-consistency, this function does not use any pressure functions
-        from the material classes, but instead finds the pressure using the
-        brentq root-finding method.
+        temperature as input to find the pressure.
+
+        In order to ensure self-consistency, this function does not
+        use any pressure functions from the material classes,
+        but instead finds the pressure using the brentq root-finding
+        method. To provide more context, even if a mineral is being
+        evaluated with a P(V, T) equation of state, there might be a
+        property modifier G_mod(P, T) added on top - which might then
+        introduce a pressure dependent V_mod(P, T).
+        Thus, we must solve for the volume iteratively:
+        V = V_i(P(V_i, T), T) + V_mod(P, T), where P(V_i, T) is
+        solved for iteratively.
+
+        This function is overloaded by the Mineral class.
 
         :param volume: The desired molar volume of the mineral [m^3].
         :type volume: float

burnman/classes/mineral.py

@@ -13,6 +13,8 @@
 from .material import Material, material_property
 from .. import eos
 from ..utils.misc import copy_documentation
+from ..utils.math import bracket
+from scipy.optimize import brentq
 
 
 class Mineral(Material):
@@ -145,6 +147,73 @@ def set_state(self, pressure, temperature):
                 "no method set for mineral, or equation_of_state given in mineral.params"
             )
 
+    def set_state_with_volume(
+        self, volume, temperature, pressure_guesses=[0.0e9, 10.0e9]
+    ):
+        """
+        This function acts similarly to set_state, but takes volume and
+        temperature as input to find the pressure.
+
+        If the mineral is an endmember that has a pressure(T, V)
+        function and does not have property modifiers,
+        this function evaluates the pressure directly.
+        Otherwise, it finds the pressure using the brentq root-finding
+        method on the molar_volume attribute of the mineral.
+        This ensures self-consistency even if the mineral has a
+        pressure-dependent volume modifier.
+
+        :param volume: The desired molar volume of the mineral [m^3].
+        :type volume: float
+
+        :param temperature: The desired temperature of the mineral [K].
+        :type temperature: float
+
+        :param pressure_guesses: A list of floats denoting the initial
+            low and high guesses for bracketing of the pressure [Pa].
+            These guesses should preferably bound the correct pressure,
+            but do not need to do so. More importantly,
+            they should not lie outside the valid region of
+            the equation of state. Defaults to [0.e9, 10.e9].
+        :type pressure_guesses: list
+        """
+
+        def _delta_volume(pressure, volume, temperature):
+            # Try to set the state with this pressure,
+            # but if the pressure is too low most equations of state
+            # fail. In this case, treat the molar_volume as infinite
+            # and brentq will try a larger pressure.
+            try:
+                self.set_state(pressure, temperature)
+                return volume - self.molar_volume
+            except Exception:
+                return -np.inf
+
+        if (
+            type(self.method) == str  # require an endmember method
+            or not self.method.pressure  # that has a P(V, T) function
+            or self.property_modifiers  # and that doesn't have modifiers
+        ):
+            # we need to have a sign change in [a,b] to find a zero.
+            args = (volume, temperature)
+            try:
+                sol = bracket(
+                    _delta_volume, pressure_guesses[0], pressure_guesses[1], args
+                )
+            except ValueError:
+                try:  # Try again with 0 Pa lower bound on the pressure
+                    sol = bracket(_delta_volume, 0.0e9, pressure_guesses[1], args)
+                except ValueError:
+                    raise Exception(
+                        "Cannot find a pressure, perhaps the volume or starting pressures "
+                        "are outside the range of validity for the equation of state?"
+                    )
+            pressure = brentq(_delta_volume, sol[0], sol[1], args=args)
+            self.set_state(pressure, temperature)
+        else:
+            self.set_state(
+                self.method.pressure(temperature, volume, self.params), temperature
+            )
+
     """
     Properties from equations of state
     We choose the P, T properties (e.g. Gibbs(P, T) rather than Helmholtz(V, T)),

burnman/eos/slb.py

@@ -254,41 +254,47 @@ def _K_T(V, T, params):
                         "the equation of state?"
                     )
 
-        return opt.brentq(_delta_pressure, sol[0], sol[1], args=args)
+        return opt.brentq(_delta_pressure, sol[0], sol[1], args=args, xtol=1.0e-24)
 
     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 = _grueneisen_parameter_slb(
-            params["V_0"], volume, params["grueneisen_0"], params["q_0"]
-        )
-        # does not depend on pressure
-        # thermal energy at temperature T
-        E_th = debye.thermal_energy(temperature, debye_T, params["n"])
-        # thermal energy at reference temperature
-        E_th_ref = debye.thermal_energy(params["T_0"], debye_T, params["n"])
+        T_0 = params["T_0"]
+        Debye_0 = params["Debye_0"]
+        V_0 = params["V_0"]
+        n = params["n"]
+
+        a1_ii = 6.0 * params["grueneisen_0"]  # EQ 47
+        a2_iikk = (
+            -12.0 * params["grueneisen_0"]
+            + 36.0 * pow(params["grueneisen_0"], 2.0)
+            - 18.0 * params["q_0"] * params["grueneisen_0"]
+        )  # EQ 47
 
         b_iikk = 9.0 * params["K_0"]  # EQ 28
-        b_iikkmm = 27.0 * params["K_0"] * (params["Kprime_0"] - 4.0)  # EQ 29
-        f = 0.5 * (pow(params["V_0"] / volume, 2.0 / 3.0) - 1.0)  # EQ 24
-        P = (1.0 / 3.0) * (pow(1.0 + 2.0 * f, 5.0 / 2.0)) * (
-            (b_iikk * f) + (0.5 * b_iikkmm * pow(f, 2.0))
-        ) + gr * (
-            E_th - E_th_ref
-        ) / volume  # EQ 21
+        b_iikkmm = 27.0 * params["K_0"] * (params["Kprime_0"] - 4.0)  # EQ 29z
+
+        bel_0, gel = 0.0, 1.0
         if self.conductive:
-            P += el.pressure(
-                temperature,
-                volume,
-                params["T_0"],
-                params["V_0"],
-                params["bel_0"],
-                params["gel"],
-            )
-        return P
+            bel_0, gel = params["bel_0"], params["gel"]
+
+        return _delta_pressure(
+            volume,
+            0.0,
+            temperature,
+            V_0,
+            T_0,
+            Debye_0,
+            n,
+            a1_ii,
+            a2_iikk,
+            b_iikk,
+            b_iikkmm,
+            bel_0,
+            gel,
+        )
 
     def isothermal_bulk_modulus_reuss(self, pressure, temperature, volume, params):
         """

burnman/tools/eos.py

@@ -318,9 +318,6 @@ def equilibration_function(mineral):
     m.set_state(P + 0.5 * dP, T + 0.5 * dT)
     equilibration_function(m)
 
-    beta0 = -(logm(F3) - logm(F2)) / dP
-    alpha0 = (logm(F1) - logm(F0)) / dT
-
     Q = m.deformed_coordinate_frame
     beta1 = m.isothermal_compressibility_tensor
     alpha1 = m.thermal_expansivity_tensor
@@ -329,6 +326,9 @@ def equilibration_function(mineral):
     alpha1 = np.einsum("mi, nj, ij->mn", Q, Q, alpha1)
 
     if m.orthotropic:
+        beta0 = -(logm(F3) - logm(F2)) / dP
+        alpha0 = (logm(F1) - logm(F0)) / dT
+
         expr.extend(
             [f"SI = -d(lnm(F))/dP ({i}{j})" for i in range(3) for j in range(i, 3)]
         )

misc/benchmarks/ab_initio_solids.py

@@ -58,7 +58,8 @@
     pressures = np.empty_like(volumes)
     for temperature in temperatures:
         for i, volume in enumerate(volumes):
-            pressures[i] = phase.method.pressure(temperature, volume, phase.params)
+            phase.set_state_with_volume(volume, temperature)
+            pressures[i] = phase.pressure
         ax_P.plot(
             volumes * 1e6,
             pressures / 1e9,
@@ -77,8 +78,7 @@
     energies = np.empty_like(volumes)
     for temperature in temperatures:
         for i, volume in enumerate(volumes):
-            P = phase.method.pressure(temperature, volume, phase.params)
-            phase.set_state(P, temperature)
+            phase.set_state_with_volume(volume, temperature)
             energies[i] = phase.molar_internal_energy
         ax_E.plot(
             volumes * 1e6,

test.sh

@@ -24,7 +24,7 @@ $PYTHON -m pip install -q -e .[dev]
 echo ""
 
 echo "Dependency tree:"
-$PYTHON -m pip install -q pipdeptree .
+$PYTHON -m pip install -q pipdeptree
 $PYTHON -m pipdeptree -p burnman -d 1 2> /dev/null
 pycddlib_version=`pip freeze | grep "pycddlib=" | awk -F"==" '{print $2}'`
 if [ ! -z "${pycddlib_version}" ]
@@ -64,7 +64,7 @@ with open('$t') as f:
 EOF
 ) >$t.tmp 2>$t.tmp.error
 ret=$?
-cat $t.tmp.error >>$t.tmp
+grep -v "^$" $t.tmp.error >>$t.tmp #append non-empty error messages to output
 rm -f $t.tmp.error
 
 grep -v "which is a non-GUI backend" $t.tmp | grep -v "plt.show()" > tmpfile

tests/test_composite.py

@@ -316,6 +316,20 @@ def test_query_properties(self):
         self.assertFloatEqual(rock1.molar_heat_capacity_v, min1.C_v)
         self.assertFloatEqual(rock1.molar_heat_capacity_p, min1.C_p)
 
+    def test_set_state_with_volume(self):
+        min1 = minerals.SLB_2011.periclase()
+        rock1 = burnman.Composite([min1, min1], [0.5, 0.5])
+        P = 1.0e9
+        P2 = 2.0e9
+        T = 1000.0
+        rock1.set_state(P, T)
+        V = rock1.molar_volume
+        rock1.set_state(P2, T)  # reset state
+        _ = rock1.molar_volume
+        rock1.set_state_with_volume(V, T)  # try to find original state
+
+        self.assertFloatEqual(P, rock1.pressure)
+
     def test_stoichiometric_matrix_single_mineral(self):
         rock = burnman.Composite([minerals.SLB_2011.quartz()])