# 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
# GPL v2 or later.
import numpy as np
from scipy.linalg import expm, logm
from numpy.linalg import cond
from .mineral import Mineral
from .material import Material, material_property
from .anisotropy import AnisotropicMaterial
from ..utils.misc import copy_documentation
from ..utils.unitcell import cell_parameters_to_vectors
from ..utils.unitcell import cell_vectors_to_parameters
from ..utils.anisotropy import (
voigt_notation_to_compliance_tensor,
voigt_notation_to_stiffness_tensor,
contract_compliances,
)
[docs]class AnisotropicMineral(Mineral, AnisotropicMaterial):
"""
A class implementing the anisotropic mineral equation of state described
in :cite:`Myhill2022`.
This class is derived from both Mineral and AnisotropicMaterial,
and inherits most of the methods from these classes.
Instantiation of an AnisotropicMineral takes three required arguments;
a reference Mineral (i.e. a standard isotropic mineral which provides
volume as a function of pressure and temperature), cell_parameters,
which give the lengths of the molar cell vectors and the angles between
them (see :func:`~burnman.utils.unitcell.cell_parameters_to_vectors`),
and an anisotropic parameters object, which should be either a
4D array of anisotropic parameters or a dictionary of parameters which
describe the anisotropic behaviour of the mineral.
For a description of the physical meaning of the parameters in the
4D array, please refer to the code
or to the original paper.
If the user chooses to define their parameters as a dictionary,
they must also provide a function to the psi_function argument
that describes how to compute the tensors Psi, dPsidf and dPsidPth
(in Voigt form). The function arguments should be f, Pth and params,
in that order. The output variables Psi, dPsidf and dPsidth
must be returned in that order in a tuple. The user should
also explicitly state whether the material is orthotropic or not
by supplying a boolean to the orthotropic argument.
States of the mineral can only be queried after setting the
pressure and temperature using set_state().
This class is available as ``burnman.AnisotropicMineral``.
All the material parameters are expected to be in plain SI units. This
means that the elastic moduli should be in Pascals and NOT Gigapascals.
Additionally, the cell parameters should be in m/(mol formula unit)
and not in unit cell lengths. To convert unit cell lengths given in
Angstrom to molar cell parameters you should multiply by 10^(-10) *
(N_a / Z)^1/3, where N_a is Avogadro's number
and Z is the number of formula units per unit cell.
You can look up Z in many places, including www.mindat.org.
Finally, it is assumed that the unit cell of the anisotropic material
is aligned in a particular way relative to the coordinate axes
(the anisotropic_parameters are defined relative to the coordinate axes).
The crystallographic a-axis is assumed to be parallel to the first
spatial coordinate axis, and the crystallographic b-axis is assumed to
be perpendicular to the third spatial coordinate axis.
"""
def __init__(
self,
isotropic_mineral,
cell_parameters,
anisotropic_parameters,
psi_function=None,
orthotropic=None,
):
if psi_function is None:
self.check_standard_parameters(anisotropic_parameters)
self.anisotropic_params = {"c": anisotropic_parameters}
self.psi_function = self.standard_psi_function
else:
if not isinstance(orthotropic, bool):
raise Exception(
"If the Psi function is provided, "
"you must specify whether your material is "
"orthotropic as a boolean."
)
self.orthotropic = orthotropic
self.anisotropic_params = anisotropic_parameters
self.psi_function = psi_function
self.cell_vectors_0 = cell_parameters_to_vectors(cell_parameters)
if (
np.abs(np.linalg.det(self.cell_vectors_0) - isotropic_mineral.params["V_0"])
> np.finfo(float).eps
):
factor = np.cbrt(
isotropic_mineral.params["V_0"] / np.linalg.det(self.cell_vectors_0)
)
raise Exception(
"The standard state unit vectors are inconsistent "
"with the volume. Suggest multiplying each "
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"]
Mineral.__init__(
self, isotropic_mineral.params, isotropic_mineral.property_modifiers
)
[docs] def standard_psi_function(self, f, Pth, params):
# Compute Psi, dPsidPth, dPsidf, needed by most anisotropic properties
c = params["c"]
ns = np.arange(c.shape[-1])
x = c[:, :, 0, :] + c[:, :, 1, :] * f
dPsidf = c[:, :, 1, :]
for i in list(range(2, c.shape[2])):
# non-intuitively, the += operator doesn't simply add in-place,
# so here we overwrite the arrays with new ones
x = x + c[:, :, i, :] * np.power(f, float(i)) / float(i)
dPsidf = dPsidf + c[:, :, i, :] * np.power(f, float(i) - 1.0)
Psi = np.einsum("ikn, n->ik", x, np.power(Pth, ns))
dPsidPth = np.einsum(
"ikn, n->ik", x[:, :, 1:], ns[1:] * np.power(Pth, ns[1:] - 1)
)
dPsidf = np.einsum("ikn, n->ik", dPsidf, np.power(Pth, ns))
return (Psi, dPsidf, dPsidPth)
@copy_documentation(Material.set_state)
def set_state(self, pressure, temperature):
# 1) Compute dPthdf|T
# relatively large dP needed for accurate estimate of dPthdf
self.isotropic_mineral.set_state(pressure, temperature)
V2 = self.isotropic_mineral.V
KT2 = self.isotropic_mineral.K_T
self.isotropic_mineral.set_state_with_volume(V2, self.params["T_0"])
P1 = self.isotropic_mineral.pressure
KT1 = self.isotropic_mineral.K_T
self.dPthdf = KT1 - KT2
self.Pth = pressure - P1
self.isotropic_mineral.set_state(pressure, temperature)
Mineral.set_state(self, pressure, temperature)
# 2) Compute other properties needed for anisotropic equation of state
V = self.V
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)
)
@material_property
def deformation_gradient_tensor(self):
"""
:returns: The deformation gradient tensor describing the deformation of the
mineral from its undeformed state
(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
@material_property
def unrotated_cell_vectors(self):
"""
:returns: The vectors of the cell [m] constructed from one mole
of formula units after deformation of the mineral from its
undeformed state (i.e. the state at the reference
pressure and temperature). See the documentation for the function
:func:`~burnman.utils.unitcell.cell_parameters_to_vectors`
for the assumed relationships between the cell vectors and
spatial coordinate axes.
:rtype: numpy.array (2D)
"""
return self.deformation_gradient_tensor
@material_property
def deformed_coordinate_frame(self):
"""
:returns: The orientations of the three spatial coordinate axes
after deformation of the mineral [m]. For orthotropic minerals,
this is equal to the identity matrix, as hydrostatic stresses only
induce rotations in monoclinic and triclinic crystals.
:rtype: numpy.array (2D)
"""
if self.orthotropic:
return np.eye(3)
else:
M = 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[1] = np.cross(Q[2], Q[0])
return Q
@material_property
def rotation_matrix(self):
"""
:returns: The matrix required to rotate the properties of the deformed
mineral into the deformed coordinate frame. For orthotropic
minerals, this is equal to the identity matrix.
:rtype: numpy.array (2D)
"""
return self.deformed_coordinate_frame.T
@material_property
def cell_vectors(self):
"""
:returns: The vectors of the cell constructed from one mole
of formula units [m]. See the documentation for the function
:func:`~burnman.utils.unitcell.cell_parameters_to_vectors`
for the assumed relationships between the cell vectors and
spatial coordinate axes.
:rtype: numpy.array (2D)
"""
if self.orthotropic:
return self.unrotated_cell_vectors
else:
return np.einsum(
"ij, jk->ik", self.unrotated_cell_vectors, self.rotation_matrix
)
@material_property
def cell_parameters(self):
"""
:returns: The molar cell parameters of the mineral, given in standard form:
[:math:`a`, :math:`b`, :math:`c`,
:math:`\\alpha`, :math:`\\beta`, :math:`\\gamma`],
where the first three floats are the lengths of the vectors in [m]
defining the cell constructed from one mole of formula units.
The last three floats are angles between vectors
(given in radians). See the documentation for the function
:func:`~burnman.utils.unitcell.cell_parameters_to_vectors`
for the assumed relationships between the cell vectors and
spatial coordinate axes.
:rtype: numpy.array (1D)
"""
return cell_vectors_to_parameters(self.cell_vectors)
@material_property
def shear_modulus(self):
"""
Anisotropic minerals do not (in general) have a single shear modulus.
This function returns a NotImplementedError. Users should instead
consider directly querying the elements in the
isothermal_stiffness_tensor or isentropic_stiffness_tensor.
"""
raise NotImplementedError(
"Anisotropic minerals do not have a shear "
"modulus property. Query "
"the isentropic or isothermal stiffness "
"tensors directory, or use"
"isentropic_shear_modulus_reuss or "
"isentropic_shear_modulus_voigt."
)
@material_property
def isothermal_bulk_modulus(self):
"""
Anisotropic minerals do not have a single isothermal bulk modulus.
This function returns a NotImplementedError. Users should instead
consider either using isothermal_bulk_modulus_reuss,
isothermal_bulk_modulus_voigt,
or directly querying the elements in the isothermal_stiffness_tensor.
"""
raise NotImplementedError(
"isothermal_bulk_modulus is not "
"sufficiently explicit for an "
"anisotropic mineral. Did you mean "
"isothermal_bulk_modulus_reuss?"
)
@material_property
def isentropic_bulk_modulus(self):
"""
Anisotropic minerals do not have a single isentropic bulk modulus.
This function returns a NotImplementedError. Users should instead
consider either using isentropic_bulk_modulus_reuss,
isentropic_bulk_modulus_voigt (both derived from AnisotropicMineral),
or directly querying the elements in the isentropic_stiffness_tensor.
"""
raise NotImplementedError(
"isentropic_bulk_modulus is not "
"sufficiently explicit for an "
"anisotropic mineral. Did you mean "
"isentropic_bulk_modulus_reuss?"
)
isothermal_bulk_modulus_reuss = Mineral.isothermal_bulk_modulus
@material_property
def isothermal_compressibility(self):
"""
Anisotropic minerals do not have a single isentropic compressibility.
This function returns a NotImplementedError. Users should instead
consider either using isothermal_compressibility_reuss,
isothermal_compressibility_voigt (both derived from AnisotropicMineral),
or directly querying the elements in the isothermal_compliance_tensor.
"""
raise NotImplementedError(
"isothermal_compressibility is not "
"sufficiently explicit for an "
"anisotropic mineral. Did you mean "
"isothermal_compressibility_reuss?"
)
@material_property
def isentropic_compressibility(self):
"""
Anisotropic minerals do not have a single isentropic compressibility.
This function returns a NotImplementedError. Users should instead
consider either using isentropic_compressibility_reuss,
isentropic_compressibility_voigt (both derived from AnisotropicMineral),
or directly querying the elements in the isentropic_compliance_tensor.
"""
raise NotImplementedError(
"isentropic_compressibility is not "
"sufficiently explicit for an "
"anisotropic mineral. Did you mean "
"isentropic_compressibility_reuss?"
)
@material_property
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
@material_property
def isothermal_compressibility_reuss(self):
"""
:returns: The Reuss bound on the isothermal compressibility in [1/Pa].
:rtype: float
"""
return 1.0 / self.isothermal_bulk_modulus_reuss
beta_T = isothermal_compressibility_reuss
@material_property
def isothermal_compressibility_voigt(self):
"""
:returns: The Voigt bound on the isothermal compressibility in [1/Pa].
:rtype: float
"""
return 1.0 / self.isothermal_bulk_modulus_voigt
@material_property
def isentropic_compressibility_reuss(self):
"""
:returns: The Reuss bound on the isentropic compressibility in [1/Pa].
:rtype: float
"""
return 1.0 / self.isentropic_bulk_modulus_reuss
beta_S = isentropic_compressibility_reuss
@material_property
def isentropic_compressibility_voigt(self):
"""
:returns: The Voigt bound on the isentropic compressibility in [1/Pa].
:rtype: float
"""
return 1.0 / self.isentropic_bulk_modulus_voigt
@material_property
def isothermal_compliance_tensor(self):
"""
:returns: The isothermal compliance tensor [1/Pa]
in Voigt form (:math:`\\mathbb{S}_{\\text{T} pq}`).
:rtype: numpy.array (2D)
"""
S_T = -self.dPsidP_Voigt
if self.orthotropic:
return S_T
else:
R = self.rotation_matrix
S = voigt_notation_to_compliance_tensor(S_T)
S_rotated = np.einsum("mi, nj, ok, pl, ijkl->mnop", R, R, R, R, S)
return contract_compliances(S_rotated)
@material_property
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
else:
R = self.rotation_matrix
return np.einsum("mi, nj, ij->mn", R, R, alpha)
# Derived properties start here
@material_property
def isothermal_stiffness_tensor(self):
"""
:returns: The isothermal stiffness tensor [Pa]
in Voigt form (:math:`\\mathbb{C}_{\\text{T} pq}`).
:rtype: numpy.array (2D)
"""
return np.linalg.inv(self.isothermal_compliance_tensor)
@material_property
def full_isothermal_compliance_tensor(self):
"""
:returns: The isothermal compliance tensor [1/Pa]
in standard form (:math:`\\mathbb{S}_{\\text{T} ijkl}`).
:rtype: numpy.array (4D)
"""
S_Voigt = self.isothermal_compliance_tensor
return voigt_notation_to_compliance_tensor(S_Voigt)
@material_property
def full_isothermal_stiffness_tensor(self):
"""
:returns: The isothermal stiffness tensor [Pa]
in standard form (:math:`\\mathbb{C}_{\\text{T} ijkl}`).
:rtype: numpy.array (4D)
"""
CT = self.isothermal_stiffness_tensor
return voigt_notation_to_stiffness_tensor(CT)
@material_property
def full_isentropic_compliance_tensor(self):
"""
:returns: The isentropic compliance tensor [1/Pa]
in standard form (:math:`\\mathbb{S}_{\\text{N} ijkl}`).
:rtype: numpy.array (4D)
"""
return (
self.full_isothermal_compliance_tensor
- np.einsum(
"ij, kl->ijkl",
self.thermal_expansivity_tensor,
self.thermal_expansivity_tensor,
)
* self.V
* self.temperature
/ self.C_p
)
@material_property
def isentropic_compliance_tensor(self):
"""
:returns: The isentropic compliance tensor [1/Pa]
in Voigt form (:math:`\\mathbb{S}_{\\text{N} pq}`).
:rtype: numpy.array (2D)
"""
S_full = self.full_isentropic_compliance_tensor
return contract_compliances(S_full)
@material_property
def isentropic_stiffness_tensor(self):
"""
:returns: The isentropic stiffness tensor [Pa]
in Voigt form (:math:`\\mathbb{C}_{\\text{N} pq}`).
:rtype: numpy.array (2D)
"""
return np.linalg.inv(self.isentropic_compliance_tensor)
@material_property
def full_isentropic_stiffness_tensor(self):
"""
:returns: The isentropic stiffness tensor [Pa]
in standard form (:math:`\\mathbb{C}_{\\text{N} ijkl}`).
:rtype: numpy.array (4D)
"""
C_Voigt = self.isentropic_stiffness_tensor
return voigt_notation_to_stiffness_tensor(C_Voigt)
@material_property
def grueneisen_tensor(self):
"""
:returns: The grueneisen tensor [unitless].
This is defined by :cite:`BarronMunn1967` as
:math:`\\mathbb{C}_{\\text{N} ijkl} \\alpha_{kl} V/C_{P}`.
:rtype: numpy.array (2D)
"""
return (
np.einsum(
"ijkl, kl->ij",
self.full_isentropic_stiffness_tensor,
self.thermal_expansivity_tensor,
)
* self.molar_volume
/ self.molar_heat_capacity_p
)
@material_property
def grueneisen_parameter(self):
"""
:returns: The scalar grueneisen parameter [unitless].
:rtype: float
"""
return (
self.thermal_expansivity
* self.V
/ (self.isentropic_compressibility_reuss * self.molar_heat_capacity_p)
)
@material_property
def isothermal_compressibility_tensor(self):
"""
:returns: The isothermal compressibility tensor [1/Pa].
:rtype: numpy.array (2D)
"""
return np.einsum(
"ijkl, kl->ij", self.full_isothermal_compliance_tensor, np.eye(3)
)
@material_property
def isentropic_compressibility_tensor(self):
"""
:returns: The isentropic compressibility tensor [1/Pa].
:rtype: numpy.array (2D)
"""
return np.einsum(
"ijkl, kl->ij", self.full_isentropic_compliance_tensor, np.eye(3)
)
@material_property
def thermal_stress_tensor(self):
"""
:returns: The change in stress with temperature at constant strain [Pa/K].
:rtype: numpy.array (2D)
"""
pi = -np.einsum(
"ijkl, kl",
self.full_isothermal_stiffness_tensor,
self.thermal_expansivity_tensor,
)
return pi
@material_property
def molar_isometric_heat_capacity(self):
"""
:returns: The molar heat capacity at constant strain [J/K/mol].
:rtype: float
"""
alpha = self.thermal_expansivity_tensor
pi = self.thermal_stress_tensor
C_isometric = (
self.molar_heat_capacity_p
+ self.V * self.temperature * np.einsum("ij, ij", alpha, pi)
)
return C_isometric
[docs] def check_standard_parameters(self, anisotropic_parameters):
if not np.all(anisotropic_parameters[:, :, 0, 0] == 0):
raise Exception(
"anisotropic_parameters_pqmn should be set to " "zero for all m = n = 0"
)
sum_ijij_block = np.sum(anisotropic_parameters[:3, :3, :, :], axis=(0, 1))
if np.abs(sum_ijij_block[1, 0] - 1.0) > 1.0e-5:
raise Exception(
"The sum of the upper 3x3 pq-block of "
"anisotropic_parameters_pqmn must equal "
"1 for m=1, n=0 for consistency with the volume. "
f"Value is {sum_ijij_block[1, 0]}"
)
for m in range(2, len(sum_ijij_block)):
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"
f"m={m}, n=0 for consistency with the volume. "
f"Value is {sum_ijij_block[m, 0]}"
)
for m in range(len(sum_ijij_block)):
for n in range(1, len(sum_ijij_block[0])):
if np.abs(sum_ijij_block[m, n]) > 1.0e-10:
raise Exception(
"The sum of the upper 3x3 pq-block of "
"anisotropic_parameters_pqmn must equal "
f"0 for m={m}, n={n} for "
"consistency with the volume. "
f"Value is {sum_ijij_block[m, n]}"
)
if cond(anisotropic_parameters[:, :, 1, 0]) > 1 / np.finfo(float).eps:
raise Exception("anisotropic_parameters[:, :, 1, 0] is singular")
sum_lower_left_block = np.sum(anisotropic_parameters[3:, :3, :, :], axis=1)
self.orthotropic = True
for i, s in enumerate(sum_lower_left_block):
if not np.all(np.abs(s) < 1.0e-10):
self.orthotropic = False