diff --git a/burnman/__init__.py b/burnman/__init__.py index f8e00f3b..54b606d0 100644 --- a/burnman/__init__.py +++ b/burnman/__init__.py @@ -228,7 +228,7 @@ from .classes.perplex import PerplexMaterial from .classes.mineral import Mineral from .classes.combinedmineral import CombinedMineral -from .classes.solution import Solution, SolidSolution +from .classes.solution import Solution, SolidSolution, RelaxedSolution from .classes.elasticsolutionmodel import ElasticSolutionModel from .classes.elasticsolution import ElasticSolution, ElasticSolidSolution from .classes.composite import Composite @@ -237,7 +237,7 @@ from .classes.anisotropicmineral import AnisotropicMineral from .classes.anisotropicmineral import cell_parameters_to_vectors from .classes.anisotropicmineral import cell_vectors_to_parameters -from .classes.anisotropicsolution import AnisotropicSolution +from .classes.anisotropicsolution import AnisotropicSolution, RelaxedAnisotropicSolution from .classes.mineral_helpers import HelperLowHighPressureRockTransition from .classes.mineral_helpers import HelperSpinTransition from .classes.mineral_helpers import HelperRockSwitcher diff --git a/burnman/classes/anisotropicsolution.py b/burnman/classes/anisotropicsolution.py index 6d8ca465..9470b651 100644 --- a/burnman/classes/anisotropicsolution.py +++ b/burnman/classes/anisotropicsolution.py @@ -5,6 +5,7 @@ import numpy as np import copy from scipy.linalg import expm, logm +from scipy.optimize import minimize from .solution import Solution from .anisotropicmineral import ( AnisotropicMineral, @@ -16,10 +17,8 @@ contract_stresses, expand_stresses, voigt_notation_to_compliance_tensor, - voigt_notation_to_stiffness_tensor, - contract_compliances, - contract_stiffnesses, ) +from ..utils.misc import copy_documentation class AnisotropicSolution(Solution, AnisotropicMineral): @@ -89,6 +88,19 @@ def set_state(self, pressure, temperature): raise Exception("To use this EoS, you must first set the composition") ss.set_state(pressure, temperature) + + try: + self.compute_base_properties() + except AttributeError: + raise Exception("You must use set_composition() before set_state().") + + Material.set_state(self, pressure, temperature) + + def compute_base_properties(self): + pressure = self._scalar_solution.pressure + temperature = self._scalar_solution.temperature + + ss = self._scalar_solution V = ss.molar_volume KT_at_T = ss.isothermal_bulk_modulus_reuss f = np.log(V) @@ -153,10 +165,562 @@ def set_state(self, pressure, temperature): def set_composition(self, molar_fractions): self._scalar_solution.set_composition(molar_fractions) - Solution.set_composition(self, molar_fractions) self.cell_vectors_0 = expm( np.einsum("p, pij->ij", molar_fractions, self._logm_M_T_0_mbr) ) # Calculate all other required properties - if self.pressure is not None: - self.set_state(self.pressure, self.temperature) + Solution.set_composition(self, molar_fractions) + if self._scalar_solution.pressure is not None: + self.compute_base_properties() + + @cached_property + def ones(self): + return np.ones(self.n_endmembers) + + @cached_property + def eye(self): + return np.eye(self.n_endmembers) + + @material_property + def _dM0dp_fixed_VT(self): + """ + Gradient in the standard state cell tensor + with respect to molar proportions at constant + volume and temperature under hydrostatic conditions. + """ + lnM0_ones = np.einsum("ij, k->kij", logm(self.cell_vectors_0.T), self.ones) + dp = 1.0e-5 + + dlnM0 = np.einsum("ijk->kij", self._logm_M0_mbr) * dp / 2.0 + logmM0_0 = lnM0_ones - dlnM0 + logmM0_1 = lnM0_ones + dlnM0 + return np.einsum("kij->ijk", (expm(logmM0_1) - expm(logmM0_0)) / dp) + + @material_property + def _unrotated_dFdp_fixed_VT(self): + """ + Gradient in the unrotated deformation gradient tensor + with respect to molar proportions at constant + volume and temperature under hydrostatic conditions. + """ + dp = 1.0e-5 + PsiI_ones = np.einsum("ij, k->kij", self._PsiI, self.ones) + dPsiI = np.einsum("ijk->kij", self._PsiI_mbr + self._dPsiIdp_xs) * dp / 2.0 + PsiI_0 = PsiI_ones - dPsiI + PsiI_1 = PsiI_ones + dPsiI + + return np.einsum("kij->ijk", (expm(PsiI_1) - expm(PsiI_0)) / dp) + + @material_property + def _unrotated_dMdn_fixed_VT(self): + """ + Gradient in unrotated cell tensor with respect to composition + at fixed volume and temperature under hydrostatic conditions. + """ + F = self._unrotated_F + M0 = self.cell_vectors_0.T + fr = self.molar_fractions + n_dpdn = self.eye - np.einsum("l, m->lm", self.ones, fr) + n = np.sum(self.molar_fractions) + n23 = np.power(n, 2.0 / 3.0) + + dM0dn = ( + np.einsum("kj, l->kjl", M0, self.ones / 3.0) + + np.einsum("kjm,lm->kjl", self._dM0dp_fixed_VT, n_dpdn) + ) / n23 + n_dVmoldn = -self.molar_volume * self.ones + + dFdV = self._unrotated_dFdf_T / self.molar_volume + dFdn = ( + np.einsum("ij, l->ijl", dFdV, n_dVmoldn) + + np.einsum("ikm,lm->ikl", self._unrotated_dFdp_fixed_VT, n_dpdn) + ) / n + + dFdnM0 = np.einsum("ikl,kj->ijl", dFdn, M0) + FdM0dn = np.einsum("ik,kjl->ijl", F, dM0dn) + return dFdnM0 + FdM0dn + + @material_property + def _dMdn_fixed_VT(self): + """ + Gradient in cell tensor with respect to composition + at fixed volume and temperature under hydrostatic conditions. + """ + R = self.rotation_matrix + return np.einsum("mi, nj, ijk->mnk", R, R, self._unrotated_dMdn_fixed_VT) + + @material_property + def depsdn_fixed_VT(self): + """ + Gradient in strain with respect to composition + at fixed volume and temperature under hydrostatic conditions. + """ + invM = np.linalg.inv(self.cell_vectors.T) + Ln = np.einsum("ijm,kj->ikm", self._dMdn_fixed_VT, invM) + LnT = np.einsum("ikm->kim", Ln) + return 0.5 * (Ln + LnT) + + +depsdeps = np.einsum("pm, qn->pqmn", np.eye(3), np.eye(3)) + + +class RelaxedAnisotropicSolution(AnisotropicSolution): + """ + A class implementing the relaxed anisotropic solution + model described in :cite:`Myhill2024b`. + This class is derived from AnisotropicSolution, + and inherits most of the methods from that class. + + Instantiation of a RelaxedAnisotropicSolution involves passing + an AnisotropicSolution, plus a set of vectors that represent rapid + deformation modes. For example, a solution of MgO, FeHSO and FeLSO + (high and low spin wuestite) can rapidly change proportion of + high spin and low spin iron, and so a single vector should be passed: + np.array([[0., -1., 1.]]) or some multiple thereof. + + States of the mineral can only be queried after setting the + pressure and temperature using set_state() and the composition using + set_composition(). + + This class is available as ``burnman.RelaxedAnisotropicSolution``. + """ + + def __init__( + self, + anisotropic_solution, + relaxation_vectors, + unrelaxed_vectors, + ): + # Make an attribute with the unrelaxed solution + self.unrelaxed = anisotropic_solution + + # The relaxation vectors + self.n_relaxation_vectors = len(relaxation_vectors) + self.n_unrelaxed_vectors = len(unrelaxed_vectors) + + self.dndq = np.array(relaxation_vectors).T + assert len(self.dndq.shape) == 2 + assert len(self.dndq) == self.unrelaxed.n_endmembers + + self.dndx = np.array(unrelaxed_vectors).T + assert len(self.dndx.shape) == 2 + assert len(self.dndx) == self.unrelaxed.n_endmembers + assert ( + len(unrelaxed_vectors) + len(relaxation_vectors) + == self.unrelaxed.n_endmembers + ) + + self.q_initial = np.zeros(len(relaxation_vectors)) + 0.001 + + try: + molar_fractions = anisotropic_solution.molar_fractions + except AttributeError: + molar_fractions = None + + # Give the relaxed solution the same base properties as the + # unrelaxed solution + AnisotropicSolution.__init__( + self, + name=anisotropic_solution.name, + solution_model=anisotropic_solution.solution_model, + psi_excess_function=anisotropic_solution.psi_excess_function, + anisotropic_parameters=anisotropic_solution.anisotropic_params, + molar_fractions=molar_fractions, + ) + + def set_state(self, pressure, temperature, relaxed=True): + """ + Sets the state of the solution. Also relaxes the + structure parameters if set_composition() has already + been used and if the relaxed argument has been set to + True. + + :param pressure: The pressure of the solution [Pa] + :type pressure: float + :param temperature: The temperature of the solution [K] + :type temperature: float + :param relaxed: Whether to minimize the Gibbs energy + of the material by changing the values of the structure + parameters. Defaults to True. + :type relaxed: bool, optional + """ + self.unrelaxed.set_state(pressure, temperature) + + if hasattr(self.unrelaxed, "molar_fractions") and relaxed: + self._relax_at_PTX() + AnisotropicSolution.set_state(self, pressure, temperature) + + def set_composition(self, molar_fractions, q_initial=None, relaxed=True): + """ + Sets the composition of the model. Also relaxes the + structure parameters if set_state() has already + been used and if the relaxed argument has been set to + True. + + :param molar_fractions: Molar fractions of the + independent endmembers corresponding to the + unrelaxed vectors specified during initialisation. + :type molar_fractions: 1D numpy array + :param q_initial: Initial values of the structure parameters. + Defaults to None, in which case the preexisting + initial values are used + (first set to 0.001 during initialisation). + :type q_initial: 1D numpy array, optional + :param relaxed: Whether to minimize the Gibbs energy + of the material by changing the values of the structure + parameters. Defaults to True. + :type relaxed: bool, optional + """ + self.unrelaxed_vectors = molar_fractions + if q_initial is not None: + self.q_initial = np.array(q_initial) + n = np.einsum("ij, j", self.dndq, self.q_initial) + np.einsum( + "ij, j", self.dndx, molar_fractions + ) + + self.unrelaxed.set_composition(n) + + if self.unrelaxed.pressure is not None and relaxed: + self._relax_at_PTX() + AnisotropicSolution.set_composition(self, n) + self.unrelaxed._scalar_solution.set_composition(self.molar_fractions) + self.unrelaxed.set_composition(self.molar_fractions) + + def _relax_at_PTX(self): + """ + Minimizes the Gibbs energy at constant pressure and + temperature by changing the structural parameters. + + Run during set_state() and set_composition(), as long as both + state and composition have already been set. This function + should not generally be needed by the user. + """ + n0 = copy.copy(self.unrelaxed._scalar_solution.molar_fractions) + + def G_func(dq): + n = n0 + np.einsum("ij, j", self.dndq, dq) + self.unrelaxed._scalar_solution.set_composition(n) + return self.unrelaxed._scalar_solution.molar_gibbs + + sol = minimize(G_func, np.zeros(len(self.dndq[0])), method="Nelder-Mead") + assert sol.success + n = n0 + np.einsum("ij, j", self.dndq, sol.x) + self.unrelaxed.set_composition(n) + AnisotropicSolution.set_composition(self, n) + + def set_state_with_volume( + self, volume, temperature, pressure_guesses=[0.0e9, 10.0e9] + ): + """ + This function acts similarly to set_state, but takes volume and + temperature as input to find the pressure. In order to ensure + self-consistency, this function does not use any pressure functions + from the material classes, but instead finds the pressure using the + brentq root-finding and Nelder-Mead minimization methods. + + Composition should have been set before this function is used. + + :param volume: The desired molar volume of the mineral [m^3]. + :type volume: float + + :param temperature: The desired temperature of the mineral [K]. + :type temperature: float + + :param pressure_guesses: A list of floats denoting the initial + low and high guesses for bracketing of the pressure [Pa]. + 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 [0.e9, 10.e9]. + :type pressure_guesses: list + """ + ss = self.unrelaxed._scalar_solution + n0 = copy.copy(ss.molar_fractions) + + # Initial set_state without changing structural parameters + # to estimate pressure + ss.set_state_with_volume(volume, temperature, pressure_guesses) + # Store in a mutable so that it can be updated in the loop + pressure = [ss.pressure] + + def F_func(dq): + n = n0 + np.einsum("ij, j", self.dndq, dq) + ss.set_composition(n) + ss.set_state_with_volume( + volume, temperature, [pressure[0] - 1.0e9, pressure[0] + 1.0e9] + ) + pressure[0] = ss.pressure + return ss.molar_helmholtz + + sol = minimize(F_func, np.zeros(len(self.dndq[0])), method="Nelder-Mead") + assert sol.success + + # Run one more time with root + F_func(sol.x) + self.unrelaxed.set_composition(ss.molar_fractions) + AnisotropicSolution.set_composition(self, ss.molar_fractions) + self.set_state(pressure[0], temperature) + + @material_property + def _depsdq_fixed_VT(self): + """ + Gradient in strain with respect to structural state + at fixed volume and temperature under hydrostatic conditions. + """ + return np.einsum("mnq, qu->mnu", self.depsdn_fixed_VT, self.dndq) + + @material_property + def _depsdT_fixed_volume(self): + """ + Gradient in strain with respect to temperature + at fixed volume and structural state under hydrostatic conditions. + """ + beta_T = self.unrelaxed.isothermal_compressibility_tensor + beta_TR = self.unrelaxed.isothermal_compressibility_reuss + alpha = self.unrelaxed.thermal_expansivity_tensor + alpha_V = self.unrelaxed.alpha + + return (alpha / alpha_V - beta_T / beta_TR) * alpha_V + + @material_property + def _dbarepsdeps(self): + """ + Gradient in non-hydrostatic strain with + respect to total strain at fixed temperature + and structural state under hydrostatic conditions. + """ + beta_T = self.unrelaxed.isothermal_compressibility_tensor + beta_TR = self.unrelaxed.isothermal_compressibility_reuss + return depsdeps - np.einsum("pq, mn->pqmn", beta_T / beta_TR, np.eye(3)) + + @material_property + def _dSdq(self): + """ + Gradient in strain with respect to structural state + at constant volume and temperature under hydrostatic + conditions. + + This function converts the partial entropies + from Gibbs natural variables (P, T) + to Helmholtz natural variables (V, T). + """ + dSdn = ( + self._scalar_solution.partial_entropies + - self._scalar_solution.alpha + * self._scalar_solution.isothermal_bulk_modulus_reuss + * self._scalar_solution.partial_volumes + ) + return np.einsum("i, ij->j", dSdn, self.dndq) + + @material_property + def _dPdq(self): + """ + Pressure gradient with respect to structural state + calculated at constant volume and temperature + under hydrostatic conditions. + + This function converts from Gibbs natural variables + (partial volumes as a function of P and T) + to Helmholtz natural variables + (partial pressures as a function of V and T). + """ + dPdn = ( + self._scalar_solution.isothermal_bulk_modulus_reuss + / self._scalar_solution.V + * self._scalar_solution.partial_volumes + ) + return np.einsum("i, ij->j", dPdn, self.dndq) + + @material_property + def _d2Fdqdq_fixed_VT(self): + """ + Helmholtz structural hessian + calculated at constant volume and temperature + under hydrostatic conditions. + + This function converts from Gibbs natural variables + (the Gibbs compositional hessian at constant P and T) + to Helmholtz natural variables + (the Helmholtz compositional hessian + at constant V and T). + """ + d2Fdndn = ( + self._scalar_solution.gibbs_hessian + + np.einsum( + "i, j->ij", + self._scalar_solution.partial_volumes, + self._scalar_solution.partial_volumes, + ) + * self._scalar_solution.isothermal_bulk_modulus_reuss + / self._scalar_solution.V + ) + return np.einsum("ij, ik, jl->kl", d2Fdndn, self.dndq, self.dndq) + + @material_property + def _d2Fdqdq_fixed_epsT(self): + """ + Second structure parameter derivative of the Helmholtz energy + at constant strain and temperature. + """ + C_T = self.unrelaxed.full_isothermal_stiffness_tensor + return self._d2Fdqdq_fixed_VT + self.V * np.einsum( + "kli, klmn, mnj->ij", self._depsdq_fixed_VT, C_T, self._depsdq_fixed_VT + ) + + @material_property + def _d2Fdqdz(self): + """ + Second derivatives of the Helmholtz energy with + respect to composition and z (strain or temperature). + """ + C_T = self.unrelaxed.full_isothermal_stiffness_tensor + d2FdqdT = -self._dSdq + self.V * np.einsum( + "kli, klmn, mn->i", self._depsdq_fixed_VT, C_T, self._depsdT_fixed_volume + ) + + d2Fdqdeps = -self.V * ( + np.einsum("i, mn->imn", self._dPdq, np.eye(3)) + + np.einsum( + "kli, klpq, pqmn->imn", self._depsdq_fixed_VT, C_T, self._dbarepsdeps + ) + ) + d2Fdqdeps = np.array([contract_stresses(arr) for arr in d2Fdqdeps]) + + return np.concatenate((d2Fdqdeps, d2FdqdT[:, np.newaxis]), axis=1) + + @material_property + def _d2Fdqdq_fixed_epsT_pinv(self): + """ + The second derivative of the Helmholtz energy + with respect to the structure parameters at constant + strain and temperature. Often referred to as the + susceptibility matrix. + """ + return np.linalg.pinv(self._d2Fdqdq_fixed_epsT) + + @material_property + def dqdz_relaxed(self): + """ + The change of the structure parameters with respect to + strain and temperature that minimizes the Helmholtz + energy. + """ + return -np.einsum("kl, lj->kj", self._d2Fdqdq_fixed_epsT_pinv, self._d2Fdqdz) + + @material_property + def _d2Fdzdz_Q(self): + """ + Block matrix of V*C_T, V*pi, -c_eps/T + at fixed Q + """ + CT = self.unrelaxed.isothermal_stiffness_tensor + pi = contract_stresses(self.unrelaxed.thermal_stress_tensor) + c_eps = self.unrelaxed.molar_isometric_heat_capacity + V = self.molar_volume + T = self.temperature + return np.block([[V * CT, V * pi[:, np.newaxis]], [V * pi, -c_eps / T]]) + + @material_property + def _d2Fdzdz(self): + """ + Block matrix of V*C_T, V*pi, -c_eps/T + under Helmholtz-minimizing Q + """ + return self._d2Fdzdz_Q + np.einsum( + "ki, kj->ij", self._d2Fdqdz, self.dqdz_relaxed + ) + + # The next three functions provide the three relaxed + # second derivative properties from which other properties + # may be derived. + @material_property + @copy_documentation(AnisotropicMineral.isothermal_stiffness_tensor) + def isothermal_stiffness_tensor(self): + return self._d2Fdzdz[:6, :6] / self.V + + @material_property + @copy_documentation(AnisotropicMineral.thermal_stress_tensor) + def thermal_stress_tensor(self): + return expand_stresses(self._d2Fdzdz[6, :6] / self.V) + + @material_property + @copy_documentation(AnisotropicMineral.molar_isometric_heat_capacity) + def molar_isometric_heat_capacity(self): + return -self._d2Fdzdz[6, 6] * self.temperature + + # The following tensor properties are all + # second derivatives that are calculated in AnisotropicMineral + # (rather than being derived from other properties), + # and must therefore be redefined in relaxed materials + + # The full and Voigt isothermal stiffness tensors, + # full and contracted isentropic stiffness and compliance tensors + # Grueneisen tensor are derived properties + # in AnisotropicMineral + # and so do not need to be redefined. + + @material_property + @copy_documentation(AnisotropicMineral.isothermal_compliance_tensor) + def isothermal_compliance_tensor(self): + return np.linalg.pinv(self.isothermal_stiffness_tensor) + + @material_property + @copy_documentation(AnisotropicMineral.full_isothermal_compliance_tensor) + def full_isothermal_compliance_tensor(self): + S_Voigt = self.isothermal_compliance_tensor + return voigt_notation_to_compliance_tensor(S_Voigt) + + @material_property + @copy_documentation(AnisotropicMineral.thermal_expansivity_tensor) + def thermal_expansivity_tensor(self): + alpha = -np.einsum( + "ijkl, kl", + self.full_isothermal_compliance_tensor, + self.thermal_stress_tensor, + ) + return alpha + + # The following scalar properties are all + # second derivatives that are calculated in AnisotropicMineral + # or Solution (rather than being derived from other properties), + # and must therefore be redefined in relaxed materials + + # The volumetric molar heat capacity and grueneisen parameter are + # derived properties in AnisotropicMineral + # and so do not need to be redefined. + + @material_property + @copy_documentation(AnisotropicMineral.isothermal_compressibility_reuss) + def isothermal_compressibility_reuss(self): + return np.trace(self.isothermal_compressibility_tensor) + + @material_property + @copy_documentation(AnisotropicMineral.isothermal_bulk_modulus_reuss) + def isothermal_bulk_modulus_reuss(self): + return 1.0 / self.isothermal_compressibility_reuss + + @material_property + @copy_documentation(AnisotropicMineral.isentropic_compressibility_reuss) + def isentropic_compressibility_reuss(self): + return np.trace(self.isentropic_compressibility_tensor) + + @material_property + @copy_documentation(AnisotropicMineral.isentropic_bulk_modulus_reuss) + def isentropic_bulk_modulus_reuss(self): + return 1.0 / self.isentropic_compressibility_reuss + + @material_property + @copy_documentation(AnisotropicMineral.thermal_expansivity) + def thermal_expansivity(self): + return np.trace(self.thermal_expansivity_tensor) + + @material_property + @copy_documentation(AnisotropicMineral.molar_heat_capacity_p) + def molar_heat_capacity_p(self): + alpha = self.thermal_expansivity_tensor + pi = self.thermal_stress_tensor + C_p = ( + self.molar_isometric_heat_capacity + - self.V * self.temperature * np.einsum("ij, ij", alpha, pi) + ) + return C_p diff --git a/burnman/classes/solution.py b/burnman/classes/solution.py index 393ef096..4ad6d100 100644 --- a/burnman/classes/solution.py +++ b/burnman/classes/solution.py @@ -18,6 +18,11 @@ from ..utils.reductions import independent_row_indices from ..utils.chemistry import sum_formulae, sort_element_list_to_IUPAC_order +from ..utils.misc import copy_documentation +from ..utils.math import complete_basis + +import copy +from scipy.optimize import minimize class Solution(Mineral): @@ -619,7 +624,7 @@ def reaction_basis(self): to the number of moles of endmember[j] involved in reaction[i]. """ reaction_basis = np.array( - [v[:] for v in self.stoichiometric_matrix.T.nullspace()] + [v[:] for v in self.stoichiometric_matrix.T.nullspace()], dtype=float ) if len(reaction_basis) == 0: @@ -634,6 +639,15 @@ def n_reactions(self): """ return len(self.reaction_basis[:, 0]) + @cached_property + def compositional_basis(self): + """_summary_ + + :return: _description_ + :rtype: _type_ + """ + return complete_basis(self.reaction_basis)[self.n_reactions :] + @cached_property def independent_element_indices(self): """ @@ -705,3 +719,319 @@ def elements(self): SolidSolution = Solution + + +class RelaxedSolution(Solution): + """ + A class implementing a relaxed solution + model. + This class is derived from Solution, + and inherits most of the methods from that class. + + Instantiation of a RelaxedSolution involves passing + an Solution, plus a set of vectors that represent rapid + deformation modes. For example, a solution of MgO, FeHSO and FeLSO + (high and low spin wuestite) can rapidly change proportion of + high spin and low spin iron, and so a single vector should be passed: + np.array([[0., -1., 1.]]) or some multiple thereof. + + States of the mineral can only be queried after setting the + pressure and temperature using set_state() and the composition using + set_composition(). + + This class is available as ``burnman.RelaxedSolution``. + """ + + def __init__( + self, + solution, + relaxation_vectors, + unrelaxed_vectors, + ): + # Make an attribute with the unrelaxed solution + self.unrelaxed = solution + + # The relaxation vectors + self.n_relaxation_vectors = len(relaxation_vectors) + self.n_unrelaxed_vectors = len(unrelaxed_vectors) + + self.dndq = np.array(relaxation_vectors).T + assert len(self.dndq.shape) == 2 + assert len(self.dndq) == self.unrelaxed.n_endmembers + + self.dndx = np.array(unrelaxed_vectors).T + assert len(self.dndx.shape) == 2 + assert len(self.dndx) == self.unrelaxed.n_endmembers + assert ( + len(unrelaxed_vectors) + len(relaxation_vectors) + == self.unrelaxed.n_endmembers + ) + + self.q_initial = np.zeros(len(relaxation_vectors)) + 0.001 + + try: + molar_fractions = solution.molar_fractions + except AttributeError: + molar_fractions = None + + # Give the relaxed solution the same base properties as the + # unrelaxed solution + Solution.__init__( + self, + name=solution.name, + solution_model=solution.solution_model, + molar_fractions=molar_fractions, + ) + + def set_state(self, pressure, temperature, relaxed=True): + """ + Sets the state of the solution. Also relaxes the + structure parameters if set_composition() has already + been used and if the relaxed argument has been set to + True. + + :param pressure: The pressure of the solution [Pa] + :type pressure: float + :param temperature: The temperature of the solution [K] + :type temperature: float + :param relaxed: Whether to minimize the Gibbs energy + of the material by changing the values of the structure + parameters. Defaults to True. + :type relaxed: bool, optional + """ + self.unrelaxed.set_state(pressure, temperature) + + if hasattr(self.unrelaxed, "molar_fractions") and relaxed: + self._relax_at_PTX() + Solution.set_state(self, pressure, temperature) + + def set_composition(self, molar_fractions, q_initial=None, relaxed=True): + """ + Sets the composition of the model. Also relaxes the + structure parameters if set_state() has already + been used and if the relaxed argument has been set to + True. + + :param molar_fractions: Molar fractions of the + independent endmembers corresponding to the + unrelaxed vectors specified during initialisation. + :type molar_fractions: 1D numpy array + :param q_initial: Initial values of the structure parameters. + Defaults to None, in which case the preexisting + initial values are used + (first set to 0.001 during initialisation). + :type q_initial: 1D numpy array, optional + :param relaxed: Whether to minimize the Gibbs energy + of the material by changing the values of the structure + parameters. Defaults to True. + :type relaxed: bool, optional + """ + self.unrelaxed_vectors = molar_fractions + if q_initial is not None: + self.q_initial = np.array(q_initial) + n = np.einsum("ij, j", self.dndq, self.q_initial) + np.einsum( + "ij, j", self.dndx, molar_fractions + ) + + self.unrelaxed.set_composition(n) + + if self.unrelaxed.pressure is not None and relaxed: + self._relax_at_PTX() + Solution.set_composition(self, n) + self.unrelaxed.set_composition(self.molar_fractions) + + def _relax_at_PTX(self): + """ + Minimizes the Gibbs energy at constant pressure and + temperature by changing the structural parameters. + + Run during set_state() and set_composition(), as long as both + state and composition have already been set. This function + should not generally be needed by the user. + """ + n0 = copy.copy(self.unrelaxed.molar_fractions) + + def G_func(dq): + n = n0 + np.einsum("ij, j", self.dndq, dq) + self.unrelaxed.set_composition(n) + return self.unrelaxed.molar_gibbs + + sol = minimize(G_func, np.zeros(len(self.dndq[0])), method="Nelder-Mead") + assert sol.success + n = n0 + np.einsum("ij, j", self.dndq, sol.x) + self.unrelaxed.set_composition(n) + Solution.set_composition(self, n) + + def set_state_with_volume( + self, volume, temperature, pressure_guesses=[0.0e9, 10.0e9] + ): + """ + This function acts similarly to set_state, but takes volume and + temperature as input to find the pressure. In order to ensure + self-consistency, this function does not use any pressure functions + from the material classes, but instead finds the pressure using the + brentq root-finding and Nelder-Mead minimization methods. + + Composition should have been set before this function is used. + + :param volume: The desired molar volume of the mineral [m^3]. + :type volume: float + + :param temperature: The desired temperature of the mineral [K]. + :type temperature: float + + :param pressure_guesses: A list of floats denoting the initial + low and high guesses for bracketing of the pressure [Pa]. + 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 [0.e9, 10.e9]. + :type pressure_guesses: list + """ + ss = self.unrelaxed + n0 = copy.copy(ss.molar_fractions) + + # Initial set_state without changing structural parameters + # to estimate pressure + ss.set_state_with_volume(volume, temperature, pressure_guesses) + # Store in a mutable so that it can be updated in the loop + pressure = [ss.pressure] + + def F_func(dq): + n = n0 + np.einsum("ij, j", self.dndq, dq) + ss.set_composition(n) + ss.set_state_with_volume( + volume, temperature, [pressure[0] - 1.0e9, pressure[0] + 1.0e9] + ) + pressure[0] = ss.pressure + return ss.molar_helmholtz + + sol = minimize(F_func, np.zeros(len(self.dndq[0])), method="Nelder-Mead") + assert sol.success + + # Run one more time with root + F_func(sol.x) + self.unrelaxed.set_composition(ss.molar_fractions) + Solution.set_composition(self, ss.molar_fractions) + self.set_state(pressure[0], temperature) + + @material_property + def _d2Gdqdq_fixed_PT(self): + """ + Gibbs structural hessian + calculated at constant pressure and temperature. + """ + return np.einsum( + "ij, ik, jl->kl", self.unrelaxed.gibbs_hessian, self.dndq, self.dndq + ) + + @material_property + def _d2Gdqdz(self): + """ + Second derivatives of the Helmholtz energy with + respect to composition and z (pressure or temperature). + """ + dVdq = np.einsum("i, ij->j", self.unrelaxed.partial_volumes, self.dndq) + dSdq = np.einsum("i, ij->j", self.unrelaxed.partial_entropies, self.dndq) + + return np.array([dVdq, -dSdq]).T + + @material_property + def _d2Gdqdq_fixed_PT_pinv(self): + """ + The second derivative of the Helmholtz energy + with respect to the structure parameters at constant + strain and temperature. Often referred to as the + susceptibility matrix. + """ + return np.linalg.pinv(self._d2Gdqdq_fixed_PT) + + @material_property + def dqdz_relaxed(self): + """ + The change of the structure parameters with respect to + strain and temperature that minimizes the Helmholtz + energy. + """ + return -np.einsum("kl, lj->kj", self._d2Gdqdq_fixed_PT_pinv, self._d2Gdqdz) + + @material_property + def _d2Gdzdz_Q(self): + """ + Block matrix of -V*beta_TR, V*alpha, -c_p/T + at fixed Q + """ + beta_TR = self.unrelaxed.isothermal_compressibility_reuss + alpha = self.unrelaxed.thermal_expansivity + c_p = self.unrelaxed.molar_heat_capacity_p + V = self.molar_volume + T = self.temperature + return np.array([[-V * beta_TR, V * alpha], [V * alpha, -c_p / T]]) + + @material_property + def _d2Gdzdz(self): + """ + Block matrix of -V*beta_TR, V*alpha, -c_p/T + under Gibbs-minimizing Q + """ + return self._d2Gdzdz_Q + np.einsum( + "ki, kj->ij", self._d2Gdqdz, self.dqdz_relaxed + ) + + # The following scalar properties are all + # second derivatives that are calculated in AnisotropicMineral + # or Solution (rather than being derived from other properties), + # and must therefore be redefined in relaxed materials + + # The volumetric molar heat capacity and grueneisen parameter are + # derived properties in AnisotropicMineral + # and so do not need to be redefined. + + @material_property + @copy_documentation(Solution.isothermal_compressibility_reuss) + def isothermal_compressibility_reuss(self): + return -self._d2Gdzdz[0, 0] / self.V + + @material_property + @copy_documentation(Solution.thermal_expansivity) + def thermal_expansivity(self): + return self._d2Gdzdz[0, 1] / self.V + + @material_property + @copy_documentation(Solution.molar_heat_capacity_p) + def molar_heat_capacity_p(self): + return -self._d2Gdzdz[1, 1] * self.T + + @material_property + @copy_documentation(Solution.isothermal_bulk_modulus_reuss) + def isothermal_bulk_modulus_reuss(self): + return 1.0 / self.isothermal_compressibility_reuss + + @material_property + @copy_documentation(Solution.molar_heat_capacity_v) + def molar_heat_capacity_v(self): + return ( + self.molar_heat_capacity_p + - self.molar_volume + * self.temperature + * self.thermal_expansivity + * self.thermal_expansivity + * self.isothermal_bulk_modulus_reuss + ) + + @material_property + @copy_documentation(Solution.isentropic_bulk_modulus_reuss) + def isentropic_bulk_modulus_reuss(self): + if self.temperature < 1.0e-10: + return self.isothermal_bulk_modulus_reuss + else: + return ( + self.isothermal_bulk_modulus_reuss + * self.molar_heat_capacity_p + / self.molar_heat_capacity_v + ) + + @material_property + @copy_documentation(Solution.isentropic_compressibility_reuss) + def isentropic_compressibility_reuss(self): + return 1.0 / self.isentropic_bulk_modulus_reuss diff --git a/burnman/tools/solution.py b/burnman/tools/solution.py index bf58c4a7..269e7db1 100644 --- a/burnman/tools/solution.py +++ b/burnman/tools/solution.py @@ -17,6 +17,7 @@ SubregularSolution, ) from ..utils.chemistry import site_occupancies_to_strings +from ..utils.math import complete_basis def _decompose_3D_matrix(W): @@ -170,23 +171,6 @@ def _subregular_matrix_conversion( return (new_endmember_excesses, new_binary_terms, new_ternary_terms) -def complete_basis(basis): - """ - Creates a full basis by filling remaining rows with - rows of the identity matrix with row indices not - in the column pivot list of the basis RREF - """ - - n, m = basis.shape - if n < m: - pivots = list(Matrix(basis).rref()[1]) - return np.concatenate( - (basis, np.identity(m)[[i for i in range(m) if i not in pivots], :]), axis=0 - ) - else: - return basis - - def transform_solution_to_new_basis( solution, new_basis, diff --git a/burnman/utils/math.py b/burnman/utils/math.py index 3a9a258f..d142646d 100644 --- a/burnman/utils/math.py +++ b/burnman/utils/math.py @@ -495,6 +495,37 @@ def array_to_rational_matrix(array): return Matrix([[Rational(v).limit_denominator(1000) for v in row] for row in array]) +def complete_basis(basis, flip_columns=True): + """ + Creates a full basis by filling remaining rows with + selected rows of the identity matrix that have indices + not in the column pivot list of the basis RREF + + :param basis: A partial basis + :type basis: numpy array + :param flip_columns: whether to calculate the RREF + from a flipped basis. This can be useful, for example, + when dependent columns are known to be further to the + right in the basis. defaults to True + :type flip_columns: bool, optional + :return: basis + :rtype: numpy array + """ + + n, m = basis.shape + if n < m: + if flip_columns: + pivots = list(m - 1 - np.array(Matrix(np.flip(basis, axis=1)).rref()[1])) + else: + pivots = list(Matrix(basis)).rref()[1] + + return np.concatenate( + (basis, np.identity(m)[[i for i in range(m) if i not in pivots], :]), axis=0 + ) + else: + return basis + + def generate_complete_basis(incomplete_basis, array): """ Given a 2D array with independent rows and a second 2D array that spans a diff --git a/tests/test_solidsolution.py b/tests/test_solidsolution.py index 73683975..223214f9 100644 --- a/tests/test_solidsolution.py +++ b/tests/test_solidsolution.py @@ -1218,6 +1218,20 @@ def test_evaluate(self): self.assertArraysAlmostEqual(G1, G2) + def test_relaxed_solution(self): + opx = burnman.minerals.JH_2015.orthopyroxene() + ropx = burnman.RelaxedSolution(opx, opx.reaction_basis, opx.compositional_basis) + + ropx.set_composition([0.2, 0.2, 0.1, 0.1, 0.2, 0.2]) + + for T in [300.0, 1000.0, 2000.0]: + ropx.set_state(1.0e9, T) + mu = np.einsum("ij, j", opx.reaction_basis, ropx.partial_gibbs)[0] + + # Check that the reaction affinity is essentially zero + # Give 5 J tolerance. + self.assertTrue(np.abs(mu) < 5.0) + if __name__ == "__main__": unittest.main()