From 7a650cc5d4a60050f2ab2ee8f042a4a5def682e5 Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Sat, 10 Feb 2024 18:41:58 +0000 Subject: [PATCH] updated test --- burnman/tools/eos.py | 35 ++++++++++++--- test.sh | 1 + tests/test_anisotropicmineral.py | 73 ++++++++++++++++++++------------ 3 files changed, 78 insertions(+), 31 deletions(-) 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/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 e02cc3e80..9048a9753 100644 --- a/tests/test_anisotropicmineral.py +++ b/tests/test_anisotropicmineral.py @@ -8,28 +8,44 @@ 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], ] ) - m = AnisotropicMineral(fo, cell_parameters, constants, orthotropic=True) + m = AnisotropicMineral(fo, cell_parameters, constants) return m @@ -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() - self.assertTrue(check_anisotropic_eos_consistency(m, verbose=True)) + 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__":