Skip to content

Commit

Permalink
updated test
Browse files Browse the repository at this point in the history
  • Loading branch information
@bobmyhill
bobmyhill committed Feb 11, 2024
1 parent bb81cbe commit 7a650cc
Show file tree
Hide file tree
Showing 3 changed files with 78 additions and 31 deletions.
35 changes: 30 additions & 5 deletions burnman/tools/eos.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down
1 change: 1 addition & 0 deletions test.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
73 changes: 47 additions & 26 deletions tests/test_anisotropicmineral.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand Down Expand Up @@ -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__":
Expand Down

0 comments on commit 7a650cc

Please sign in to comment.