added anisotropic relaxation

burnman/__init__.py

@@ -1,6 +1,6 @@
 # This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for
 # the Earth and Planetary Sciences
-# Copyright (C) 2012 - 2022 by the BurnMan team, released under the GNU
+# Copyright (C) 2012 - 2024 by the BurnMan team, released under the GNU
 # GPL v2 or later.
 
 

burnman/classes/anisotropicmineral.py

@@ -83,9 +83,11 @@ def deformation_gradient_alpha_and_compliance(
         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)
+    :returns: The deformation gradient tensor (3x3),
+        derivative with respect to lnV (3x3),
+        the thermal expansivity tensor (3x3),
+        and the unrotated isothermal compliance tensor in Voigt form (6x6).
+    :rtype: Tuple of four objects of type numpy.array (2D)
     """
     # Numerical derivatives with respect to f and T
     df = 1.0e-7
@@ -128,7 +130,7 @@ def deformation_gradient_alpha_and_compliance(
         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
+    return F, dFdf_T, alpha, S_T
 
 
 class AnisotropicMineral(Mineral, AnisotropicMaterial):
@@ -309,7 +311,12 @@ def set_state(self, pressure, temperature):
                 self._dPsiIdf_T,
                 self._dPsiIdT_f,
             )
-            self._unrotated_F, self._unrotated_alpha, self._unrotated_S_T_Voigt = out
+            (
+                self._unrotated_F,
+                self._unrotated_dFdf_T,
+                self._unrotated_alpha,
+                self._unrotated_S_T_Voigt,
+            ) = out
 
     @material_property
     def deformation_gradient_tensor(self):
@@ -458,8 +465,9 @@ def beta_T(self):
         """
         Anisotropic minerals do not have a single isentropic compressibility.
         This function returns a NotImplementedError. Users should instead
-        consider either using isothermal_compressibility_reuss,
-        isothermal_compressibility_voigt (both derived from AnisotropicMineral),
+        consider using isothermal_compressibility_tensor,
+        isothermal_compressibility_reuss or isothermal_compressibility_voigt
+        (all derived from AnisotropicMineral),
         or directly querying the elements in the isothermal_compliance_tensor.
         """
         raise NotImplementedError(
@@ -474,8 +482,9 @@ def beta_S(self):
         """
         Anisotropic minerals do not have a single isentropic compressibility.
         This function returns a NotImplementedError. Users should instead
-        consider either using isentropic_compressibility_reuss,
-        isentropic_compressibility_voigt (both derived from AnisotropicMineral),
+        consider either using isentropic_compressibility_tensor,
+        isentropic_compressibility_reuss or isentropic_compressibility_voigt
+        (all derived from AnisotropicMineral),
         or directly querying the elements in the isentropic_compliance_tensor.
         """
         raise NotImplementedError(
@@ -511,6 +520,14 @@ def isothermal_compressibility_voigt(self):
         """
         return 1.0 / self.isothermal_bulk_modulus_voigt
 
+    @material_property
+    def isentropic_bulk_modulus_voigt(self):
+        """
+        :returns: The Voigt bound on the isentropic bulk modulus in [Pa].
+        :rtype: float
+        """
+        return np.sum(self.isentropic_stiffness_tensor[:3, :3]) / 9.0
+
     @material_property
     def isentropic_compressibility_reuss(self):
         """

burnman/classes/anisotropicsolution.py

@@ -24,7 +24,7 @@
 class AnisotropicSolution(Solution, AnisotropicMineral):
     """
     A class implementing the anisotropic solution model described
-    in :cite:`Myhill2024`.
+    in :cite:`Myhill2024a`.
     This class is derived from Solution and AnisotropicMineral,
     and inherits most of the methods from those classes.
 
@@ -35,9 +35,11 @@ class AnisotropicSolution(Solution, AnisotropicMineral):
     excess contributions to the anisotropic state tensor (Psi_xs)
     and its derivatives with respect to volume and temperature.
     The function arguments should be ln(V), Pth,
-    X (a vector) and the parameter dictionary, in that order.
+    p (a vector of proportions) and the parameter dictionary,
+    in that order.
     The output variables Psi_xs_Voigt, dPsidf_Pth_Voigt_xs and
-    dPsidPth_f_Voigt_xs (all 6x6 matrices) must be returned in that
+    dPsidPth_f_Voigt_xs (all 6x6 matrices) and dPsiIdp_xs
+    (a 3 x 3 x n_endmember matrix) must be returned in that
     order in a tuple.
     States of the mineral can only be queried after setting the
     pressure and temperature using set_state() and the composition using
@@ -63,7 +65,7 @@ def __init__(
 
         # Initialise as both Material and Solution object
         Material.__init__(self)
-        Solution.__init__(self, name, solution_model)
+        Solution.__init__(self, name, solution_model, molar_fractions)
 
         # Store a scalar copy of the solution model to speed up set_state
         scalar_solution_model = copy.copy(solution_model)
@@ -72,21 +74,20 @@ def __init__(
         ]
         self._scalar_solution = Solution(name, scalar_solution_model, molar_fractions)
 
-        self._logm_M_T_0_mbr = np.array(
-            [logm(m[0].cell_vectors_0) for m in self.endmembers]
+        self._logm_M0_mbr = np.einsum(
+            "kij->ijk", np.array([logm(m[0].cell_vectors_0.T) for m in self.endmembers])
         )
+
         self.anisotropic_params = anisotropic_parameters
         self.psi_excess_function = psi_excess_function
 
+        # Finally, set the composition
         if molar_fractions is not None:
             self.set_composition(molar_fractions)
 
     def set_state(self, pressure, temperature):
         # Set solution conditions
         ss = self._scalar_solution
-        if not hasattr(ss, "molar_fractions"):
-            raise Exception("To use this EoS, you must first set the composition")
-
         ss.set_state(pressure, temperature)
 
         try:
@@ -105,7 +106,6 @@ def compute_base_properties(self):
         KT_at_T = ss.isothermal_bulk_modulus_reuss
         f = np.log(V)
         self._f = f
-
         # Evaluate endmember properties at V, T
         # Here we explicitly manipulate each of the anisotropic endmembers
         pressure_guesses = [np.max([0.0e9, pressure - 2.0e9]), pressure + 2.0e9]
@@ -114,13 +114,13 @@ def compute_base_properties(self):
             mbr[0].set_state_with_volume(V, temperature, pressure_guesses)
 
         # Endmember cell vectors and other functions of Psi (all unrotated)
-        PsiI_mbr = np.array([m[0]._PsiI for m in mbrs])
+        self._PsiI_mbr = np.einsum("kij->ijk", np.array([m[0]._PsiI for m in mbrs]))
         dPsidf_T_Voigt_mbr = np.array([m[0]._dPsidf_T_Voigt for m in mbrs])
         dPsiIdf_T_mbr = np.array([m[0]._dPsiIdf_T for m in mbrs])
         dPsiIdT_f_mbr = np.array([m[0]._dPsiIdT_f for m in mbrs])
 
         fr = self.molar_fractions
-        PsiI_ideal = np.einsum("p,pij->ij", fr, PsiI_mbr)
+        PsiI_ideal = np.einsum("ijp, p->ij", self._PsiI_mbr, fr)
         dPsidf_T_Voigt_ideal = np.einsum("p,pij->ij", fr, dPsidf_T_Voigt_mbr)
         dPsiIdf_T_ideal = np.einsum("p,pij->ij", fr, dPsiIdf_T_mbr)
         dPsiIdT_f_ideal = np.einsum("p,pij->ij", fr, dPsiIdT_f_mbr)
@@ -142,7 +142,8 @@ def compute_base_properties(self):
         out = self.psi_excess_function(
             f, self.Pth, self.molar_fractions, self.anisotropic_params
         )
-        Psi_xs_Voigt, dPsidf_Pth_Voigt_xs, dPsidPth_f_Voigt_xs = out
+        Psi_xs_Voigt, dPsidf_Pth_Voigt_xs, dPsidPth_f_Voigt_xs = out[:3]
+        self._dPsiIdp_xs = out[3]
         Psi_xs_full = voigt_notation_to_compliance_tensor(Psi_xs_Voigt)
         PsiI_xs = np.einsum("ijkl, kl", Psi_xs_full, np.eye(3))
 
@@ -153,20 +154,27 @@ def compute_base_properties(self):
         )
         dPsidf_T_Voigt_xs, dPsiIdf_T_xs, dPsiIdT_f_xs = out
 
+        self._PsiI = PsiI_ideal + PsiI_xs
+
         out = deformation_gradient_alpha_and_compliance(
-            self.alpha,
-            self.isothermal_compressibility_reuss,
-            PsiI_ideal + PsiI_xs,
+            ss.alpha,
+            ss.isothermal_compressibility_reuss,
+            self._PsiI,
             dPsidf_T_Voigt_ideal + dPsidf_T_Voigt_xs,
             dPsiIdf_T_ideal + dPsiIdf_T_xs,
             dPsiIdT_f_ideal + dPsiIdT_f_xs,
         )
-        self._unrotated_F, self._unrotated_alpha, self._unrotated_S_T_Voigt = out
+        (
+            self._unrotated_F,
+            self._unrotated_dFdf_T,
+            self._unrotated_alpha,
+            self._unrotated_S_T_Voigt,
+        ) = out
 
     def set_composition(self, molar_fractions):
         self._scalar_solution.set_composition(molar_fractions)
         self.cell_vectors_0 = expm(
-            np.einsum("p, pij->ij", molar_fractions, self._logm_M_T_0_mbr)
+            np.einsum("ijk, k ->ji", self._logm_M0_mbr, molar_fractions)
         )
         # Calculate all other required properties
         Solution.set_composition(self, molar_fractions)

tests/test_anisotropicsolution.py

@@ -5,6 +5,7 @@
 
 from burnman.constants import Avogadro
 from burnman import AnisotropicMineral, AnisotropicSolution
+from burnman import RelaxedAnisotropicSolution
 from burnman.classes.solutionmodel import SymmetricRegularSolution
 from burnman.tools.eos import check_eos_consistency
 from burnman.tools.eos import check_anisotropic_eos_consistency
@@ -45,7 +46,7 @@ def make_nonorthotropic_mineral(a, b, c, alpha, beta, gamma, d, e, f):
     return m
 
 
-def make_nonorthotropic_solution():
+def make_nonorthotropic_solution(two_fos=False):
 
     cell_lengths_A = np.array([4.7646, 10.2296, 5.9942])
     lth = cell_lengths_A * 1.0e-10 * np.cbrt(Avogadro / 4.0)
@@ -57,15 +58,25 @@ def make_nonorthotropic_solution():
         lth[0] * 1.2, lth[1] * 1.4, lth[2], 90.0, 90.0, 90.0, 0.4, -1.0, -0.6
     )
 
-    solution_model = SymmetricRegularSolution(
-        endmembers=[[m1, "[Mg]2SiO4"], [m2, "[Fe]2SiO4"]],
-        energy_interaction=[[10.0e3]],
-        volume_interaction=[[1.0e-6]],
-    )
+    if two_fos:
+        n_mbrs = 3
+        solution_model = SymmetricRegularSolution(
+            endmembers=[[m1, "[Mg]2SiO4"], [m1, "[Mg]2SiO4"], [m2, "[Fe]2SiO4"]],
+            energy_interaction=[[-1.0e3, 10.0e3], [10.0e3]],
+            volume_interaction=[[0.0, 1.0e-6], [1.0e-6]],
+        )
+    else:
+        n_mbrs = 2
+        solution_model = SymmetricRegularSolution(
+            endmembers=[[m1, "[Mg]2SiO4"], [m2, "[Fe]2SiO4"]],
+            energy_interaction=[[10.0e3]],
+            volume_interaction=[[1.0e-6]],
+        )
 
     def fn(lnV, Pth, X, params):
         z = np.zeros((6, 6))
-        return (z, z, z)
+        f = np.zeros((3, 3, n_mbrs))
+        return (z, z, z, f)
 
     prm = {}
 
@@ -101,6 +112,16 @@ def test_non_orthotropic_solution_consistency(self):
             check_anisotropic_eos_consistency(ss, P=2.0e10, T=1000.0, tol=1.0e-6)
         )
 
+    def test_relaxed_non_orthotropic_solution_consistency(self):
+        ss = make_nonorthotropic_solution(two_fos=True)
+        ss = RelaxedAnisotropicSolution(
+            ss, [[1.0, -1.0, 0.0]], [[0.5, 0.5, 0.0], [0.0, 0.0, 1.0]]
+        )
+        ss.set_composition([0.8, 0.2], relaxed=False)
+        self.assertTrue(
+            check_anisotropic_eos_consistency(ss, P=2.0e10, T=1000.0, tol=1.0e-6)
+        )
+
     def test_non_orthotropic_solution_clone(self):
         ss = make_nonorthotropic_solution()
         ss1 = ss._scalar_solution