diff --git a/burnman/__init__.py b/burnman/__init__.py index 54b606d0..e4d41f1f 100644 --- a/burnman/__init__.py +++ b/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. diff --git a/burnman/classes/anisotropicmineral.py b/burnman/classes/anisotropicmineral.py index d9484b02..99bfa947 100644 --- a/burnman/classes/anisotropicmineral.py +++ b/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): """ diff --git a/burnman/classes/anisotropicsolution.py b/burnman/classes/anisotropicsolution.py index 9470b651..c04c026d 100644 --- a/burnman/classes/anisotropicsolution.py +++ b/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) diff --git a/tests/test_anisotropicsolution.py b/tests/test_anisotropicsolution.py index d0800a56..922ebbf4 100644 --- a/tests/test_anisotropicsolution.py +++ b/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