diff --git a/burnman/__init__.py b/burnman/__init__.py index afda8f73..f8e00f3b 100644 --- a/burnman/__init__.py +++ b/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 diff --git a/burnman/classes/anisotropicsolution.py b/burnman/classes/anisotropicsolution.py new file mode 100644 index 00000000..6d8ca465 --- /dev/null +++ b/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) diff --git a/tests/test_anisotropicsolution.py b/tests/test_anisotropicsolution.py new file mode 100644 index 00000000..d0800a56 --- /dev/null +++ b/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()