Merge pull request #570 from bobmyhill/anisotropicsolution_new

New: AnisotropicSolution class

burnman/__init__.py

@@ -237,6 +237,7 @@
 from .classes.anisotropicmineral import AnisotropicMineral
 from .classes.anisotropicmineral import cell_parameters_to_vectors
 from .classes.anisotropicmineral import cell_vectors_to_parameters
+from .classes.anisotropicsolution import AnisotropicSolution
 from .classes.mineral_helpers import HelperLowHighPressureRockTransition
 from .classes.mineral_helpers import HelperSpinTransition
 from .classes.mineral_helpers import HelperRockSwitcher

burnman/classes/anisotropicsolution.py

@@ -0,0 +1,162 @@
+# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit
+# for the Earth and Planetary Sciences
+# Copyright (C) 2012 - 2024 by the BurnMan team, released under the GNU
+# GPL v2 or later.
+import numpy as np
+import copy
+from scipy.linalg import expm, logm
+from .solution import Solution
+from .anisotropicmineral import (
+    AnisotropicMineral,
+    convert_f_Pth_to_f_T_derivatives,
+    deformation_gradient_alpha_and_compliance,
+)
+from .material import material_property, cached_property, Material
+from ..utils.anisotropy import (
+    contract_stresses,
+    expand_stresses,
+    voigt_notation_to_compliance_tensor,
+    voigt_notation_to_stiffness_tensor,
+    contract_compliances,
+    contract_stiffnesses,
+)
+
+
+class AnisotropicSolution(Solution, AnisotropicMineral):
+    """
+    A class implementing the anisotropic solution model described
+    in :cite:`Myhill2024`.
+    This class is derived from Solution and AnisotropicMineral,
+    and inherits most of the methods from those classes.
+
+    Instantiation of an AnisotropicSolution is similar to that of a scalar
+    Solution, except that each of the endmembers must be an instance of
+    an AnisotropicMineral, and an additional function and parameter dictionary
+    are passed to the constructor of AnisotropicSolution that describe
+    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.
+    The output variables Psi_xs_Voigt, dPsidf_Pth_Voigt_xs and
+    dPsidPth_f_Voigt_xs (all 6x6 matrices) 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
+    set_composition().
+
+    This class is available as ``burnman.AnisotropicSolution``.
+    """
+
+    def __init__(
+        self,
+        name,
+        solution_model,
+        psi_excess_function,
+        anisotropic_parameters,
+        molar_fractions=None,
+    ):
+        # Always set orthotropic as false to ensure correct
+        # derivatives are taken.
+        # To do: relax this to allow faster calculations
+        self.orthotropic = False
+
+        self.T_0 = 298.15
+
+        # Initialise as both Material and Solution object
+        Material.__init__(self)
+        Solution.__init__(self, name, solution_model)
+
+        # Store a scalar copy of the solution model to speed up set_state
+        scalar_solution_model = copy.copy(solution_model)
+        scalar_solution_model.endmembers = [
+            [mbr.isotropic_mineral, formula] for mbr, formula in self.endmembers
+        ]
+        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.anisotropic_params = anisotropic_parameters
+        self.psi_excess_function = psi_excess_function
+
+        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)
+        V = ss.molar_volume
+        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]
+        mbrs = self.endmembers
+        for mbr in mbrs:
+            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])
+        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)
+        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)
+
+        # Now calculate the thermal pressure terms by
+        # evaluating the scalar solution at V, T_0
+        ss.set_state_with_volume(V, self.T_0)
+        P_at_T0 = ss.pressure
+        KT_at_T0 = ss.isothermal_bulk_modulus_reuss
+        self.dPthdf = KT_at_T0 - KT_at_T
+        self.Pth = pressure - P_at_T0
+
+        # And we're done with calculating endmember properties at V
+        # Now we can set state of this object and the scalar one
+        ss.set_state(pressure, temperature)
+        Solution.set_state(self, pressure, temperature)
+
+        # Calculate excess properties at the given f and Pth
+        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_full = voigt_notation_to_compliance_tensor(Psi_xs_Voigt)
+        PsiI_xs = np.einsum("ijkl, kl", Psi_xs_full, np.eye(3))
+
+        # Change of variables: (f, Pth) -> Psi(f, T)
+        aK_T = ss.alpha * ss.isothermal_bulk_modulus_reuss
+        out = convert_f_Pth_to_f_T_derivatives(
+            dPsidf_Pth_Voigt_xs, dPsidPth_f_Voigt_xs, aK_T, self.dPthdf
+        )
+        dPsidf_T_Voigt_xs, dPsiIdf_T_xs, dPsiIdT_f_xs = out
+
+        out = deformation_gradient_alpha_and_compliance(
+            self.alpha,
+            self.isothermal_compressibility_reuss,
+            PsiI_ideal + PsiI_xs,
+            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
+
+    def set_composition(self, molar_fractions):
+        self._scalar_solution.set_composition(molar_fractions)
+        Solution.set_composition(self, molar_fractions)
+        self.cell_vectors_0 = expm(
+            np.einsum("p, pij->ij", molar_fractions, self._logm_M_T_0_mbr)
+        )
+        # Calculate all other required properties
+        if self.pressure is not None:
+            self.set_state(self.pressure, self.temperature)

tests/test_anisotropicsolution.py

@@ -0,0 +1,135 @@
+from __future__ import absolute_import
+import unittest
+from util import BurnManTest
+import numpy as np
+
+from burnman.constants import Avogadro
+from burnman import AnisotropicMineral, AnisotropicSolution
+from burnman.classes.solutionmodel import SymmetricRegularSolution
+from burnman.tools.eos import check_eos_consistency
+from burnman.tools.eos import check_anisotropic_eos_consistency
+from burnman.minerals.SLB_2011 import forsterite
+from burnman.utils.unitcell import cell_parameters_to_vectors
+
+
+def make_nonorthotropic_mineral(a, b, c, alpha, beta, gamma, d, e, f):
+    fo = forsterite()
+    cell_lengths = np.array([a, b, c])
+    cell_parameters = np.array(
+        [cell_lengths[0], cell_lengths[1], cell_lengths[2], alpha, beta, gamma]
+    )
+    fo.params["V_0"] = np.linalg.det(cell_parameters_to_vectors(cell_parameters))
+    constants = np.zeros((6, 6, 3, 1))
+    constants[:, :, 1, 0] = np.array(
+        [
+            [0.44, -0.12, -0.1, d, e, f],
+            [-0.12, 0.78, -0.22, 0.0, 0.0, 0.0],
+            [-0.1, -0.22, 0.66, 0.0, 0.0, 0.0],
+            [d, 0.0, 0.0, 1.97, 0.0, 0.0],
+            [e, 0.0, 0.0, 0.0, 1.61, 0.0],
+            [f, 0.0, 0.0, 0.0, 0.0, 1.55],
+        ]
+    )
+    constants[:, :, 2, 0] = np.array(
+        [
+            [0.24, -0.12, -0.1, 0.0, 0.0, 0.0],
+            [-0.12, 0.38, -0.22, d, d, d],
+            [-0.1, -0.22, 0.26, f, e, d],
+            [0.0, d, f, 0.0, 0.0, 0.0],
+            [0.0, d, e, 0.0, 0.0, 0.0],
+            [0.0, d, d, 0.0, 0.0, 0.0],
+        ]
+    )
+
+    m = AnisotropicMineral(fo, cell_parameters, constants)
+    return m
+
+
+def make_nonorthotropic_solution():
+
+    cell_lengths_A = np.array([4.7646, 10.2296, 5.9942])
+    lth = cell_lengths_A * 1.0e-10 * np.cbrt(Avogadro / 4.0)
+
+    m1 = make_nonorthotropic_mineral(
+        lth[0], lth[1], lth[2], 85.0, 80.0, 87.0, 0.4, -1.0, -0.6
+    )
+    m2 = make_nonorthotropic_mineral(
+        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]],
+    )
+
+    def fn(lnV, Pth, X, params):
+        z = np.zeros((6, 6))
+        return (z, z, z)
+
+    prm = {}
+
+    ss = AnisotropicSolution(
+        name="double forsterite",
+        solution_model=solution_model,
+        psi_excess_function=fn,
+        anisotropic_parameters=prm,
+    )
+
+    return ss
+
+
+class test_two_member_solution(BurnManTest):
+    def test_volume(self):
+        ss = make_nonorthotropic_solution()
+        ps = [0.2, 0.8]
+        ss.set_composition(ps)
+        Vs = np.array([np.linalg.det(mbr[0].cell_vectors_0) for mbr in ss.endmembers])
+        V1 = np.power(Vs[0], ps[0]) * np.power(Vs[1], ps[1])
+        V2 = np.linalg.det(ss.cell_vectors_0)
+        self.assertFloatEqual(V1, V2)
+
+    def test_non_orthotropic_endmember_consistency(self):
+        ss = make_nonorthotropic_solution()
+        ss.set_composition([1.0, 0.0])
+        self.assertTrue(check_anisotropic_eos_consistency(ss, P=2.0e10, tol=1.0e-6))
+
+    def test_non_orthotropic_solution_consistency(self):
+        ss = make_nonorthotropic_solution()
+        ss.set_composition([0.8, 0.2])
+        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
+        ss1.set_composition([0.8, 0.2])
+        self.assertTrue(check_eos_consistency(ss1, P=2.0e10, T=1000.0, tol=1.0e-6))
+
+        ss.set_composition([0.8, 0.2])
+        ss.set_state(2.0e10, 1000.0)
+        self.assertFloatEqual(
+            ss.isothermal_bulk_modulus_reuss, ss1.isothermal_bulk_modulus_reuss
+        )
+
+    def test_stiffness(self):
+        ss = make_nonorthotropic_solution()
+        ss.set_composition([0.8, 0.2])
+        ss.set_state(1.0e9, 300.0)
+        Cijkl = ss.full_isothermal_stiffness_tensor
+        Cij = ss.isothermal_stiffness_tensor
+
+        self.assertFloatEqual(Cij[0, 0], Cijkl[0, 0, 0, 0])
+        self.assertFloatEqual(Cij[1, 1], Cijkl[1, 1, 1, 1])
+        self.assertFloatEqual(Cij[2, 2], Cijkl[2, 2, 2, 2])
+        self.assertFloatEqual(Cij[0, 1], Cijkl[0, 0, 1, 1])
+        self.assertFloatEqual(Cij[0, 2], Cijkl[0, 0, 2, 2])
+        self.assertFloatEqual(Cij[1, 2], Cijkl[1, 1, 2, 2])
+        self.assertFloatEqual(Cij[3, 3], Cijkl[1, 2, 1, 2])
+        self.assertFloatEqual(Cij[4, 4], Cijkl[0, 2, 0, 2])
+        self.assertFloatEqual(Cij[5, 5], Cijkl[0, 1, 0, 1])
+
+
+if __name__ == "__main__":
+    unittest.main()