diff --git a/burnman/classes/anisotropicmineral.py b/burnman/classes/anisotropicmineral.py index 7d73217b9..872ca6b5c 100644 --- a/burnman/classes/anisotropicmineral.py +++ b/burnman/classes/anisotropicmineral.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 - 2021 by the BurnMan team, released under the GNU +# Copyright (C) 2012 - 2024 by the BurnMan team, released under the GNU # GPL v2 or later. import numpy as np from scipy.linalg import expm, logm @@ -91,6 +91,7 @@ def __init__( self.anisotropic_params = anisotropic_parameters self.psi_function = psi_function + # cell_vectors is the transpose of the cell tensor M self.cell_vectors_0 = cell_parameters_to_vectors(cell_parameters) if ( @@ -102,14 +103,10 @@ def __init__( ) raise Exception( "The standard state unit vectors are inconsistent " - "with the volume. Suggest multiplying each " + "with the volume. Suggest multiplying each vector length" f"by {factor}." ) - # Note, Psi_0 may be asymmetric, in which case the Voigt contraction - # cannot be applied - self.Psi_0 = np.einsum("ij, kl", logm(self.cell_vectors_0), np.eye(3) / 3.0) - self.isotropic_mineral = isotropic_mineral if "name" in isotropic_mineral.params: self.name = isotropic_mineral.params["name"] @@ -162,22 +159,91 @@ def set_state(self, pressure, temperature): V_0 = self.params["V_0"] Vrel = V / V_0 f = np.log(Vrel) - self._Vrel = Vrel - self._f = f out = self.psi_function(f, self.Pth, self.anisotropic_params) - Psi_Voigt, self.dPsidf_Voigt, self.dPsidPth_Voigt = out - self.Psi = voigt_notation_to_compliance_tensor(Psi_Voigt) + self.Psi_0 - - # Convert to (f, T) variables - self.dPsidP_Voigt = -self.isothermal_compressibility_reuss * ( - self.dPsidf_Voigt + self.dPsidPth_Voigt * self.dPthdf - ) - self.dPsidT_Voigt = self.alpha * ( - self.dPsidf_Voigt - + self.dPsidPth_Voigt * (self.dPthdf + self.isothermal_bulk_modulus_reuss) + Psi_Voigt, dPsidf_Voigt, dPsidPth_Voigt = out + Psi_full = voigt_notation_to_compliance_tensor(Psi_Voigt) + dPsidf_full = voigt_notation_to_compliance_tensor(dPsidf_Voigt) + dPsidPth_full = voigt_notation_to_compliance_tensor(dPsidPth_Voigt) + + PsiI = np.einsum("ijkl, kl", Psi_full, np.eye(3)) + dPsidfI = np.einsum("ijkl, kl", dPsidf_full, np.eye(3)) + dPsidPthI = np.einsum("ijkl, kl", dPsidPth_full, np.eye(3)) + + # dPsidP_Voigt is needed for both orthotropic and nonorthotropic + # materials + dPsidP_Voigt = -self.isothermal_compressibility_reuss * ( + dPsidf_Voigt + dPsidPth_Voigt * self.dPthdf ) + # Calculate F, dFdP, dFdT + self._F = expm(PsiI) + + if self.orthotropic: + self._S_T_unrotated_Voigt = -dPsidP_Voigt + + self._unrotated_alpha = self.alpha * ( + dPsidfI + dPsidPthI * (self.dPthdf + self.isothermal_bulk_modulus_reuss) + ) + else: + # Numerical derivatives with respect to f and Pth + df = f * 1.0e-5 + Psi_full0 = voigt_notation_to_compliance_tensor( + Psi_Voigt - dPsidf_Voigt * df / 2.0 + ) + Psi_full1 = voigt_notation_to_compliance_tensor( + Psi_Voigt + dPsidf_Voigt * df / 2.0 + ) + F0 = expm(np.einsum("ijkl, kl", Psi_full0, np.eye(3))) + F1 = expm(np.einsum("ijkl, kl", Psi_full1, np.eye(3))) + dFdf = (F1 - F0) / df + + dPth = self.Pth * 1.0e-5 + Psi_full0 = voigt_notation_to_compliance_tensor( + Psi_Voigt - dPsidPth_Voigt * dPth / 2.0 + ) + Psi_full1 = voigt_notation_to_compliance_tensor( + Psi_Voigt + dPsidPth_Voigt * dPth / 2.0 + ) + F0 = expm(np.einsum("ijkl, kl", Psi_full0, np.eye(3))) + F1 = expm(np.einsum("ijkl, kl", Psi_full1, np.eye(3))) + dFdPth = (F1 - F0) / dPth + + # Convert to pressure and temperature derivatives + dFdP = -self.isothermal_compressibility_reuss * ( + dFdf + dFdPth * self.dPthdf + ) + dFdT = self.alpha * ( + dFdf + dFdPth * (self.dPthdf + self.isothermal_bulk_modulus_reuss) + ) + + # Calculate the unrotated isothermal compressibility + # and unrotated thermal expansivity tensors + invF = np.linalg.inv(self._F) + LP = np.einsum("ij,kj->ik", dFdP, invF) + beta_T = -0.5 * (LP + LP.T) + LT = np.einsum("ij,kj->ik", dFdT, invF) + self._unrotated_alpha = 0.5 * (LT + LT.T) + + # Calculate the unrotated isothermal compliance + # tensor in Voigt form. + S = -dPsidP_Voigt + for i, j, k in [[0, 1, 2], [1, 2, 0], [0, 2, 1]]: + S[i][j] = 0.5 * ( + -S[i][i] + - S[j][j] + + S[k][k] + + beta_T[i][i] + + beta_T[j][j] + - beta_T[k][k] + ) + S[j][i] = S[i][j] + + S[i][i + 3] = 2.0 * beta_T[j][k] - S[j][i + 3] - S[k][i + 3] + S[i + 3][i] = S[i][i + 3] + + self._S_T_unrotated_Voigt = S + @material_property def deformation_gradient_tensor(self): """ @@ -186,8 +252,7 @@ def deformation_gradient_tensor(self): (i.e. the state at the reference pressure and temperature). :rtype: numpy.array (2D) """ - F = expm(np.einsum("ijkl, kl", self.Psi, np.eye(3))) - return F + return self._F @material_property def unrotated_cell_vectors(self): @@ -201,7 +266,9 @@ def unrotated_cell_vectors(self): spatial coordinate axes. :rtype: numpy.array (2D) """ - return self.deformation_gradient_tensor + return np.einsum( + "ij,kj->ki", self.deformation_gradient_tensor, self.cell_vectors_0 + ) @material_property def deformed_coordinate_frame(self): @@ -215,10 +282,10 @@ def deformed_coordinate_frame(self): if self.orthotropic: return np.eye(3) else: - M = self.unrotated_cell_vectors + M_T = self.unrotated_cell_vectors Q = np.empty((3, 3)) - Q[0] = M[0] / np.linalg.norm(M[0]) - Q[2] = np.cross(M[0], M[1]) / np.linalg.norm(np.cross(M[0], M[1])) + Q[0] = M_T[0] / np.linalg.norm(M_T[0]) + Q[2] = np.cross(M_T[0], M_T[1]) / np.linalg.norm(np.cross(M_T[0], M_T[1])) Q[1] = np.cross(Q[2], Q[0]) return Q @@ -409,12 +476,11 @@ def isothermal_compliance_tensor(self): in Voigt form (:math:`\\mathbb{S}_{\\text{T} pq}`). :rtype: numpy.array (2D) """ - S_T = -self.dPsidP_Voigt if self.orthotropic: - return S_T + return self._S_T_unrotated_Voigt else: R = self.rotation_matrix - S = voigt_notation_to_compliance_tensor(S_T) + S = voigt_notation_to_compliance_tensor(self._S_T_unrotated_Voigt) S_rotated = np.einsum("mi, nj, ok, pl, ijkl->mnop", R, R, R, R, S) return contract_compliances(S_rotated) @@ -424,17 +490,11 @@ def thermal_expansivity_tensor(self): :returns: The tensor of thermal expansivities [1/K]. :rtype: numpy.array (2D) """ - alpha = np.einsum( - "ijkl, kl", - voigt_notation_to_compliance_tensor(self.dPsidT_Voigt), - np.eye(3), - ) - if self.orthotropic: - return alpha + return self._unrotated_alpha else: R = self.rotation_matrix - return np.einsum("mi, nj, ij->mn", R, R, alpha) + return np.einsum("mi, nj, ij->mn", R, R, self._unrotated_alpha) # Derived properties start here @material_property @@ -611,7 +671,7 @@ def check_standard_parameters(self, anisotropic_parameters): if np.abs(sum_ijij_block[m, 0]) > 1.0e-10: raise Exception( "The sum of the upper 3x3 pq-block of " - "anisotropic_parameters_pqmn must equal 0 for" + "anisotropic_parameters_pqmn must equal 0 for " f"m={m}, n=0 for consistency with the volume. " f"Value is {sum_ijij_block[m, 0]}" ) @@ -627,6 +687,14 @@ def check_standard_parameters(self, anisotropic_parameters): f"Value is {sum_ijij_block[m, n]}" ) + for m in range(len(sum_ijij_block)): + for n in range(len(sum_ijij_block[0])): + a = anisotropic_parameters[:, :, m, n] + if not np.allclose(a, a.T, rtol=1.0e-8, atol=1.0e-8): + raise Exception( + f"The anisotropic_parameters_pq{m}{n} must be symmetric." + ) + if cond(anisotropic_parameters[:, :, 1, 0]) > 1 / np.finfo(float).eps: raise Exception("anisotropic_parameters[:, :, 1, 0] is singular") diff --git a/burnman/classes/averaging_schemes.py b/burnman/classes/averaging_schemes.py index 0a8a3edf5..71971e960 100644 --- a/burnman/classes/averaging_schemes.py +++ b/burnman/classes/averaging_schemes.py @@ -9,7 +9,6 @@ class AveragingScheme(object): - """ Base class defining an interface for determining average elastic properties of a rock. Given a list of volume @@ -144,7 +143,6 @@ def average_heat_capacity_p(self, fractions, c_p): class VoigtReussHill(AveragingScheme): - """ Class for computing the Voigt-Reuss-Hill average for elastic properties. This derives from :class:`burnman.averaging_schemes.averaging_scheme`, @@ -218,7 +216,6 @@ def average_shear_moduli(self, volumes, bulk_moduli, shear_moduli): class Voigt(AveragingScheme): - """ Class for computing the Voigt (iso-strain) bound for elastic properties. This derives from :class:`burnman.averaging_schemes.averaging_scheme`, @@ -277,7 +274,6 @@ def average_shear_moduli(self, volumes, bulk_moduli, shear_moduli): class Reuss(AveragingScheme): - """ Class for computing the Reuss (iso-stress) bound for elastic properties. This derives from :class:`burnman.averaging_schemes.averaging_scheme`, @@ -336,7 +332,6 @@ def average_shear_moduli(self, volumes, bulk_moduli, shear_moduli): class HashinShtrikmanUpper(AveragingScheme): - """ Class for computing the upper Hashin-Shtrikman bound for elastic properties. This derives from :class:`burnman.averaging_schemes.averaging_scheme`, @@ -423,7 +418,6 @@ def average_shear_moduli(self, volumes, bulk_moduli, shear_moduli): class HashinShtrikmanLower(AveragingScheme): - """ Class for computing the lower Hashin-Shtrikman bound for elastic properties. This derives from :class:`burnman.averaging_schemes.averaging_scheme`, @@ -510,7 +504,6 @@ def average_shear_moduli(self, volumes, bulk_moduli, shear_moduli): class HashinShtrikmanAverage(AveragingScheme): - """ Class for computing arithmetic mean of the Hashin-Shtrikman bounds on elastic properties. diff --git a/burnman/classes/combinedmineral.py b/burnman/classes/combinedmineral.py index 47aead982..091a38549 100644 --- a/burnman/classes/combinedmineral.py +++ b/burnman/classes/combinedmineral.py @@ -14,7 +14,6 @@ class CombinedMineral(Mineral): - """ This is the base class for endmembers constructed from a linear combination of other minerals. diff --git a/burnman/classes/composite.py b/burnman/classes/composite.py index f992afefe..7b18b9b95 100644 --- a/burnman/classes/composite.py +++ b/burnman/classes/composite.py @@ -38,7 +38,6 @@ def check_pairs(phases, fractions): # static composite of minerals/composites class Composite(Material): - """ Base class for a composite material. The static phases can be minerals or materials, diff --git a/burnman/classes/elasticsolutionmodel.py b/burnman/classes/elasticsolutionmodel.py index caae99a70..bbebe0e88 100644 --- a/burnman/classes/elasticsolutionmodel.py +++ b/burnman/classes/elasticsolutionmodel.py @@ -26,7 +26,6 @@ class ElasticSolutionModel(object): - """ This is the base class for an Elastic solution model, intended for use in defining solutions and performing thermodynamic calculations @@ -216,7 +215,6 @@ def excess_partial_pressures(self, volume, temperature, molar_fractions): class ElasticMechanicalSolution(ElasticSolutionModel): - """ An extremely simple class representing a mechanical solution model. A mechanical solution experiences no interaction between endmembers. @@ -251,7 +249,6 @@ def excess_partial_entropies(self, volume, temperature, molar_fractions): class ElasticIdealSolution(ElasticSolutionModel): - """ A class representing an ideal solution model. Calculates the excess Helmholtz energy and entropy due to configurational @@ -382,7 +379,6 @@ def _ideal_activities(self, molar_fractions): class ElasticAsymmetricRegularSolution(ElasticIdealSolution): - """ Solution model implementing the asymmetric regular solution model formulation as described in :cite:`HP2003`. @@ -510,7 +506,6 @@ def pressure_hessian(self, volume, temperature, molar_fractions): class ElasticSymmetricRegularSolution(ElasticAsymmetricRegularSolution): - """ Solution model implementing the symmetric regular solution model. This is a special case of the @@ -536,7 +531,6 @@ def __init__( class ElasticSubregularSolution(ElasticIdealSolution): - """ Solution model implementing the subregular solution model formulation as described in :cite:`HW1989`. The excess conconfigurational diff --git a/burnman/classes/material.py b/burnman/classes/material.py index 71f8088e5..98f30c0ee 100644 --- a/burnman/classes/material.py +++ b/burnman/classes/material.py @@ -14,7 +14,6 @@ # TODO: When we require Python 3.8+, replace with # functools.cached_property decorator class cached_property(property): - """A decorator that converts a function into a lazy property. The function wrapped is called the first time to retrieve the result and then that calculated result is used the next time you access @@ -87,7 +86,6 @@ def get(self, obj): class Material(object): - """ Base class for all materials. The main functionality is unroll() which returns a list of objects of type :class:`~burnman.Mineral` and their molar diff --git a/burnman/classes/mineral.py b/burnman/classes/mineral.py index fdeb39fbb..f7507402d 100644 --- a/burnman/classes/mineral.py +++ b/burnman/classes/mineral.py @@ -15,7 +15,6 @@ class Mineral(Material): - """ This is the base class for all minerals. States of the mineral can only be queried after setting the pressure and temperature diff --git a/burnman/classes/mineral_helpers.py b/burnman/classes/mineral_helpers.py index 1a7e5c29e..5fc4c8919 100644 --- a/burnman/classes/mineral_helpers.py +++ b/burnman/classes/mineral_helpers.py @@ -168,7 +168,6 @@ def debug_print(self, indent=""): class HelperSpinTransition(Composite): - """ Helper class that makes a mineral that switches between two materials (for low and high spin) based on some transition pressure [Pa] diff --git a/burnman/classes/solutionmodel.py b/burnman/classes/solutionmodel.py index cd97682ef..e92cabee7 100644 --- a/burnman/classes/solutionmodel.py +++ b/burnman/classes/solutionmodel.py @@ -128,7 +128,6 @@ def inverseish(x, eps=1.0e-5): class SolutionModel(object): - """ This is the base class for a solution model, intended for use in defining solutions and performing thermodynamic calculations @@ -326,7 +325,6 @@ def VoverKT_excess(self): class MechanicalSolution(SolutionModel): - """ An extremely simple class representing a mechanical solution model. A mechanical solution experiences no interaction between endmembers. @@ -372,7 +370,6 @@ def activities(self, pressure, temperature, molar_fractions): class IdealSolution(SolutionModel): - """ A class representing an ideal solution model. Calculates the excess gibbs free energy and entropy due to configurational @@ -511,7 +508,6 @@ def activities(self, pressure, temperature, molar_fractions): class AsymmetricRegularSolution(IdealSolution): - """ Solution model implementing the asymmetric regular solution model formulation as described in :cite:`HP2003`. @@ -657,7 +653,6 @@ def activities(self, pressure, temperature, molar_fractions): class SymmetricRegularSolution(AsymmetricRegularSolution): - """ Solution model implementing the symmetric regular solution model. This is a special case of the @@ -683,7 +678,6 @@ def __init__( class SubregularSolution(IdealSolution): - """ Solution model implementing the subregular solution model formulation as described in :cite:`HW1989`. The excess nonconfigurational diff --git a/burnman/eos/birch_murnaghan.py b/burnman/eos/birch_murnaghan.py index d4bf98a8b..14f2758a8 100644 --- a/burnman/eos/birch_murnaghan.py +++ b/burnman/eos/birch_murnaghan.py @@ -114,7 +114,6 @@ def shear_modulus_third_order(volume, params): class BirchMurnaghanBase(eos.EquationOfState): - """ Base class for the isothermal Birch Murnaghan equation of state. This is third order in strain, and has no temperature dependence. However, the shear modulus is sometimes fit to a second order @@ -258,7 +257,6 @@ def validate_parameters(self, params): class BM3(BirchMurnaghanBase): - """ Third order Birch Murnaghan isothermal equation of state. This uses the third order expansion for shear modulus. @@ -269,7 +267,6 @@ def __init__(self): class BM2(BirchMurnaghanBase): - """ Third order Birch Murnaghan isothermal equation of state. This uses the second order expansion for shear modulus. diff --git a/burnman/eos/birch_murnaghan_4th.py b/burnman/eos/birch_murnaghan_4th.py index b658783af..3d04aab7d 100644 --- a/burnman/eos/birch_murnaghan_4th.py +++ b/burnman/eos/birch_murnaghan_4th.py @@ -81,7 +81,6 @@ def birch_murnaghan_fourth(x, params): class BM4(eos.EquationOfState): - """ Base class for the isothermal Birch Murnaghan equation of state. This is fourth order in strain, and has no temperature dependence. diff --git a/burnman/eos/cork.py b/burnman/eos/cork.py index c71954f96..cbd3f51f6 100644 --- a/burnman/eos/cork.py +++ b/burnman/eos/cork.py @@ -34,7 +34,6 @@ def cork_variables(cork, cork_P, cork_T, temperature): class CORK(eos.EquationOfState): - """ Class for the CoRK equation of state detailed in :cite:`HP1991`. The CoRK EoS is a simple virial-type extension to the modified Redlich-Kwong diff --git a/burnman/eos/dks_solid.py b/burnman/eos/dks_solid.py index d4a737ceb..09b9ed8d5 100644 --- a/burnman/eos/dks_solid.py +++ b/burnman/eos/dks_solid.py @@ -61,7 +61,6 @@ def _delta_pressure( class DKS_S(eos.EquationOfState): - """ Base class for the finite strain solid equation of state detailed in :cite:`deKoker2013` (supplementary materials). diff --git a/burnman/eos/equation_of_state.py b/burnman/eos/equation_of_state.py index e488a1700..901c3515d 100644 --- a/burnman/eos/equation_of_state.py +++ b/burnman/eos/equation_of_state.py @@ -5,7 +5,6 @@ class EquationOfState(object): - """ This class defines the interface for an equation of state that a mineral uses to determine its properties at a diff --git a/burnman/eos/hp.py b/burnman/eos/hp.py index 252ca57dc..c5021a781 100644 --- a/burnman/eos/hp.py +++ b/burnman/eos/hp.py @@ -17,7 +17,6 @@ class HP_TMT(eos.EquationOfState): - """ Base class for the thermal equation of state based on the generic modified Tait equation of state (class MT), @@ -395,7 +394,6 @@ def validate_parameters(self, params): class HP_TMTL(eos.EquationOfState): - """ Base class for the thermal equation of state described in :cite:`HP1998`, but with the Modified Tait as the static part, @@ -668,7 +666,6 @@ def validate_parameters(self, params): class HP98(eos.EquationOfState): - """ Base class for the thermal equation of state described in :cite:`HP1998`. diff --git a/burnman/eos/mie_grueneisen_debye.py b/burnman/eos/mie_grueneisen_debye.py index d29116712..0eab09989 100644 --- a/burnman/eos/mie_grueneisen_debye.py +++ b/burnman/eos/mie_grueneisen_debye.py @@ -16,7 +16,6 @@ class MGDBase(eos.EquationOfState): - """ Base class for a generic finite-strain Mie-Grueneisen-Debye equation of state. References for this can be found in many @@ -304,7 +303,6 @@ def validate_parameters(self, params): class MGD3(MGDBase): - """ MGD equation of state with third order finite strain expansion for the shear modulus (this should be preferred, as it is more thermodynamically @@ -316,7 +314,6 @@ def __init__(self): class MGD2(MGDBase): - """ MGD equation of state with second order finite strain expansion for the shear modulus. In general, this should not be used, but sometimes diff --git a/burnman/eos/modified_tait.py b/burnman/eos/modified_tait.py index a2d9fb89b..b5290fbe9 100644 --- a/burnman/eos/modified_tait.py +++ b/burnman/eos/modified_tait.py @@ -92,7 +92,6 @@ def intVdP(pressure, params): class MT(eos.EquationOfState): - """ Base class for the generic modified Tait equation of state. References for this can be found in :cite:`HC1974` diff --git a/burnman/eos/morse_potential.py b/burnman/eos/morse_potential.py index bc6daf336..e80b865b1 100644 --- a/burnman/eos/morse_potential.py +++ b/burnman/eos/morse_potential.py @@ -75,7 +75,6 @@ def volume(pressure, params): class Morse(eos.EquationOfState): - """ Class for the isothermal Morse Potential equation of state detailed in :cite:`Stacey1981`. diff --git a/burnman/eos/murnaghan.py b/burnman/eos/murnaghan.py index b3c9485c9..3d77c9a2e 100644 --- a/burnman/eos/murnaghan.py +++ b/burnman/eos/murnaghan.py @@ -39,7 +39,6 @@ def intVdP(pressure, V_0, K_0, Kprime_0): class Murnaghan(eos.EquationOfState): - """ Base class for the isothermal Murnaghan equation of state, as described in :cite:`Murnaghan1944`. diff --git a/burnman/eos/reciprocal_kprime.py b/burnman/eos/reciprocal_kprime.py index ea60439b7..201f839a9 100644 --- a/burnman/eos/reciprocal_kprime.py +++ b/burnman/eos/reciprocal_kprime.py @@ -124,7 +124,6 @@ def shear_modulus(pressure, params): class RKprime(eos.EquationOfState): - """ Class for the isothermal reciprocal K-prime equation of state detailed in :cite:`StaceyDavis2004`. This equation of state is diff --git a/burnman/eos/slb.py b/burnman/eos/slb.py index a7e515940..b58bd46c1 100644 --- a/burnman/eos/slb.py +++ b/burnman/eos/slb.py @@ -63,7 +63,6 @@ def _delta_pressure( class SLBBase(eos.EquationOfState): - """ Base class for the finite strain-Mie-Grueneiesen-Debye equation of state detailed in :cite:`Stixrude2005`. For the most part the equations are @@ -468,7 +467,6 @@ def validate_parameters(self, params): class SLB3(SLBBase): - """ SLB equation of state with third order finite strain expansion for the shear modulus (this should be preferred, as it is more thermodynamically @@ -480,7 +478,6 @@ def __init__(self): class SLB2(SLBBase): - """ SLB equation of state with second order finite strain expansion for the shear modulus. In general, this should not be used, but sometimes diff --git a/burnman/eos/vinet.py b/burnman/eos/vinet.py index 829a2eac0..910100cbc 100644 --- a/burnman/eos/vinet.py +++ b/burnman/eos/vinet.py @@ -56,7 +56,6 @@ def volume(pressure, params): class Vinet(eos.EquationOfState): - """ Base class for the isothermal Vinet equation of state. References for this equation of state are :cite:`vinet1986` diff --git a/burnman/minerals/Murakami_2013.py b/burnman/minerals/Murakami_2013.py index c8c497a40..8fd6feff5 100644 --- a/burnman/minerals/Murakami_2013.py +++ b/burnman/minerals/Murakami_2013.py @@ -36,7 +36,6 @@ def __init__(self): class wuestite(Mineral): - """ Murakami 2013 and references therein """ diff --git a/burnman/minerals/other.py b/burnman/minerals/other.py index 7ac0aecdb..c094648a7 100644 --- a/burnman/minerals/other.py +++ b/burnman/minerals/other.py @@ -130,7 +130,6 @@ def __init__(self): class Speziale_fe_periclase_HS(Mineral): - """ Speziale et al. 2007, Mg#=83 """ @@ -151,7 +150,6 @@ def __init__(self): class Speziale_fe_periclase_LS(Mineral): - """ Speziale et al. 2007, Mg#=83 """ @@ -172,7 +170,6 @@ def __init__(self): class Liquid_Fe_Anderson(Mineral): - """ Anderson & Ahrens, 1994 JGR """ @@ -190,7 +187,6 @@ def __init__(self): class Fe_Dewaele(Mineral): - """ Dewaele et al., 2006, Physical Review Letters """ diff --git a/burnman/tools/eos.py b/burnman/tools/eos.py index 1d6aeb802..64d88f8d9 100644 --- a/burnman/tools/eos.py +++ b/burnman/tools/eos.py @@ -260,12 +260,37 @@ def check_anisotropic_eos_consistency(m, P=1.0e9, T=2000.0, tol=1.0e-4, verbose= beta1 = np.einsum("mi, nj, ij->mn", Q, Q, beta1) alpha1 = np.einsum("mi, nj, ij->mn", Q, Q, alpha1) - expr.extend([f"SI = -d(lnm(F))/dP ({i}{j})" for i in range(3) for j in range(i, 3)]) - eq.extend([[beta0[i, j], beta1[i, j]] for i in range(3) for j in range(i, 3)]) + if m.orthotropic: + expr.extend( + [f"SI = -d(lnm(F))/dP ({i}{j})" for i in range(3) for j in range(i, 3)] + ) + expr.extend( + [f"alpha = d(lnm(F))/dT ({i}{j})" for i in range(3) for j in range(i, 3)] + ) + else: + expr.extend( + [ + f"SI = -0.5(dF/dP*F^-1 + (dF/dP*F^-1)^T) ({i}{j})" + for i in range(3) + for j in range(i, 3) + ] + ) + expr.extend( + [ + f"alpha = 0.5(dF/dT*F^-1 + (dF/dT*F^-1)^T) ({i}{j})" + for i in range(3) + for j in range(i, 3) + ] + ) + invF = np.linalg.inv(m.deformation_gradient_tensor) + dFdP = (F3 - F2) / dP + LP = np.einsum("ij,kj->ik", dFdP, invF) + beta0 = -0.5 * (LP + LP.T) + dFdT = (F1 - F0) / dT + LT = np.einsum("ij,kj->ik", dFdT, invF) + alpha0 = 0.5 * (LT + LT.T) - expr.extend( - [f"alpha = d(lnm(F))/dT ({i}{j})" for i in range(3) for j in range(i, 3)] - ) + eq.extend([[beta0[i, j], beta1[i, j]] for i in range(3) for j in range(i, 3)]) eq.extend([[alpha0[i, j], alpha1[i, j]] for i in range(3) for j in range(i, 3)]) expr.extend( diff --git a/burnman/tools/equilibration.py b/burnman/tools/equilibration.py index ec55330a4..ec59c4146 100644 --- a/burnman/tools/equilibration.py +++ b/burnman/tools/equilibration.py @@ -644,9 +644,11 @@ def get_equilibration_parameters(assemblage, composition, free_compositional_vec prm.free_compositional_vectors = np.array( [ [ - free_compositional_vectors[i][e] - if e in free_compositional_vectors[i] - else 0.0 + ( + free_compositional_vectors[i][e] + if e in free_compositional_vectors[i] + else 0.0 + ) for e in assemblage.elements ] for i in range(n_free_compositional_vectors) diff --git a/contrib/CHRU2014/helper_solid_solution.py b/contrib/CHRU2014/helper_solid_solution.py index eef58957f..e243b1280 100644 --- a/contrib/CHRU2014/helper_solid_solution.py +++ b/contrib/CHRU2014/helper_solid_solution.py @@ -12,7 +12,6 @@ class HelperSolidSolution(burnman.Mineral): - """ This material is deprecated! diff --git a/contrib/CHRU2014/paper_uncertain.py b/contrib/CHRU2014/paper_uncertain.py index 2fbc7e0a4..60cb67741 100644 --- a/contrib/CHRU2014/paper_uncertain.py +++ b/contrib/CHRU2014/paper_uncertain.py @@ -28,7 +28,6 @@ class my_perovskite(burnman.Mineral): - """ based on Stixrude & Lithgow-Bertelloni 2011 and references therein """ diff --git a/contrib/anisotropic_eos/tools.py b/contrib/anisotropic_eos/tools.py index eab6050cb..0da41a1cf 100644 --- a/contrib/anisotropic_eos/tools.py +++ b/contrib/anisotropic_eos/tools.py @@ -42,9 +42,11 @@ def print_table_for_mineral_constants_2(mineral, param_list, indices): row = [f"{i}", f"{j}"] row.extend( [ - f"{constants[ci, i]:.4e}" - if constants[ci, i] != 0 and constants[ci, i] != 1 - else "-" + ( + f"{constants[ci, i]:.4e}" + if constants[ci, i] != 0 and constants[ci, i] != 1 + else "-" + ) for i in range(len(param_list)) ] ) diff --git a/examples/example_user_input_material.py b/examples/example_user_input_material.py index 2f49ddfc3..aeae81b29 100644 --- a/examples/example_user_input_material.py +++ b/examples/example_user_input_material.py @@ -81,7 +81,7 @@ def __init__(self): # See Stixrude & Lithgow-Bertelloni, 2005 for values "q_0": 0.917, # isotropic strain derivative of gruneisen # parameter. Values in Stixrude & Lithgow-Bertelloni, 2005 - "eta_s_0": 3.0 # full strain derivative of gruneisen parameter + "eta_s_0": 3.0, # full strain derivative of gruneisen parameter # parameter. Values in Stixrude & Lithgow-Bertelloni, 2005 } burnman.Mineral.__init__(self) diff --git a/misc/gen_doc.py b/misc/gen_doc.py index 74bdf8107..8d7b4b0a2 100644 --- a/misc/gen_doc.py +++ b/misc/gen_doc.py @@ -1,4 +1,5 @@ """ generates a list with the examples """ + from __future__ import absolute_import from __future__ import print_function diff --git a/test.sh b/test.sh index de3a87529..4d4968b6d 100755 --- a/test.sh +++ b/test.sh @@ -76,6 +76,7 @@ else sed -i.bak -e '/UserWarning: findfont: Font family/d' $t.tmp #remove font warning crap sed -i.bak -e '/tight_layout : falling back to Agg renderer/d' $t.tmp #remove font warning crap sed -i.bak -e '/cannot be converted with the encoding. Glyph may be wrong/d' $t.tmp #remove font warning crap + sed -i.bak -e '/UserWarning: FigureCanvasTemplate/d' $t.tmp #remove plotting nonsense sed -i.bak -e '/time old .* time new/d' $t.tmp #remove timing from tests/debye.py sed -i.bak -e '/ relative central core pressure error.*/d' $t.tmp #remove residuals from examples/example_build_planet diff --git a/tests/test_anisotropicmineral.py b/tests/test_anisotropicmineral.py index c2c746096..9048a9753 100644 --- a/tests/test_anisotropicmineral.py +++ b/tests/test_anisotropicmineral.py @@ -8,24 +8,40 @@ from burnman.minerals.SLB_2011 import periclase, forsterite -def make_simple_forsterite(): +def make_forsterite(orthotropic=True): fo = forsterite() cell_lengths = np.array([4.7646, 10.2296, 5.9942]) cell_lengths *= np.cbrt(fo.params["V_0"] / np.prod(cell_lengths)) - cell_lengths *= 1.0710965273664157 + + if orthotropic: + alpha, beta, gamma, a, b, c = [90.0, 90.0, 90.0, 0.0, 0.0, 0.0] + else: + alpha, beta, gamma, a, b, c = [85.0, 80.0, 87.0, 0.4, -1.0, -0.6] + cell_lengths *= 1.006635538793111 + cell_parameters = np.array( - [cell_lengths[0], cell_lengths[1], cell_lengths[2], 60, 70, 80] + [cell_lengths[0], cell_lengths[1], cell_lengths[2], alpha, beta, gamma] ) - constants = np.zeros((6, 6, 2, 1)) + constants = np.zeros((6, 6, 3, 1)) constants[:, :, 1, 0] = np.array( [ - [0.44, -0.12, -0.1, 1.0, 0.5, 0.5], + [0.44, -0.12, -0.1, a, b, c], [-0.12, 0.78, -0.22, 0.0, 0.0, 0.0], [-0.1, -0.22, 0.66, 0.0, 0.0, 0.0], - [1.0, 0.0, 0.0, 1.97, 0.0, 0.0], - [0.5, 0.0, 0.0, 0.0, 1.61, 0.0], - [0.5, 0.0, 0.0, 0.0, 0.0, 1.55], + [a, 0.0, 0.0, 1.97, 0.0, 0.0], + [b, 0.0, 0.0, 0.0, 1.61, 0.0], + [c, 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, a, a, a], + [-0.1, -0.22, 0.26, c, b, a], + [0.0, a, c, 0.0, 0.0, 0.0], + [0.0, a, b, 0.0, 0.0, 0.0], + [0.0, a, a, 0.0, 0.0, 0.0], ] ) @@ -61,25 +77,30 @@ def test_isotropic_grueneisen(self): gr = per.grueneisen_parameter self.assertArraysAlmostEqual(np.diag(per2.grueneisen_tensor), [gr, gr, gr]) - def test_orthorhombic_consistency(self): - m = make_simple_forsterite() + def test_orthotropic_consistency(self): + m = make_forsterite(orthotropic=True) self.assertTrue(check_anisotropic_eos_consistency(m)) + def test_non_orthotropic_consistency(self): + m = make_forsterite(orthotropic=False) + self.assertTrue(check_anisotropic_eos_consistency(m, P=2.0e10, tol=1.0e-6)) + def test_stiffness(self): - m = make_simple_forsterite() - m.set_state(1.0e9, 300.0) - Cijkl = m.full_isothermal_stiffness_tensor - Cij = m.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]) + for orthotropic in [True, False]: + m = make_forsterite(orthotropic) + m.set_state(1.0e9, 300.0) + Cijkl = m.full_isothermal_stiffness_tensor + Cij = m.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__": diff --git a/tests/test_averaging.py b/tests/test_averaging.py index 7d6f95467..a257ccd33 100644 --- a/tests/test_averaging.py +++ b/tests/test_averaging.py @@ -9,7 +9,6 @@ class mypericlase(burnman.Mineral): - """ Stixrude & Lithgow-Bertelloni 2005 and references therein """ diff --git a/tests/test_debye.py b/tests/test_debye.py index a02fb5932..0ed8436e8 100644 --- a/tests/test_debye.py +++ b/tests/test_debye.py @@ -6,7 +6,6 @@ class mypericlase(burnman.Mineral): - """ Stixrude & Lithgow-Bertelloni 2005 and references therein """ diff --git a/tests/test_eos.py b/tests/test_eos.py index 86f024f29..c86ec9c03 100644 --- a/tests/test_eos.py +++ b/tests/test_eos.py @@ -11,7 +11,6 @@ class mypericlase(burnman.Mineral): - """ Stixrude & Lithgow-Bertelloni 2005 and references therein """ @@ -36,7 +35,6 @@ def __init__(self): class Fe_Dewaele(burnman.Mineral): - """ Dewaele et al., 2006, Physical Review Letters """ @@ -54,7 +52,6 @@ def __init__(self): class Liquid_Fe_Anderson(burnman.Mineral): - """ Anderson & Ahrens, 1994 JGR """ @@ -72,7 +69,6 @@ def __init__(self): class outer_core_rkprime(burnman.Mineral): - """ Stacey and Davis, 2004 PEPI (Table 5) """ @@ -90,7 +86,6 @@ def __init__(self): class periclase_morse(burnman.Mineral): - """ Periclase parameters from SLB dataset (which uses BM3) """ diff --git a/tests/test_geotherm.py b/tests/test_geotherm.py index 34d17d4db..c82557bb6 100644 --- a/tests/test_geotherm.py +++ b/tests/test_geotherm.py @@ -6,7 +6,6 @@ class mypericlase(burnman.Mineral): - """ Stixrude & Lithgow-Bertelloni 2005 and references therein """ diff --git a/tests/test_material.py b/tests/test_material.py index 7e0ea89db..8f142b27d 100644 --- a/tests/test_material.py +++ b/tests/test_material.py @@ -7,11 +7,9 @@ class test_material_name(BurnManTest): - """test Material.name and that we can edit and override it in Mineral""" class min_no_name(burnman.Mineral): - """ Stixrude & Lithgow-Bertelloni 2005 and references therein """ @@ -36,7 +34,6 @@ def __init__(self): burnman.Mineral.__init__(self) class min_with_name(burnman.Mineral): - """ Stixrude & Lithgow-Bertelloni 2005 and references therein """ @@ -62,7 +59,6 @@ def __init__(self): burnman.Mineral.__init__(self) class min_with_name_manually(burnman.Mineral): - """ Stixrude & Lithgow-Bertelloni 2005 and references therein """