diff --git a/burnman/classes/anisotropicmineral.py b/burnman/classes/anisotropicmineral.py index 872ca6b5..6c1d1c9a 100644 --- a/burnman/classes/anisotropicmineral.py +++ b/burnman/classes/anisotropicmineral.py @@ -18,6 +18,119 @@ ) +def convert_f_Pth_to_f_T_derivatives( + dPsidf_Pth_Voigt, dPsidPth_f_Voigt, aK_T, dPthdf_T +): + """ + Convenience function converting Psi and its derivatives with respect to + f and Pth into derivatives of f and T using the chain rule. + + :param dPsidf_Pth_Voigt: The first derivative of the anisotropic tensor + Psi with respect to f at constant Pth (Voigt form). + :type dPsiIdf_Pth: numpy array (6x6) + :param dPsidPth_f_Voigt: The first derivative of the anisotropic tensor + Psi with respect to Pth at constant f (Voigt form). + :type dPsiIdPth_f: numpy array (6x6) + :param aK_T: The volumetric thermal expansivity multiplied by + the Reuss isothermal bulk modulus. + :type aK_T: float + :param dPthdf: The change in thermal pressure with respect to volume + at constant temperature. + :type dPthdf: float + + :returns: The first derivative of Psi with respect to f in Voigt form, + and the first derivatives of PsiI with respect to f and T. + :rtype: Tuple of three objects of type numpy.array (2D) + """ + # The following manipulation uses the chain rule to convert + # from Psi(f, Pth) and derivatives to Psi(f, T) and derivatives + dPsidf_full = voigt_notation_to_compliance_tensor(dPsidf_Pth_Voigt) + dPsidPth_full = voigt_notation_to_compliance_tensor(dPsidPth_f_Voigt) + + dPsiIdf_Pth = np.einsum("ijkl, kl", dPsidf_full, np.eye(3)) + dPsiIdPth_f = np.einsum("ijkl, kl", dPsidPth_full, np.eye(3)) + + dPsidf_T_Voigt = dPsidf_Pth_Voigt + dPsidPth_f_Voigt * dPthdf_T + dPsiIdf_T = dPsiIdf_Pth + dPsiIdPth_f * dPthdf_T + dPsiIdT_f = aK_T * dPsiIdPth_f + return (dPsidf_T_Voigt, dPsiIdf_T, dPsiIdT_f) + + +def deformation_gradient_alpha_and_compliance( + alpha_V, beta_TR, PsiI, dPsidf_T_Voigt, dPsiIdf_T, dPsiIdT_f +): + """ + Convenience function converting Psi and its derivatives with respect to + f and T into the deformation gradient, thermal expansivity and + isothermal compliance tensors. This function is appropriate for + both orthotropic and non-orthotropic materials. + + :param alpha_V: The volumetric thermal expansivity. + :type alpha_V: float + :param beta_TR: The Reuss isothermal compressibility. + :type beta_TR: float + :param psiI: The anisotropic tensor Psi matrix multiplied with + the identity matrix. + :type psiI: numpy array (3x3) + :param dPsidf_T_Voigt: Voigt-form matrix of the first + derivative of the anisotropic tensor Psi with respect to f + at constant T. + :type dPsidf_T_Voigt: numpy array (6x6) + :param dPsiIdf_T: The first derivative of PsiI + with respect to f at constant T. + :type dPsiIdf_T: numpy array (3x3) + :param dPsiIdT_f: The first derivative of PsiI + 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) + """ + # Numerical derivatives with respect to f and T + df = 1.0e-7 + F0 = expm(PsiI - dPsiIdf_T * df / 2.0) + F1 = expm(PsiI + dPsiIdf_T * df / 2.0) + dFdf_T = (F1 - F0) / df + + dT = 0.1 + F0 = expm(PsiI - dPsiIdT_f * dT / 2.0) + F1 = expm(PsiI + dPsiIdT_f * dT / 2.0) + dFdT_f = (F1 - F0) / dT + + # Convert to pressure and temperature derivatives + dFdP_T = -beta_TR * dFdf_T + dFdT_P = dFdT_f + alpha_V * dFdf_T + + # Calculate the unrotated isothermal compressibility + # and unrotated thermal expansivity tensors + F = expm(PsiI) + invF = np.linalg.inv(F) + LP = np.einsum("ij,kj->ik", dFdP_T, invF) + beta_T = -0.5 * (LP + LP.T) + LT = np.einsum("ij,kj->ik", dFdT_P, invF) + alpha = 0.5 * (LT + LT.T) + + # Calculate the unrotated isothermal compliance + # tensor in Voigt form. + S_T = beta_TR * dPsidf_T_Voigt + for i, j, k in [[0, 1, 2], [1, 2, 0], [0, 2, 1]]: + S_T[i][j] = 0.5 * ( + -S_T[i][i] + - S_T[j][j] + + S_T[k][k] + + beta_T[i][i] + + beta_T[j][j] + - beta_T[k][k] + ) + S_T[j][i] = S_T[i][j] + + 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 + + class AnisotropicMineral(Mineral, AnisotropicMaterial): """ A class implementing the anisotropic mineral equation of state described @@ -103,7 +216,7 @@ def __init__( ) raise Exception( "The standard state unit vectors are inconsistent " - "with the volume. Suggest multiplying each vector length" + "with the volume. Suggest multiplying each vector length " f"by {factor}." ) @@ -161,88 +274,34 @@ def set_state(self, pressure, temperature): f = np.log(Vrel) out = self.psi_function(f, self.Pth, self.anisotropic_params) - Psi_Voigt, dPsidf_Voigt, dPsidPth_Voigt = out + Psi_Voigt, dPsidf_Pth_Voigt, dPsidPth_f_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)) + self._PsiI = np.einsum("ijkl, kl", Psi_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 + # Change of variables: (f, Pth) -> Psi(f, T) + aK_T = self.alpha * self.isothermal_bulk_modulus_reuss + out = convert_f_Pth_to_f_T_derivatives( + dPsidf_Pth_Voigt, dPsidPth_f_Voigt, aK_T, self.dPthdf ) + self._dPsidf_T_Voigt, self._dPsiIdf_T, self._dPsiIdT_f = out - # Calculate F, dFdP, dFdT - self._F = expm(PsiI) - + # Calculate F, thermal expansivity and compliance, dFdT if self.orthotropic: - self._S_T_unrotated_Voigt = -dPsidP_Voigt - - self._unrotated_alpha = self.alpha * ( - dPsidfI + dPsidPthI * (self.dPthdf + self.isothermal_bulk_modulus_reuss) + self._unrotated_F = expm(self._PsiI) + self._unrotated_alpha = self._dPsiIdT_f + self.alpha * self._dPsiIdf_T + self._unrotated_S_T_Voigt = ( + self.isothermal_compressibility_reuss * self._dPsidf_T_Voigt ) 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 + out = deformation_gradient_alpha_and_compliance( + self.alpha, + self.isothermal_compressibility_reuss, + self._PsiI, + self._dPsidf_T_Voigt, + self._dPsiIdf_T, + self._dPsiIdT_f, ) - 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 + self._unrotated_F, self._unrotated_alpha, self._unrotated_S_T_Voigt = out @material_property def deformation_gradient_tensor(self): @@ -252,7 +311,7 @@ def deformation_gradient_tensor(self): (i.e. the state at the reference pressure and temperature). :rtype: numpy.array (2D) """ - return self._F + return self._unrotated_F @material_property def unrotated_cell_vectors(self): @@ -366,6 +425,8 @@ def isothermal_bulk_modulus(self): "isothermal_bulk_modulus_reuss?" ) + isothermal_bulk_modulus_reuss = Mineral.isothermal_bulk_modulus + @material_property def isentropic_bulk_modulus(self): """ @@ -382,7 +443,7 @@ def isentropic_bulk_modulus(self): "isentropic_bulk_modulus_reuss?" ) - isothermal_bulk_modulus_reuss = Mineral.isothermal_bulk_modulus + isentropic_bulk_modulus_reuss = Mineral.adiabatic_bulk_modulus @material_property def isothermal_compressibility(self): @@ -422,16 +483,7 @@ def isothermal_bulk_modulus_voigt(self): :returns: The Voigt bound on the isothermal bulk modulus in [Pa]. :rtype: float """ - K = ( - np.sum( - [ - [self.isothermal_stiffness_tensor[i][k] for k in range(3)] - for i in range(3) - ] - ) - / 9.0 - ) - return K + return np.sum(self.isothermal_stiffness_tensor[:3, :3]) / 9.0 @material_property def isothermal_compressibility_reuss(self): @@ -477,10 +529,10 @@ def isothermal_compliance_tensor(self): :rtype: numpy.array (2D) """ if self.orthotropic: - return self._S_T_unrotated_Voigt + return self._unrotated_S_T_Voigt else: R = self.rotation_matrix - S = voigt_notation_to_compliance_tensor(self._S_T_unrotated_Voigt) + S = voigt_notation_to_compliance_tensor(self._unrotated_S_T_Voigt) S_rotated = np.einsum("mi, nj, ok, pl, ijkl->mnop", R, R, R, R, S) return contract_compliances(S_rotated) diff --git a/burnman/classes/material.py b/burnman/classes/material.py index d35e9d21..951cdbfc 100644 --- a/burnman/classes/material.py +++ b/burnman/classes/material.py @@ -202,7 +202,7 @@ def set_state_with_volume( These guesses should preferably bound the correct pressure, but do not need to do so. More importantly, they should not lie outside the valid region of - the equation of state. Defaults to [5.e9, 10.e9]. + the equation of state. Defaults to [0.e9, 10.e9]. :type pressure_guesses: list """ @@ -222,10 +222,13 @@ def _delta_volume(pressure, volume, temperature): try: sol = bracket(_delta_volume, pressure_guesses[0], pressure_guesses[1], args) except ValueError: - raise Exception( - "Cannot find a pressure, perhaps the volume or starting pressures " - "are outside the range of validity for the equation of state?" - ) + try: # Try again with 0 Pa lower bound on the pressure + sol = bracket(_delta_volume, 0.0e9, pressure_guesses[1], args) + except ValueError: + raise Exception( + "Cannot find a pressure, perhaps the volume or starting pressures " + "are outside the range of validity for the equation of state?" + ) pressure = brentq(_delta_volume, sol[0], sol[1], args=args) self.set_state(pressure, temperature) @@ -613,6 +616,16 @@ def molar_heat_capacity_p(self): "need to implement molar_heat_capacity_p() in derived class!" ) + @material_property + def isentropic_thermal_gradient(self): + """ + :returns: dTdP, the change in temperature with pressure at constant entropy [Pa/K] + :rtype: float + """ + raise NotImplementedError( + "need to implement isentropic_thermal_gradient() in derived class!" + ) + # # Aliased properties @property diff --git a/burnman/classes/mineral.py b/burnman/classes/mineral.py index eedb52e4..7ccf0161 100644 --- a/burnman/classes/mineral.py +++ b/burnman/classes/mineral.py @@ -1,4 +1,5 @@ -# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for the Earth and Planetary Sciences +# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit +# for the Earth and Planetary Sciences # Copyright (C) 2012 - 2017 by the BurnMan team, released under the GNU # GPL v2 or later. @@ -197,6 +198,8 @@ def isothermal_bulk_modulus(self): - self._property_modifiers["d2GdP2"] ) + isothermal_bulk_modulus_reuss = isothermal_bulk_modulus + @material_property @copy_documentation(Material.molar_heat_capacity_p) def molar_heat_capacity_p(self): @@ -379,3 +382,10 @@ def molar_heat_capacity_v(self): * self.thermal_expansivity * self.isothermal_bulk_modulus ) + + @material_property + @copy_documentation(Material.isentropic_thermal_gradient) + def isentropic_thermal_gradient(self): + return ( + self.molar_volume * self.temperature * self.thermal_expansivity + ) / self.molar_heat_capacity_p diff --git a/burnman/tools/eos.py b/burnman/tools/eos.py index ad62fe59..4f58a21c 100644 --- a/burnman/tools/eos.py +++ b/burnman/tools/eos.py @@ -191,7 +191,7 @@ def equilibration_function(mineral): def check_anisotropic_eos_consistency( - m, P=1.0e9, T=2000.0, tol=1.0e-4, verbose=False, equilibration_function=None + m, P=1.0e9, T=300.0, tol=1.0e-4, verbose=False, equilibration_function=None ): """ Checks that numerical derivatives of the Gibbs energy of an anisotropic mineral