Merge pull request #566 from bobmyhill/non_orthotropic_anisotropy

Non orthotropic anisotropy

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")
 

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.

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.

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,

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

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

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

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]