refactored anisotropic functions

burnman/classes/anisotropicmineral.py

@@ -18,6 +18,119 @@
 )
 
 
+def convert_f_Pth_to_f_T_derivatives(
+    dPsidf_Pth_Voigt, dPsidPth_f_Voigt, aK_T, dPthdf_T
+):
+    """
+    Convenience function converting Psi and its derivatives with respect to
+    f and Pth into derivatives of f and T using the chain rule.
+
+    :param dPsidf_Pth_Voigt: The first derivative of the anisotropic tensor
+        Psi with respect to f at constant Pth (Voigt form).
+    :type dPsiIdf_Pth: numpy array (6x6)
+    :param dPsidPth_f_Voigt: The first derivative of the anisotropic tensor
+        Psi with respect to Pth at constant f (Voigt form).
+    :type dPsiIdPth_f: numpy array (6x6)
+    :param aK_T: The volumetric thermal expansivity multiplied by
+        the Reuss isothermal bulk modulus.
+    :type aK_T: float
+    :param dPthdf: The change in thermal pressure with respect to volume
+        at constant temperature.
+    :type dPthdf: float
+
+    :returns: The first derivative of Psi with respect to f in Voigt form,
+        and the first derivatives of PsiI with respect to f and T.
+    :rtype: Tuple of three objects of type numpy.array (2D)
+    """
+    # The following manipulation uses the chain rule to convert
+    # from Psi(f, Pth) and derivatives to Psi(f, T) and derivatives
+    dPsidf_full = voigt_notation_to_compliance_tensor(dPsidf_Pth_Voigt)
+    dPsidPth_full = voigt_notation_to_compliance_tensor(dPsidPth_f_Voigt)
+
+    dPsiIdf_Pth = np.einsum("ijkl, kl", dPsidf_full, np.eye(3))
+    dPsiIdPth_f = np.einsum("ijkl, kl", dPsidPth_full, np.eye(3))
+
+    dPsidf_T_Voigt = dPsidf_Pth_Voigt + dPsidPth_f_Voigt * dPthdf_T
+    dPsiIdf_T = dPsiIdf_Pth + dPsiIdPth_f * dPthdf_T
+    dPsiIdT_f = aK_T * dPsiIdPth_f
+    return (dPsidf_T_Voigt, dPsiIdf_T, dPsiIdT_f)
+
+
+def deformation_gradient_alpha_and_compliance(
+    alpha_V, beta_TR, PsiI, dPsidf_T_Voigt, dPsiIdf_T, dPsiIdT_f
+):
+    """
+    Convenience function converting Psi and its derivatives with respect to
+    f and T into the deformation gradient, thermal expansivity and
+    isothermal compliance tensors. This function is appropriate for
+    both orthotropic and non-orthotropic materials.
+
+    :param alpha_V: The volumetric thermal expansivity.
+    :type alpha_V: float
+    :param beta_TR: The Reuss isothermal compressibility.
+    :type beta_TR: float
+    :param psiI: The anisotropic tensor Psi matrix multiplied with
+        the identity matrix.
+    :type psiI: numpy array (3x3)
+    :param dPsidf_T_Voigt: Voigt-form matrix of the first
+        derivative of the anisotropic tensor Psi with respect to f
+        at constant T.
+    :type dPsidf_T_Voigt: numpy array (6x6)
+    :param dPsiIdf_T: The first derivative of PsiI
+        with respect to f at constant T.
+    :type dPsiIdf_T: numpy array (3x3)
+    :param dPsiIdT_f: The first derivative of PsiI
+        with respect to T at constant f.
+    :type dPsiIdT_f: numpy array (3x3)
+
+    :returns: The unrotated isothermal compliance tensor in Voigt form (6x6),
+        and the thermal expansivity tensor (3x3).
+    :rtype: Tuple of two objects of type numpy.array (2D)
+    """
+    # Numerical derivatives with respect to f and T
+    df = 1.0e-7
+    F0 = expm(PsiI - dPsiIdf_T * df / 2.0)
+    F1 = expm(PsiI + dPsiIdf_T * df / 2.0)
+    dFdf_T = (F1 - F0) / df
+
+    dT = 0.1
+    F0 = expm(PsiI - dPsiIdT_f * dT / 2.0)
+    F1 = expm(PsiI + dPsiIdT_f * dT / 2.0)
+    dFdT_f = (F1 - F0) / dT
+
+    # Convert to pressure and temperature derivatives
+    dFdP_T = -beta_TR * dFdf_T
+    dFdT_P = dFdT_f + alpha_V * dFdf_T
+
+    # Calculate the unrotated isothermal compressibility
+    # and unrotated thermal expansivity tensors
+    F = expm(PsiI)
+    invF = np.linalg.inv(F)
+    LP = np.einsum("ij,kj->ik", dFdP_T, invF)
+    beta_T = -0.5 * (LP + LP.T)
+    LT = np.einsum("ij,kj->ik", dFdT_P, invF)
+    alpha = 0.5 * (LT + LT.T)
+
+    # Calculate the unrotated isothermal compliance
+    # tensor in Voigt form.
+    S_T = beta_TR * dPsidf_T_Voigt
+    for i, j, k in [[0, 1, 2], [1, 2, 0], [0, 2, 1]]:
+        S_T[i][j] = 0.5 * (
+            -S_T[i][i]
+            - S_T[j][j]
+            + S_T[k][k]
+            + beta_T[i][i]
+            + beta_T[j][j]
+            - beta_T[k][k]
+        )
+        S_T[j][i] = S_T[i][j]
+
+        S_T[i][i + 3] = 2.0 * beta_T[j][k] - S_T[j][i + 3] - S_T[k][i + 3]
+        S_T[i + 3][i] = S_T[i][i + 3]
+
+    return F, alpha, S_T
+
+
 class AnisotropicMineral(Mineral, AnisotropicMaterial):
     """
     A class implementing the anisotropic mineral equation of state described
@@ -103,7 +216,7 @@ def __init__(
             )
             raise Exception(
                 "The standard state unit vectors are inconsistent "
-                "with the volume. Suggest multiplying each vector length"
+                "with the volume. Suggest multiplying each vector length "
                 f"by {factor}."
             )
 
@@ -161,88 +274,34 @@ def set_state(self, pressure, temperature):
         f = np.log(Vrel)
 
         out = self.psi_function(f, self.Pth, self.anisotropic_params)
-        Psi_Voigt, dPsidf_Voigt, dPsidPth_Voigt = out
+        Psi_Voigt, dPsidf_Pth_Voigt, dPsidPth_f_Voigt = out
         Psi_full = voigt_notation_to_compliance_tensor(Psi_Voigt)
-        dPsidf_full = voigt_notation_to_compliance_tensor(dPsidf_Voigt)
-        dPsidPth_full = voigt_notation_to_compliance_tensor(dPsidPth_Voigt)
-
-        PsiI = np.einsum("ijkl, kl", Psi_full, np.eye(3))
-        dPsidfI = np.einsum("ijkl, kl", dPsidf_full, np.eye(3))
-        dPsidPthI = np.einsum("ijkl, kl", dPsidPth_full, np.eye(3))
+        self._PsiI = np.einsum("ijkl, kl", Psi_full, np.eye(3))
 
-        # dPsidP_Voigt is needed for both orthotropic and nonorthotropic
-        # materials
-        dPsidP_Voigt = -self.isothermal_compressibility_reuss * (
-            dPsidf_Voigt + dPsidPth_Voigt * self.dPthdf
+        # Change of variables: (f, Pth) -> Psi(f, T)
+        aK_T = self.alpha * self.isothermal_bulk_modulus_reuss
+        out = convert_f_Pth_to_f_T_derivatives(
+            dPsidf_Pth_Voigt, dPsidPth_f_Voigt, aK_T, self.dPthdf
         )
+        self._dPsidf_T_Voigt, self._dPsiIdf_T, self._dPsiIdT_f = out
 
-        # Calculate F, dFdP, dFdT
-        self._F = expm(PsiI)
-
+        # Calculate F, thermal expansivity and compliance, dFdT
         if self.orthotropic:
-            self._S_T_unrotated_Voigt = -dPsidP_Voigt
-
-            self._unrotated_alpha = self.alpha * (
-                dPsidfI + dPsidPthI * (self.dPthdf + self.isothermal_bulk_modulus_reuss)
+            self._unrotated_F = expm(self._PsiI)
+            self._unrotated_alpha = self._dPsiIdT_f + self.alpha * self._dPsiIdf_T
+            self._unrotated_S_T_Voigt = (
+                self.isothermal_compressibility_reuss * self._dPsidf_T_Voigt
             )
         else:
-            # Numerical derivatives with respect to f and Pth
-            df = f * 1.0e-5
-            Psi_full0 = voigt_notation_to_compliance_tensor(
-                Psi_Voigt - dPsidf_Voigt * df / 2.0
+            out = deformation_gradient_alpha_and_compliance(
+                self.alpha,
+                self.isothermal_compressibility_reuss,
+                self._PsiI,
+                self._dPsidf_T_Voigt,
+                self._dPsiIdf_T,
+                self._dPsiIdT_f,
             )
-            Psi_full1 = voigt_notation_to_compliance_tensor(
-                Psi_Voigt + dPsidf_Voigt * df / 2.0
-            )
-            F0 = expm(np.einsum("ijkl, kl", Psi_full0, np.eye(3)))
-            F1 = expm(np.einsum("ijkl, kl", Psi_full1, np.eye(3)))
-            dFdf = (F1 - F0) / df
-
-            dPth = self.Pth * 1.0e-5
-            Psi_full0 = voigt_notation_to_compliance_tensor(
-                Psi_Voigt - dPsidPth_Voigt * dPth / 2.0
-            )
-            Psi_full1 = voigt_notation_to_compliance_tensor(
-                Psi_Voigt + dPsidPth_Voigt * dPth / 2.0
-            )
-            F0 = expm(np.einsum("ijkl, kl", Psi_full0, np.eye(3)))
-            F1 = expm(np.einsum("ijkl, kl", Psi_full1, np.eye(3)))
-            dFdPth = (F1 - F0) / dPth
-
-            # Convert to pressure and temperature derivatives
-            dFdP = -self.isothermal_compressibility_reuss * (
-                dFdf + dFdPth * self.dPthdf
-            )
-            dFdT = self.alpha * (
-                dFdf + dFdPth * (self.dPthdf + self.isothermal_bulk_modulus_reuss)
-            )
-
-            # Calculate the unrotated isothermal compressibility
-            # and unrotated thermal expansivity tensors
-            invF = np.linalg.inv(self._F)
-            LP = np.einsum("ij,kj->ik", dFdP, invF)
-            beta_T = -0.5 * (LP + LP.T)
-            LT = np.einsum("ij,kj->ik", dFdT, invF)
-            self._unrotated_alpha = 0.5 * (LT + LT.T)
-
-            # Calculate the unrotated isothermal compliance
-            # tensor in Voigt form.
-            S = -dPsidP_Voigt
-            for i, j, k in [[0, 1, 2], [1, 2, 0], [0, 2, 1]]:
-                S[i][j] = 0.5 * (
-                    -S[i][i]
-                    - S[j][j]
-                    + S[k][k]
-                    + beta_T[i][i]
-                    + beta_T[j][j]
-                    - beta_T[k][k]
-                )
-                S[j][i] = S[i][j]
-
-                S[i][i + 3] = 2.0 * beta_T[j][k] - S[j][i + 3] - S[k][i + 3]
-                S[i + 3][i] = S[i][i + 3]
-
-            self._S_T_unrotated_Voigt = S
+            self._unrotated_F, self._unrotated_alpha, self._unrotated_S_T_Voigt = out
 
     @material_property
     def deformation_gradient_tensor(self):
@@ -252,7 +311,7 @@ def deformation_gradient_tensor(self):
             (i.e. the state at the reference pressure and temperature).
         :rtype: numpy.array (2D)
         """
-        return self._F
+        return self._unrotated_F
 
     @material_property
     def unrotated_cell_vectors(self):
@@ -366,6 +425,8 @@ def isothermal_bulk_modulus(self):
             "isothermal_bulk_modulus_reuss?"
         )
 
+    isothermal_bulk_modulus_reuss = Mineral.isothermal_bulk_modulus
+
     @material_property
     def isentropic_bulk_modulus(self):
         """
@@ -382,7 +443,7 @@ def isentropic_bulk_modulus(self):
             "isentropic_bulk_modulus_reuss?"
         )
 
-    isothermal_bulk_modulus_reuss = Mineral.isothermal_bulk_modulus
+    isentropic_bulk_modulus_reuss = Mineral.adiabatic_bulk_modulus
 
     @material_property
     def isothermal_compressibility(self):
@@ -422,16 +483,7 @@ def isothermal_bulk_modulus_voigt(self):
         :returns: The Voigt bound on the isothermal bulk modulus in [Pa].
         :rtype: float
         """
-        K = (
-            np.sum(
-                [
-                    [self.isothermal_stiffness_tensor[i][k] for k in range(3)]
-                    for i in range(3)
-                ]
-            )
-            / 9.0
-        )
-        return K
+        return np.sum(self.isothermal_stiffness_tensor[:3, :3]) / 9.0
 
     @material_property
     def isothermal_compressibility_reuss(self):
@@ -477,10 +529,10 @@ def isothermal_compliance_tensor(self):
         :rtype: numpy.array (2D)
         """
         if self.orthotropic:
-            return self._S_T_unrotated_Voigt
+            return self._unrotated_S_T_Voigt
         else:
             R = self.rotation_matrix
-            S = voigt_notation_to_compliance_tensor(self._S_T_unrotated_Voigt)
+            S = voigt_notation_to_compliance_tensor(self._unrotated_S_T_Voigt)
             S_rotated = np.einsum("mi, nj, ok, pl, ijkl->mnop", R, R, R, R, S)
             return contract_compliances(S_rotated)
 

burnman/classes/material.py

@@ -202,7 +202,7 @@ def set_state_with_volume(
             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 [5.e9, 10.e9].
+            the equation of state. Defaults to [0.e9, 10.e9].
         :type pressure_guesses: list
         """
 
@@ -222,10 +222,13 @@ def _delta_volume(pressure, volume, temperature):
         try:
             sol = bracket(_delta_volume, pressure_guesses[0], 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?"
-            )
+            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)
 
@@ -613,6 +616,16 @@ def molar_heat_capacity_p(self):
             "need to implement molar_heat_capacity_p() in derived class!"
         )
 
+    @material_property
+    def isentropic_thermal_gradient(self):
+        """
+        :returns: dTdP, the change in temperature with pressure at constant entropy [Pa/K]
+        :rtype: float
+        """
+        raise NotImplementedError(
+            "need to implement isentropic_thermal_gradient() in derived class!"
+        )
+
     #
     # Aliased properties
     @property

burnman/classes/mineral.py

@@ -1,4 +1,5 @@
-# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for the Earth and Planetary Sciences
+# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit
+# for the Earth and Planetary Sciences
 # Copyright (C) 2012 - 2017 by the BurnMan team, released under the GNU
 # GPL v2 or later.
 
@@ -197,6 +198,8 @@ def isothermal_bulk_modulus(self):
             - self._property_modifiers["d2GdP2"]
         )
 
+    isothermal_bulk_modulus_reuss = isothermal_bulk_modulus
+
     @material_property
     @copy_documentation(Material.molar_heat_capacity_p)
     def molar_heat_capacity_p(self):
@@ -379,3 +382,10 @@ def molar_heat_capacity_v(self):
             * self.thermal_expansivity
             * self.isothermal_bulk_modulus
         )
+
+    @material_property
+    @copy_documentation(Material.isentropic_thermal_gradient)
+    def isentropic_thermal_gradient(self):
+        return (
+            self.molar_volume * self.temperature * self.thermal_expansivity
+        ) / self.molar_heat_capacity_p

burnman/tools/eos.py

@@ -191,7 +191,7 @@ def equilibration_function(mineral):
 
 
 def check_anisotropic_eos_consistency(
-    m, P=1.0e9, T=2000.0, tol=1.0e-4, verbose=False, equilibration_function=None
+    m, P=1.0e9, T=300.0, tol=1.0e-4, verbose=False, equilibration_function=None
 ):
     """
     Checks that numerical derivatives of the Gibbs energy of an anisotropic mineral