# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for
# the Earth and Planetary Sciences
# Copyright (C) 2012 - 2025 by the BurnMan team, released under the GNU
# GPL v2 or later.
import numpy as np
from itertools import product
from scipy.linalg import lu_factor, lu_solve
from collections import namedtuple
import logging
from ..constants import logish_eps
from ..optimize.nonlinear_solvers import damped_newton_solve, TerminationCode
from ..classes.solution import Solution
P_scaling = 1.0e9 # Pa
T_scaling = 1.0e3 # K
S_scaling = 1.0e3 # J/K
V_scaling = 1.0e-6 # m^3
G_scaling = 1.0e3 # J/mol
X_scaling = 1.0 # mole amount
[docs]
def calculate_constraints(assemblage, n_free_compositional_vectors):
"""
This function calculates the linear inequality constraints bounding
the valid parameter space for a given assemblage.
The constraints are as follows:
- Pressure and temperature must be positive
- All phase fractions must be positive
- All site-species occupancies must be positive
The constraints are stored in a vector (b) and matrix (A).
The sign convention is chosen such that the constraint is satisfied
if A.x + b < eps.
:param assemblage: The assemblage for which the constraints are calculated.
:type assemblage: :class:`burnman.Composite`
:returns: The constraints vector and matrix.
:rtype: tuple
"""
bounds = []
n_constraints = 0
for i, n in enumerate(assemblage.endmembers_per_phase):
n_constraints += 1
if n == 1:
bounds.append(np.array([[]]))
else:
bounds.append(assemblage.phases[i].solution_model.endmember_occupancies)
n_constraints += len(bounds[-1][0])
c_vector = np.zeros((n_constraints + 2))
c_matrix = np.zeros(
(n_constraints + 2, assemblage.n_endmembers + 2 + n_free_compositional_vectors)
) # includes P, T
c_matrix[0, 0] = -1 # P>0
c_matrix[1, 1] = -1 # T>0
cidx = 2 # index of current compositional constraint
pidx = 0 # starting index of current phase
for i, n in enumerate(assemblage.endmembers_per_phase):
m = len(bounds[i][0])
# The first endmember proportion is not a free variable
# (all endmembers in any solution must sum to one)
# Re-express the constraints without the first endmember
c_matrix[cidx, pidx + 2] = -1.0 # need phase proportions > 0
cidx += 1
if m != 0:
c_vector[cidx : cidx + m] = -bounds[i][0]
c_matrix[cidx : cidx + m, pidx + 1 + 2 : pidx + n + 2] = (
np.einsum("i, j", bounds[i][0], np.ones_like(bounds[i][1:, 0]))
- bounds[i].T[:, 1:]
)
cidx += m
pidx += n
return c_vector, c_matrix
[docs]
def get_parameters(assemblage, n_free_compositional_vectors=0):
"""
Gets the starting parameters vector (x) for the current equilibrium problem.
These are:
- pressure
- temperature
- absolute amount of each phase. if a phase is a solution
with >1 endmember, the following parameters are the mole fractions
of the independent endmembers in the solution, except for the first
endmember (as the mole fractions must sum to one).
:param assemblage: The assemblage to be equilibrated.
:type assemblage: :class:`burnman.Composite`
:returns: The current values of all the parameters.
:rtype: numpy.array
"""
params = np.zeros(assemblage.n_endmembers + 2 + n_free_compositional_vectors)
n_moles_phase = assemblage.number_of_moles * np.array(assemblage.molar_fractions)
try:
params[:2] = [
assemblage.pressure / P_scaling,
assemblage.temperature / T_scaling,
]
except AttributeError:
raise Exception("You need to set_state before getting parameters")
j = 2
for i, ph in enumerate(assemblage.phases):
params[j] = n_moles_phase[i]
if isinstance(ph, Solution):
params[j + 1 : j + assemblage.endmembers_per_phase[i]] = assemblage.phases[
i
].molar_fractions[1:]
j += assemblage.endmembers_per_phase[i]
return params
[docs]
def get_endmember_amounts(assemblage):
"""
Gets the absolute amounts of all the endmembers in the solution.
:param assemblage: The assemblage to be equilibrated.
:type assemblage: :class:`burnman.Composite`
:returns: The current amounts of all the endmembers.
:rtype: numpy.array
"""
phase_amounts = assemblage.number_of_moles * assemblage.molar_fractions
amounts = np.empty(assemblage.n_endmembers)
j = 0
for i, ph in enumerate(assemblage.phases):
if isinstance(ph, Solution):
amounts[j : j + assemblage.endmembers_per_phase[i]] = (
phase_amounts[i] * assemblage.phases[i].molar_fractions
)
else:
amounts[j] = phase_amounts[i]
j += assemblage.endmembers_per_phase[i]
return amounts
[docs]
def set_compositions_and_state_from_parameters(assemblage, parameters):
"""
Sets the phase compositions, amounts and state of the assemblage
from a list of parameter values.
:param assemblage: The assemblage to be equilibrated.
:type assemblage: :class:`burnman.Composite`
:param parameters: The current parameter values.
:type parameters: numpy.array
"""
i = 2
phase_amounts = np.zeros(len(assemblage.phases))
for phase_idx, ph in enumerate(assemblage.phases):
phase_amounts[phase_idx] = parameters[i]
if isinstance(ph, Solution):
n_mbrs = len(ph.endmembers)
f = [0.0] * n_mbrs
f[1:] = parameters[i + 1 : i + n_mbrs]
f[0] = 1.0 - sum(f)
ph.set_composition(f)
i += n_mbrs
else:
i += 1
assert np.all(phase_amounts > -1.0e-8)
phase_amounts = np.abs(phase_amounts)
assemblage.number_of_moles = sum(phase_amounts)
assemblage.set_fractions(phase_amounts / assemblage.number_of_moles)
# We set_state last because setting the state may depend on the
# current composition of the assemblage. This would occur, for example,
# in the case of relaxed composites, where the relaxed state depends
# on the current bulk composition.
# Note that in the case of relaxed composites, the equilibrate function
# actually equilibrates the unrelaxed assemblage, and then copies the
# state to the relaxed assemblage afterwards.
assemblage.set_state(parameters[0] * P_scaling, parameters[1] * T_scaling)
return None
[docs]
def default_F_tolerances(assemblage, equality_constraints, n_atoms):
"""
Returns the default tolerances for each component of F(x).
"""
F_tolerances = np.zeros(assemblage.n_endmembers + len(equality_constraints))
i = 0
for i, (type_c, _) in enumerate(equality_constraints):
if type_c == "P":
F_tolerances[i] = 1.0 / P_scaling # Pa
elif type_c == "T":
F_tolerances[i] = 1.0e-6 / T_scaling # K
elif type_c == "S":
F_tolerances[i] = 1.0e-8 / S_scaling # J/K
elif type_c == "V":
# Typical volume of 1 mole of atoms is ~1e-5 m^3
# We want a much smaller tolerance than this
F_tolerances[i] = 1.0e-11 * n_atoms / V_scaling # m^3
elif type_c == "X":
F_tolerances[i] = 1.0e-8 / X_scaling # [no units]
else:
raise Exception("constraint type not recognised")
i += 1
# Next set of tolerances are for reaction affinities (1 J/mol)
# and bulk composition (1e-8)
F_tolerances[i : i + assemblage.n_reactions] = 1.0 / G_scaling
F_tolerances[i + assemblage.n_reactions :] = 1.0e-8 / X_scaling
return F_tolerances
[docs]
def F(
x,
assemblage,
equality_constraints,
reduced_composition_vector,
reduced_free_composition_vectors,
):
"""
The vector-valued function for which the root is sought.
The first two vector values depend on the
equality_constraints chosen. For example, if
- eq[i][0] = 'P', F[i] = P - eq[i][1]
- eq[i][0] = 'T', F[i] = T - eq[i][1]
- eq[i][0] = 'S', F[i] = entropy - eq[i][1]
- eq[i][0] = 'V', F[i] = volume - eq[i][1]
- eq[i][0] = 'X', np.dot(eq[i][1][0], x) - eq[i][1][1]
The next set of vector values correspond to the reaction affinities.
The final set of vector values correspond to the bulk
composition constraints.
:param x: Parameter values for the equilibrium problem to be solved.
:type x: numpy array
:param assemblage: The assemblage to be equilibrated.
:type assemblage: :class:`burnman.Composite`
:param equality_constraints: A list of the equality constraints
(see above for valid formats).
:type equality_constraints: list of lists
:param reduced_composition_vector: The amounts of the independent
elements.
:type reduced_composition_vector: numpy.array
:param reduced_free_composition_vectors: The amounts of the
independent elements in each of the free_compositional_vectors.
:type reduced_free_composition_vectors: 2D numpy.array
:returns: The vector corresponding to F(x).
:rtype: numpy.array
"""
set_compositions_and_state_from_parameters(assemblage, x)
new_endmember_amounts = get_endmember_amounts(assemblage)
# We want to find the root of the following equations
n_equality_constraints = len(equality_constraints)
eqns = np.zeros((assemblage.n_endmembers + n_equality_constraints))
i = 0
for i, (type_c, eq_c) in enumerate(equality_constraints):
if type_c == "P":
eqns[i] = (assemblage.pressure - eq_c) / P_scaling
elif type_c == "T":
eqns[i] = (assemblage.temperature - eq_c) / T_scaling
elif type_c == "S":
eqns[i] = (assemblage.entropy - eq_c) / S_scaling
elif type_c == "V":
eqns[i] = (assemblage.volume - eq_c) / V_scaling
elif type_c == "X":
eqns[i] = (np.dot(eq_c[0], x) - eq_c[1]) / X_scaling # i.e. Ax = b
else:
raise Exception("constraint type not recognised")
i += 1
if n_equality_constraints > 2:
new_reduced_composition_vector = (
reduced_composition_vector
+ x[2 - n_equality_constraints :].dot(reduced_free_composition_vectors)
/ X_scaling
)
else:
new_reduced_composition_vector = reduced_composition_vector / X_scaling
eqns[i : i + assemblage.n_reactions] = assemblage.reaction_affinities / G_scaling
eqns[i + assemblage.n_reactions :] = (
np.dot(assemblage.reduced_stoichiometric_array.T, new_endmember_amounts)
- new_reduced_composition_vector
) / X_scaling
return eqns
[docs]
def jacobian(x, assemblage, equality_constraints, reduced_free_composition_vectors):
"""
The Jacobian of the vector-valued function :math:`F` for which the
root is sought (:math:`\\partial F / \\partial x`).
See documentation for :func:`F` and :func:`get_parameters`
(which return :math:`F` and :math:`x` respectively) for more details.
:param x: Parameter values for the equilibrium problem to be solved.
:type x: numpy.array
:param assemblage: The assemblage to be equilibrated.
:type assemblage: :class:`burnman.Composite`
:param equality_constraints: A list of the equality constraints
(see documentation for :func:`burnman.tools.equilbration.F`).
:type equality_constraints: list of lists
:param reduced_free_composition_vectors: The amounts of the
independent elements in each of the free_compositional_vectors.
:type reduced_free_composition_vectors: 2D numpy array
:returns: The Jacobian for the equilibrium problem.
:rtype: 2D numpy.array
"""
# The solver always calls the Jacobian with the same
# x parameters as used previously for the root functions
# Therefore we don't need to set compositions or state again here.
# First, we find out the effect of the two constraint parameters F[:2]
# on the pressure (x[0]) and temperature (x[1]):
n_equality_constraints = len(equality_constraints)
jacobian = np.zeros(
(
assemblage.n_endmembers + n_equality_constraints,
assemblage.n_endmembers + n_equality_constraints,
)
)
ic = 0
for ic, (type_c, eq_c) in enumerate(equality_constraints):
if type_c == "P": # dP/dx
jacobian[ic, 0] = 1.0 / P_scaling # jacobian[i, j!=0] = 0
elif type_c == "T": # dT/dx
jacobian[ic, 1] = 1.0 / T_scaling # jacobian[i, j!=1] = 0
elif type_c == "S": # dS/dx
# dS/dP = -aV, dS/dT = Cp/T
jacobian[ic, 0:2] = [
-assemblage.alpha * assemblage.volume / S_scaling,
assemblage.heat_capacity_p / assemblage.temperature / S_scaling,
]
j = 2
for k, n in enumerate(assemblage.endmembers_per_phase):
jacobian[ic, j] = assemblage.phases[k].molar_entropy / S_scaling
if n > 1: # for solutions with >1 endmember
jacobian[ic, j + 1 : j + n] = (
assemblage.number_of_moles
* assemblage.molar_fractions[k]
* (
assemblage.phases[k].partial_entropies[1:]
- assemblage.phases[k].partial_entropies[0]
)
/ S_scaling
)
j += n
elif type_c == "V": # dV/dx
# dV/dP = -V/K_T, dV/dT = aV
jacobian[ic, 0:2] = [
-assemblage.volume
/ assemblage.isothermal_bulk_modulus_reuss
/ V_scaling,
assemblage.volume * assemblage.alpha / V_scaling,
]
j = 2
for k, n in enumerate(assemblage.endmembers_per_phase):
jacobian[ic, j] = assemblage.phases[k].molar_volume / V_scaling
if n > 1: # for solutions with >1 stable endmember
jacobian[ic, j + 1 : j + n] = (
assemblage.number_of_moles
* assemblage.molar_fractions[k]
* (
assemblage.phases[k].partial_volumes[1:]
- assemblage.phases[k].partial_volumes[0]
)
/ V_scaling
)
j += n
elif type_c == "X":
jacobian[ic, :] = eq_c[0] / X_scaling
else:
raise Exception("constraint type not recognised")
ic += 1
# Next, let's get the effect of pressure and temperature
# on each of the independent reactions
# i.e. dF(i, reactions)/dx[0] and dF(i, reactions)/dx[1]
partial_volumes_vector = np.zeros((assemblage.n_endmembers))
partial_entropies_vector = np.zeros((assemblage.n_endmembers))
j = 0
for i, n in enumerate(assemblage.endmembers_per_phase):
if n == 1: # for endmembers
partial_volumes_vector[j] = assemblage.phases[i].molar_volume
partial_entropies_vector[j] = assemblage.phases[i].molar_entropy
else: # for solutions
partial_volumes_vector[j : j + n] = assemblage.phases[i].partial_volumes
partial_entropies_vector[j : j + n] = assemblage.phases[i].partial_entropies
j += n
reaction_volumes = np.dot(assemblage.reaction_basis, partial_volumes_vector)
reaction_entropies = np.dot(assemblage.reaction_basis, partial_entropies_vector)
# dGi/dP = deltaVi; dGi/dT = -deltaSi
jacobian[ic : ic + len(reaction_volumes), 0] = reaction_volumes / G_scaling
jacobian[ic : ic + len(reaction_volumes), 1] = -reaction_entropies / G_scaling
# Pressure and temperature have no effect on the bulk
# compositional constraints
# i.e. dF(i, bulk)/dx[0] and dF(i, bulk)/dx[1] = 0
# Finally, let's build the compositional Hessian d2G/dfidfj = dmui/dfj
# where fj is the fraction of endmember j in a phase
phase_amounts = np.array(assemblage.molar_fractions) * assemblage.number_of_moles
comp_hessian = np.zeros((assemblage.n_endmembers, assemblage.n_endmembers))
dfi_dxj = np.zeros((assemblage.n_endmembers, assemblage.n_endmembers))
dpi_dxj = np.zeros((assemblage.n_endmembers, assemblage.n_endmembers))
j = 0
for i, n in enumerate(assemblage.endmembers_per_phase):
if n == 1:
# changing the amount (p) of a pure phase
# does not change its fraction in that phase,
# so dfi_dxj remains unchanged
dpi_dxj[j, j] = 1.0
else:
comp_hessian[j : j + n, j : j + n] = assemblage.phases[i].gibbs_hessian
# x[0] = p(phase) and x[1:] = f[1:] - f[0]
# Therefore
# df[0]/dx[0] = 0
# df[0]/dx[1:] = -1
# (because changing the fraction of any endmember
# depletes the fraction of the first endmember)
# df[1:]/dx[1:] = 1 on diagonal, 0 otherwise
# (because all other fractions are independent of each other)
dfi_dxj[j : j + n, j : j + n] = np.eye(n)
dfi_dxj[j, j : j + n] -= 1.0
# Total amounts of endmembers (p) are the fractions
# multiplied by the amounts of their representative phases
dpi_dxj[j : j + n, j : j + n] = (
dfi_dxj[j : j + n, j : j + n] * phase_amounts[i]
)
# The derivative of the amount of each endmember with respect
# to the amount of each phase is equal to the molar fractions
# of the endmembers.
dpi_dxj[j : j + n, j] = assemblage.phases[i].molar_fractions
j += n
# dfi_dxj converts the endmember hessian to the parameter hessian.
reaction_hessian = assemblage.reaction_basis.dot(comp_hessian).dot(dfi_dxj)
bulk_hessian = assemblage.reduced_stoichiometric_array.T.dot(dpi_dxj)
if reaction_hessian.shape[0] > 0:
jacobian[ic:, 2 : 2 + len(reaction_hessian[0])] = np.concatenate(
(reaction_hessian / G_scaling, bulk_hessian / X_scaling)
)
else:
jacobian[ic:, 2 : 2 + len(bulk_hessian[0])] = bulk_hessian / X_scaling
if len(reduced_free_composition_vectors) > 0:
jacobian[
-reduced_free_composition_vectors.shape[1] :, 2 + len(reaction_hessian[0]) :
] = (-reduced_free_composition_vectors.T / X_scaling)
x_scaling_array = np.ones_like(x)
x_scaling_array[0] = P_scaling # Pa
x_scaling_array[1] = T_scaling # K
return jacobian * x_scaling_array[np.newaxis, :]
[docs]
def lambda_bounds(dx, x, endmembers_per_phase):
"""
Returns the lambda bounds for the damped affine invariant modification
to Newton's method for nonlinear problems (Deuflhard, 1974;1975;2004).
:param dx: The proposed newton step.
:type dx: numpy.array
:param x: Parameter values for the equilibrium problem to be solved.
:type x: numpy.array
:param endmembers_per_phase: A list of the number of endmembers in each phase.
:type endmembers_per_phase: list of int
:returns: Minimum and maximum allowed fractions of the full newton step (dx).
:rtype: tuple of floats
"""
max_steps = np.ones((len(x))) * 100000.0
# first two constraints are P and T, scaled appropriately
max_steps[0:2] = np.array(
[20.0e9 / P_scaling, 500.0 / T_scaling]
) # biggest reasonable P and T steps
j = 2
for i, n in enumerate(endmembers_per_phase):
# if the phase fraction constraint would be broken,
# set a step that is marginally smaller
if x[j] + dx[j] < 0.0:
max_steps[j] = max(x[j] * 0.999, 0.001)
max_steps[j + 1 : j + n] = [
max(xi * 0.99, 0.01) for xi in x[j + 1 : j + n]
] # maximum compositional step
j += n
max_lmda = min(
[
1.0 if step <= max_steps[i] else max_steps[i] / step
for i, step in enumerate(np.abs(dx))
]
)
if np.isnan(max_lmda):
raise ValueError(
"NaN in lambda bounds calculation.\n"
f"dx: {dx}\n"
f"x: {x}\n"
f"endmembers_per_phase: {endmembers_per_phase}."
)
return (1.0e-8, max_lmda)
[docs]
def phase_fraction_constraints(phase, assemblage, fractions, prm):
"""
Converts a phase fraction constraint into standard linear form
that can be processed by the root finding problem.
We start with a single fraction or an array of fractions
for a particular phase (:math:`n_p / \\sum n = f`).
These are then converted into the "X" form of constraint
by multiplying by :math:`\\sum n` and moving all terms to
the LHS of the equation:
:math:`-f n_0 - f n_1 - \\ldots + (1-f) n_p - \\ldots = 0`
This form is less readable, but easier to use as a constraint
in a nonlinear solve.
:param phase: The phase for which the fraction is to be constrained
:type phase: :class:`burnman.Solution` or :class:`burnman.Mineral`
:param assemblage: The assemblage to be equilibrated.
:type assemblage: :class:`burnman.Composite`
:param fractions: The phase fractions to be satified at equilibrium.
:type fractions: numpy.array
:param prm: A tuple with attributes n_parameters
(the number of parameters for the current equilibrium problem)
and phase_amount_indices (the indices of the parameters that
correspond to phase amounts).
:type prm: namedtuple
:returns: The phase fraction constraints.
:rtype: list
"""
phase_idx = assemblage.phases.index(phase)
constraints = []
for fraction in fractions:
constraints.append(["X", [np.zeros((prm.n_parameters)), 0.0]])
constraints[-1][-1][0][prm.phase_amount_indices] = -fraction
constraints[-1][-1][0][prm.phase_amount_indices[phase_idx]] += 1.0
return constraints
[docs]
def phase_composition_constraints(phase, assemblage, constraints, prm):
"""
Converts a phase composition constraint into standard linear form
that can be processed by the root finding problem.
We start with constraints in the form (site_names, n, d, v), where
:math:`(n x)/ (d x) = v` and :math:`n` and :math:`d` are fixed vectors of
site coefficients. So, one could for example choose a constraint
([Mg_A, Fe_A], [1., 0.], [1., 1.], [0.5]) which would
correspond to equal amounts Mg and Fe on the A site.
This function converts the user-defined vectors of site constraints
:math:`n` and :math:`d` into vectors of endmember proportion
constraints :math:`n'` and :math:`d'`, such that
:math:`(n' x)/ (d' x) = v`. This is done via linear transformation
using the site occupancy matrix provided by :class:`burnman.Solution`.
By multiplying by the denominator, we have the following scalar
comparison: :math:`(n' x) = v (d' x)`
The equilibration function does not use the proportion of
the first endmember (as the endmember proportions must sum to one),
and so we split :math:`x`, :math:`n'` and :math:`d'` into the first
element and following elements:
:math:`(n'_0 x_0 + n'_i x_i) = v (d'_0 x_0 + d'_i x_i)`
where :math:`i` is taken over all elements apart from the first.
With some more rearranging we can express the constraint in standard
linear form:
:math:`(n'_0 (1 - \\sum_j x_j) + n'_i x_i) = v (d'_0 (1 - \\sum_j x_j) + d'_i x_i)`
:math:`(n'_0 + (n'_i - 1_i n'_0) x_i) = v (d'_0 + (d'_i - 1_i d'_0) x_i)`
:math:`(((n'_i - 1_i n'_0) - v(d'_i - 1_i d'_0)) x_i) = (v d'_0 - n'_0)`
This form is less readable, but easier to use as a constraint
in a nonlinear solve.
:param phase: The phase for which the composition is to be constrained.
:type phase: :class:`burnman.Solution`
:param assemblage: The assemblage to be equilibrated.
:type assemblage: :class:`burnman.Composite`
:param constraints: The desired constraints in the form:
site_names (list of strings), numerator (numpy.array),
denominator (numpy.array), values (numpy.array).
:type constraints: tuple
:returns: The phase composition constraints in standard form.
:rtype: list
"""
phase_idx = assemblage.phases.index(phase)
site_names, numerator, denominator, values = constraints
site_indices = [phase.solution_model.site_names.index(name) for name in site_names]
noccs = phase.solution_model.endmember_noccupancies
# Converts site constraints into endmember constraints
# Ends up with shape (n_endmembers, 2)
endmembers = np.dot(noccs[:, site_indices], np.array([numerator, denominator]).T)
numer0, denom0 = endmembers[0]
endmembers -= endmembers[0]
numer, denom = endmembers.T[:, 1:]
# We start from the correct index
start_idx = sum(assemblage.endmembers_per_phase[:phase_idx]) + 3
n_indices = assemblage.endmembers_per_phase[phase_idx] - 1
x_constraints = []
for v in values:
f = v * denom0 - numer0
x_constraints.append(["X", [np.zeros((prm.n_parameters)), f]])
x_constraints[-1][1][0][start_idx : start_idx + n_indices] = numer - v * denom
return x_constraints
[docs]
def get_equilibration_parameters(assemblage, composition, free_compositional_vectors):
"""
Builds a named tuple containing the parameter names and
various other parameters needed by the equilibrium solve.
:param assemblage: The assemblage to be equilibrated.
:type assemblage: :class:`burnman.Composite`
:param composition: The bulk composition for the equilibrium problem.
:type composition: dict
:param free_compositional_vectors: The bulk compositional
degrees of freedom for the equilibrium problem.
:type free_compositional_vectors: list of dictionaries
:returns: A tuple with attributes:
- parameter_names: list of strings, the names of the parameters
for the current equilibrium problem
- default_tolerances: numpy.array, the default tolerances for each parameter
- n_parameters: int, the number of parameters for the current equilibrium problem
- phase_amount_indices: list of int, the indices of the parameters that
correspond to phase amounts
- bulk_composition_vector: numpy.array, the bulk composition vector
- free_compositional_vectors: 2D numpy.array, the free compositional vectors
- reduced_composition_vector: numpy.array, the bulk composition vector
reduced to the independent elements
- reduced_free_composition_vectors: 2D numpy.array, the free compositional vectors
reduced to the independent elements
- constraint_vector: numpy.array, vector for the linear constraints bounding the
valid parameter space (b in Ax + b < eps)
- constraint_matrix: 2D numpy.array, matrix for the linear constraints bounding the
valid parameter space (A in Ax + b < eps)
:rtype: namedtuple
"""
# Initialize a named tuple for the equilibration parameters
prm = namedtuple("assemblage_parameters", [])
# Process parameter names
prm.parameter_names = ["Pressure (Pa)", "Temperature (K)"]
prm.default_tolerances = [
1.0 / P_scaling,
1.0e-6 / T_scaling,
] # Tolerances scaled to match F scaling
for i, n in enumerate(assemblage.endmembers_per_phase):
prm.parameter_names.append("x({0})".format(assemblage.phases[i].name))
prm.default_tolerances.append(1.0e-9)
if n > 1:
p_names = [
"p({0} in {1})".format(n, assemblage.phases[i].name)
for n in assemblage.phases[i].endmember_names[1:]
]
prm.parameter_names.extend(p_names)
prm.default_tolerances.extend([1.0e-9] * len(p_names))
n_free_compositional_vectors = len(free_compositional_vectors)
for i in range(n_free_compositional_vectors):
prm.parameter_names.append(f"v_{i}")
prm.default_tolerances.append(1.0e-9)
prm.default_tolerances = np.array(prm.default_tolerances)
prm.n_parameters = len(prm.parameter_names)
prm.phase_amount_indices = [
i for i in range(len(prm.parameter_names)) if "x(" in prm.parameter_names[i]
]
# Find the bulk composition vector
prm.bulk_composition_vector = np.array(
[composition[e] for e in assemblage.elements]
)
if n_free_compositional_vectors > 0:
prm.free_compositional_vectors = np.array(
[
[
(
free_compositional_vectors[i][e]
if e in free_compositional_vectors[i]
else 0.0
)
for e in assemblage.elements
]
for i in range(n_free_compositional_vectors)
]
)
else:
prm.free_compositional_vectors = np.empty((0, len(assemblage.elements)))
if assemblage.compositional_null_basis.shape[0] != 0:
if (
np.abs(
assemblage.compositional_null_basis.dot(prm.bulk_composition_vector)[0]
)
> 1.0e-12
):
raise Exception(
"The bulk composition is not within the "
"compositional space of the assemblage"
)
prm.reduced_composition_vector = prm.bulk_composition_vector[
assemblage.independent_element_indices
]
prm.reduced_free_composition_vectors = prm.free_compositional_vectors[
:, assemblage.independent_element_indices
]
prm.constraint_vector, prm.constraint_matrix = calculate_constraints(
assemblage, n_free_compositional_vectors
)
return prm
[docs]
def process_eq_constraints(equality_constraints, assemblage, prm):
"""
A function that processes the equality constraints
into a form that can be processed by the F and jacobian functions.
This function has two main tasks: it turns phase_fraction and
phase_composition constraints into standard linear constraints in the
solution parameters. It also turns vector-valued constraints into a
list of scalar-valued constraints.
:param equality_constraints: A list of equality constraints.
For valid types of constraints, please see the documentation for
:func:`burnman.equilibrate`.
:type equality_constraints: list
:param assemblage: The assemblage to be equilibrated.
:type assemblage: :class:`burnman.Composite`
:param prm: A tuple with attributes n_parameters
(the number of parameters for the current equilibrium problem)
and phase_amount_indices (the indices of the parameters that
correspond to phase amounts).
:type prm: namedtuple
:returns: Equality constraints in a form that can be processed
by the F and jacobian functions.
:rtype: list of lists
"""
eq_constraint_lists = []
for i in range(len(equality_constraints)):
if equality_constraints[i][0] == "phase_fraction":
phase = equality_constraints[i][1][0]
fraction = equality_constraints[i][1][1]
if isinstance(fraction, float):
fraction = np.array([fraction])
if not isinstance(fraction, np.ndarray):
raise Exception(
"The constraint fraction in equality {0} "
"should be a float or numpy array".format(i + 1)
)
eq_constraint_lists.append(
phase_fraction_constraints(phase, assemblage, fraction, prm)
)
elif equality_constraints[i][0] == "phase_composition":
phase = equality_constraints[i][1][0]
constraint = equality_constraints[i][1][1]
if isinstance(constraint[3], float):
constraint = (
constraint[0],
constraint[1],
constraint[2],
np.array([constraint[3]]),
)
if not isinstance(constraint[3], np.ndarray):
raise Exception(
"The last constraint parameter in equality "
"{0} should be a float "
"or numpy array".format(i + 1)
)
eq_constraint_lists.append(
phase_composition_constraints(phase, assemblage, constraint, prm)
)
elif equality_constraints[i][0] == "X":
constraint = equality_constraints[i][1]
if isinstance(constraint[-1], float):
constraint = (constraint[0], np.array([constraint[-1]]))
if not isinstance(constraint[-1], np.ndarray):
raise Exception(
"The last constraint parameter in "
"equality {0} should be "
"a float or numpy array".format(i + 1)
)
eq_constraint_lists.append(
[["X", [np.array(constraint[0]), p]] for p in constraint[1]]
)
elif (
equality_constraints[i][0] == "P"
or equality_constraints[i][0] == "T"
or equality_constraints[i][0] == "S"
or equality_constraints[i][0] == "V"
):
if isinstance(equality_constraints[i][1], float):
equality_constraints[i] = (
equality_constraints[i][0],
np.array([equality_constraints[i][1]]),
)
if not isinstance(equality_constraints[i][1], np.ndarray):
raise Exception(
"The last parameter in "
f"equality_constraint[{i+1}] should be a "
"float or numpy array"
)
eq_constraint_lists.append(
[[equality_constraints[i][0], p] for p in equality_constraints[i][1]]
)
else:
raise Exception(
"The type of equality_constraint is "
"not recognised for constraint {0}.\n"
"Should be one of P, T, S, V, X,\n"
"phase_fraction, or phase_composition.".format(i + 1)
)
return eq_constraint_lists
[docs]
def equilibrate(
composition,
assemblage,
equality_constraints,
free_compositional_vectors=[],
tol=None,
store_iterates=False,
store_assemblage=True,
max_iterations=100.0,
verbose=False,
):
"""
A function that finds the thermodynamic equilibrium state of an
assemblage subject to given equality constraints by solving a set of
nonlinear equations related to the chemical potentials and
other state variables of the system.
The user chooses an assemblage (e.g. olivine, garnet and orthopyroxene)
and :math:`2+n_c` equality constraints, where :math:`n_c` is the number of
bulk compositional degrees of freedom. The equilibrate function attempts to
find the remaining unknowns that satisfy those constraints.
There are a number of equality constraints implemented in burnman. These are
given as a list of lists. Each constraint should have the form:
[<constraint type>, <constraint>], where
<constraint type> is one of 'P', 'T', 'S', 'V', 'X',
'phase_fraction', or 'phase_composition'. The format of the
<constraint> object depends on the constraint type:
- P: float or numpy.array of
pressures [Pa]
- T: float or numpy.array of
temperatures [K]
- S: float or numpy.array of
entropies [J/K]
- V: float or numpy.array of
volumes [m:math:`^3`]
- phase_fraction: tuple with the form (<phase> <fraction(s)>),
where <phase> is one of the phase objects in the assemblage
and <fraction(s)> is a float or numpy.array corresponding
to the desired phase fractions.
- phase_composition: tuple with the form (<phase> <constraint>),
where <phase> is one of the phase objects in the assemblage
and <constraint> has the form (site_names, n, d, v),
where :math:`(nx)/(dx) = v`, n and d are constant vectors of
site coefficients, and v is a float or numpy.array. For example,
a constraint of the form ([Mg_A, Fe_A], [1., 0.], [1., 1.], [0.5])
would correspond to equal amounts Mg and Fe on the A site
(Mg_A / (Mg_A + Fe_A) = 0.5).
- X: list of a numpy.array and a float or numpy.array,
where the equality satisfies the linear equation arr[0] x = arr[1].
The first numpy.array is fixed, and the second is to be looped over
by the equilibrate function.
This is a generic compositional equality constraint.
:param composition: The bulk composition that the assemblage must satisfy.
:type composition: dict
:param assemblage: The assemblage to be equilibrated.
:type assemblage: :class:`burnman.Composite` or :class:`burnman.RelaxedComposite`
:param equality_constraints: The list of equality constraints. See above
for valid formats.
:type equality_constraints: list of list
:param free_compositional_vectors: A list of dictionaries containing
the compositional freedom of the solution. For example, if the list
contains the vector {'Mg': 1., 'Fe': -1}, that implies that the bulk
composition is equal to composition + :math:`a` (n_Mg - n_Fe),
where the value of :math:`a` is to be determined by the solve.
Vector given in atomic (molar) units of elements.
:type free_compositional_vectors: list of dict
:param tol: The tolerance for the nonlinear solver. This can be a single float,
which is then applied to all parameters, or a numpy array with the
same length as the number of parameters in the solve (2 + number of
phases + number of free compositional vectors). The default is None,
which sets the tolerances to 1 Pa for pressure, 1e-3 K for temperature,
1e-6 for phase amounts and 1e-6 for free compositional vectors.
One of the termination criteria for the nonlinear solver is that
the absolute change in each parameter is less than the corresponding
tolerance.
:type tol: float or numpy.array
:param store_iterates: Whether to store the parameter values for
each iteration in each solution object.
:type store_iterates: bool
:param store_assemblage: Whether to store a copy of the assemblage
object in each solution object.
:type store_assemblage: bool
:param max_iterations: The maximum number of iterations for the
nonlinear solver.
:type max_iterations: int
:param verbose: Whether to print output updating the user on the status of
equilibration.
:type verbose: bool
:returns: Solver solution object (or a list, or a 2D list of solution objects)
created by :func:`burnman.optimize.nonlinear_solvers.damped_newton_solve`,
and a namedtuple object created by
:func:`burnman.tools.equilibration.get_equilibration_parameters`. See
documentation of these functions for the return types. If store_assemblage
is True, each solution object also has an attribute called `assemblage`,
which contains a copy of the input assemblage with the equilibrated
properties. So, for a 2D grid of solution objects, one could call either
sols[0][1].x[0] or sols[0][1].assemblage.pressure to get the pressure.
:rtype: tuple
"""
logger = logging.getLogger(__name__)
outer_logger_level = logger.level
logger.setLevel(logging.INFO if verbose else outer_logger_level)
relaxed_assemblage = None
if hasattr(assemblage, "unrelaxed"):
relaxed_assemblage = assemblage
assemblage = assemblage.unrelaxed
for ph in assemblage.phases:
if isinstance(ph, Solution) and not hasattr(ph, "molar_fractions"):
raise Exception(
f"set_composition for solution {ph} before running equilibrate."
)
if assemblage.molar_fractions is None:
n_phases = len(assemblage.phases)
f = 1.0 / float(n_phases)
assemblage.set_fractions([f for i in range(n_phases)])
n_atoms = sum(composition.values())
assemblage.number_of_moles = n_atoms / sum(assemblage.formula.values())
n_equality_constraints = len(equality_constraints)
n_free_compositional_vectors = len(free_compositional_vectors)
if n_equality_constraints != n_free_compositional_vectors + 2:
raise Exception(
"The number of equality constraints "
f"(currently {n_equality_constraints}) "
"must be two more than the number of "
"free_compositional vectors "
f"(currently {n_free_compositional_vectors})."
)
for v in free_compositional_vectors:
if np.abs(sum(v.values())) > 1.0e-12:
raise Exception(
"The amounts of each free_compositional_vector" "must sum to zero"
)
# Make parameter tuple
prm = get_equilibration_parameters(
assemblage, composition, free_compositional_vectors
)
# Set default tolerances if not provided
if tol is None:
tol = prm.default_tolerances
# Check equality constraints have the correct structure
# Convert into the format readable by the function and jacobian functions
eq_constraint_lists = process_eq_constraints(equality_constraints, assemblage, prm)
# Set up solves
nc = [len(eq_constraint_list) for eq_constraint_list in eq_constraint_lists]
try:
initial_state = [assemblage.pressure, assemblage.temperature]
except AttributeError:
initial_state = [None, None]
# Reset initial state if equality constraints
# are related to pressure or temperature
for i in range(n_equality_constraints):
if eq_constraint_lists[i][0][0] == "P":
initial_state[0] = eq_constraint_lists[i][0][1]
elif eq_constraint_lists[i][0][0] == "T":
initial_state[1] = eq_constraint_lists[i][0][1]
if initial_state[0] is None:
initial_state[0] = 5.0e9
if initial_state[1] is None:
initial_state[1] = 1200.0
assemblage.set_state(*initial_state)
parameters = get_parameters(assemblage, n_free_compositional_vectors)
parameters = np.array(parameters)
# Solve the system of equations, loop over input parameters
sol_array = np.empty(shape=tuple(nc), dtype="object")
# Loop over problems
problems = list(product(*[list(range(nc[i])) for i in range(len(nc))]))
n_problems = len(problems)
for i_problem, i_c in enumerate(problems):
if logger.isEnabledFor(logging.DEBUG):
string = "Processing solution"
for i in range(len(i_c)):
string += " {0}/{1}".format(i_c[i] + 1, nc[i])
logging.info(string + ":")
equality_constraints = [eq_constraint_lists[i][i_c[i]] for i in range(len(nc))]
F_tol = default_F_tolerances(assemblage, equality_constraints, n_atoms)
def F_fn(x):
return F(
x,
assemblage,
equality_constraints,
prm.reduced_composition_vector,
prm.reduced_free_composition_vectors,
)
def J_fn(x):
return jacobian(
x,
assemblage,
equality_constraints,
prm.reduced_free_composition_vectors,
)
def lambda_bounds_fn(dx, x):
return lambda_bounds(dx, x, assemblage.endmembers_per_phase)
linear_constraints = (prm.constraint_matrix, prm.constraint_vector)
sol = damped_newton_solve(
F=F_fn,
J=J_fn,
lambda_bounds=lambda_bounds_fn,
guess=parameters,
linear_constraints=linear_constraints,
tol=tol,
F_tol=F_tol,
store_iterates=store_iterates,
max_iterations=max_iterations,
)
# If we hit max iterations, we try relaxing the logish_eps
# and solving again. This helps convergence in hard problems because
# it allows larger steps when phase compositions are close to the
# boundaries of the solution (where the gradient in excess entropy is large).
if sol.code == TerminationCode.MAX_ITERATIONS:
old_logish_eps = logish_eps[0]
logish_eps[0] = 1.0e-5
c_shifts = sol.x[-n_free_compositional_vectors:]
new_parameters = get_parameters(assemblage, n_free_compositional_vectors)
new_parameters[-n_free_compositional_vectors:] = c_shifts
set_compositions_and_state_from_parameters(assemblage, parameters) # reset
sol2 = damped_newton_solve(
F=F_fn,
J=J_fn,
lambda_bounds=lambda_bounds_fn,
guess=new_parameters,
linear_constraints=linear_constraints,
tol=tol,
F_tol=F_tol,
store_iterates=store_iterates,
max_iterations=max_iterations,
)
logish_eps[0] = old_logish_eps
c_shifts = sol2.x[-n_free_compositional_vectors:]
new_parameters = get_parameters(assemblage, n_free_compositional_vectors)
new_parameters[-n_free_compositional_vectors:] = c_shifts
set_compositions_and_state_from_parameters(assemblage, parameters) # reset
sol3 = damped_newton_solve(
F=F_fn,
J=J_fn,
lambda_bounds=lambda_bounds_fn,
guess=new_parameters,
linear_constraints=linear_constraints,
tol=tol,
F_tol=F_tol,
store_iterates=store_iterates,
max_iterations=max_iterations,
)
# combine iterates from all three solves
sol3.n_it = sol.n_it + sol2.n_it + sol3.n_it
if store_iterates:
sol3.iterates.x = np.concatenate(
(sol.iterates.x, sol2.iterates.x, sol3.iterates.x), axis=0
)
sol3.iterates.F = np.concatenate(
(sol.iterates.F, sol2.iterates.F, sol3.iterates.F), axis=0
)
sol3.iterates.lmda = np.concatenate(
(sol.iterates.lmda, sol2.iterates.lmda, sol3.iterates.lmda),
axis=0,
)
sol = sol3
if relaxed_assemblage is not None:
relaxed_assemblage.copy_state_from_unrelaxed()
assemblage = relaxed_assemblage
if sol.success and len(assemblage.reaction_affinities) > 0:
maxres = np.max(np.abs(assemblage.reaction_affinities)) + 1.0e-5
assemblage.equilibrium_tolerance = maxres
if store_assemblage:
sol.assemblage = assemblage.copy()
if sol.success and len(assemblage.reaction_affinities) > 0.0:
sol.assemblage.equilibrium_tolerance = maxres
logging.info(sol.text)
# Scale solution back to physical units
sol.x[0] = sol.x[0] * P_scaling
sol.x[1] = sol.x[1] * T_scaling
if store_iterates:
sol.iterates.x[:, 0] = sol.iterates.x[:, 0] * P_scaling
sol.iterates.x[:, 1] = sol.iterates.x[:, 1] * T_scaling
sol_array[i_c] = sol
# Next, we use the solution values and Jacobian
# to provide a starting guess for the next problem.
# First, we find the equality constraints for the next problem
if i_problem < n_problems - 1:
next_i_c = problems[i_problem + 1]
next_equality_constraints = [
eq_constraint_lists[i][next_i_c[i]] for i in range(len(nc))
]
# We use the nearest solutions as potential starting points
# to make the next guess
prev_sols = []
for i in range(len(nc)):
if next_i_c[i] != 0:
prev_i_c = np.copy(next_i_c)
prev_i_c[i] -= 1
prev_sols.append(sol_array[tuple(prev_i_c)])
updated_params = False
for s in prev_sols:
if s.success and not updated_params:
# next guess based on a Newton step
# using the old solution vector and Jacobian
# with the new constraints.
x = s.x.copy()
x[0] = x[0] / P_scaling
x[1] = x[1] / T_scaling
dF = F(
x,
assemblage,
next_equality_constraints,
prm.reduced_composition_vector,
prm.reduced_free_composition_vectors,
)
luJ = lu_factor(s.J)
new_parameters = x + lu_solve(luJ, -dF)
c = (
prm.constraint_matrix.dot(new_parameters)
+ prm.constraint_vector
)
if all(c <= 0.0): # accept new guess
parameters = new_parameters
else: # use the parameters from this step
parameters = x
exhausted_phases = [
assemblage.phases[phase_idx].name
for phase_idx, v in enumerate(
new_parameters[prm.phase_amount_indices]
)
if v < 0.0
]
if len(exhausted_phases) > 0:
logging.info(
"A phase might be exhausted before the "
f"next step: {exhausted_phases}"
)
updated_params = True
# Finally, make dimensions of sol_array equal the input dimensions
if np.prod(sol_array.shape) > 1:
sol_array = np.squeeze(sol_array)
else:
sol_array = sol_array.flatten()[0]
logger.setLevel(outer_logger_level)
return sol_array, prm