add relaxed anisotropic solution

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)

contrib/anisotropic_relaxation/data/Andrault_2003_post_stishovite_unit_cell.dat

@@ -0,0 +1,14 @@
+# P	P_err	a	a_err	b	b_err	c	c_err	V	V_err
+66.31	0.19893	3.8971	0.0007	4.0008	0.0008	2.5658	0.0002	40.005	0.013
+69.58	0.20874	3.8829	0.0007	3.9955	0.0007	2.5608	0.0003	39.729	0.013
+76.42	0.22926	3.8467	0.0011	4.0114	0.0012	2.5509	0.0004	39.363	0.021
+81.72	0.24516	3.8325	0.0009	3.9848	0.001	2.5438	0.0002	38.848	0.016
+89.39	0.26817	3.817	0.001	3.9787	0.0011	2.5339	0.0002	38.481	0.017
+94.64	0.28392	3.8014	0.001	3.975	0.0011	2.5279	0.0002	38.198	0.017
+99.41	0.29823	3.7882	0.0009	3.9709	0.001	2.5221	0.0002	37.939	0.016
+103.7	0.3111	3.7773	0.001	3.9664	0.0011	2.5169	0.0003	37.709	0.017
+106.9	0.3207	3.7675	0.0006	3.9576	0.0007	2.513	0.0002	37.469	0.011
+114.6	0.3438	3.7505	0.0007	3.9533	0.0008	2.5065	0.0003	37.164	0.013
+120.4	0.3612	3.7364	0.0007	3.9482	0.0008	2.4998	0.0002	36.876	0.013
+124.1	0.3723	3.727	0.0008	3.9449	0.0009	2.496	0.0003	36.699	0.013
+128	0.384	3.7201	0.0008	3.9422	0.0009	2.4913	0.0003	36.535	0.014

contrib/anisotropic_relaxation/data/Andrault_2003_stishovite_unit_cell.dat

@@ -0,0 +1,28 @@
+# P	P_err	a	a_err	c	c_err	V	V_err
+0.0001	0.00001	4.17755	0.00016	2.66518	0.00034	46.5126	0.0061
+1.168	0.006	4.17119	0.00027	2.66356	0.0003	46.3429	0.0053
+2.299	0.008	4.16504	0.00011	2.6618	0.00024	46.1756	0.0043
+3.137	0.008	4.16051	0.00014	2.66062	0.00029	46.055	0.0051
+4.252	0.008	4.15488	0.00013	2.65868	0.00026	45.8969	0.0045
+5.037	0.011	4.15084	0.00015	2.65767	0.0003	45.7902	0.0053
+5.851	0.009	4.14669	0.00012	2.65612	0.00022	45.6721	0.0038
+6.613	0.008	4.14323	0.00014	2.6547	0.00027	45.5715	0.005
+7.504	0.009	4.13888	0.00012	2.6534	0.00022	45.4536	0.0041
+8.264	0.011	4.13555	0.00015	2.65226	0.0003	45.3609	0.0056
+9.635	0.012	4.12962	0.00011	2.64976	0.00023	45.1885	0.0042
+11.69	0.03507	4.1217	0.0001	2.6458	0.0001	44.947	0.002
+17.67	0.05301	4.0962	0.0001	2.6381	0.0001	44.264	0.002
+22.38	0.06714	4.0781	0.0001	2.6321	0.0002	43.776	0.003
+29.38	0.08814	4.0552	0.0003	2.6192	0.0005	43.073	0.009
+37.71	0.11313	4.0287	0.0003	2.6048	0.0004	42.278	0.008
+46.03	0.13809	4.0035	0.0004	2.5919	0.001	41.544	0.017
+52.73	0.15819	3.9859	0.0004	2.5806	0.0004	40.999	0.009
+26.32	0.07896	4.0576	0.0002	2.6217	0.0003	43.164	0.006
+30.98	0.09294	4.044	0.0002	2.6155	0.0003	42.772	0.005
+34.21	0.10263	4.0308	0.0001	2.6101	0.0002	42.407	0.003
+38.45	0.11535	4.0207	0.0002	2.6037	0.0002	42.093	0.004
+43.37	0.13011	4.003	0.0002	2.5966	0.0002	41.61	0.004
+47.49	0.14247	3.9907	0.0003	2.5921	0.0004	41.28	0.007
+50.47	0.15141	3.9836	0.0002	2.5837	0.0004	41	0.008
+55.91	0.16773	3.9755	0.0003	2.5767	0.0008	40.724	0.013
+57.56	0.17268	3.9725	0.0003	2.5788	0.0009	40.695	0.015