diff --git a/burnman/optimize/nonlinear_solvers.py b/burnman/optimize/nonlinear_solvers.py index 8e1a5e8b..1a5290b8 100644 --- a/burnman/optimize/nonlinear_solvers.py +++ b/burnman/optimize/nonlinear_solvers.py @@ -3,8 +3,82 @@ import numpy.typing as npt from typing import Any from scipy.linalg import lu_factor, lu_solve -from types import SimpleNamespace from collections import namedtuple +from types import SimpleNamespace + + +class Solution: + """ + Container for the results of a nonlinear solver. + Attributes: + - x: Final solution vector. + - n_it: Number of iterations performed. + - F: Function evaluation at the solution, F(x). + - F_norm: Euclidean (L2) norm of F(x). + - J: Jacobian matrix at the solution, J(x). + - code: Exit status code (0=success, >0=failure). + - text: Description of the exit status. + - success: Boolean indicating whether the solver converged. + - iterates: Object storing all intermediate iterates, function evaluations, + and lambda values (if enabled). + """ + + def __init__( + self, + x=None, + n_it=None, + F=None, + F_norm=None, + J=None, + code=None, + text=None, + success=False, + iterates=None, + ): + self.x = x + self.n_it = n_it + self.F = F + self.F_norm = F_norm + self.J = J + self.code = code + self.text = text + self.success = success + + if iterates is None: + self.iterates = Iterates() + else: + self.iterates = iterates + + def __repr__(self): + s = f"{self.text}\n" + s += f"x = {self.x}\n" + s += f"F = {self.F}\n" + s += f"F_norm = {self.F_norm}\n" + s += f"J = {self.J}\n" + if len(self.iterates.x) > 0: + s += f"{self.iterates}" + return s + + def __str__(self): + return self.__repr__() + + +class Iterates: + def __init__(self): + self.x = [] + self.F = [] + self.lmda = [] + + def append(self, x, F, lmda): + self.x.append(x) + self.F.append(F) + self.lmda.append(lmda) + + def __repr__(self): + s = "Iterates:\n" + for i in range(len(self.x)): + s += f"Iter {i}: x={self.x[i]}, F={self.F[i]}, lmda={self.lmda[i]}\n" + return s class DampedNewtonSolver: @@ -46,7 +120,8 @@ def __init__( F, J, guess, - tol: float = 1.0e-6, + tol: float | npt.NDArray = 1.0e-6, + F_tol: float | npt.NDArray = 1.0e-8, max_iterations: int = 100, lambda_bounds=lambda dx, x: (1.0e-8, 1.0), linear_constraints=(0.0, np.array([-1.0])), @@ -69,8 +144,13 @@ def __init__( :param guess: Initial guess for the solution vector x. :type guess: np.ndarray - :param tol: Convergence tolerance for Newton iterations, defaults to 1.0e-6. - :type tol: float, optional + :param tol: Convergence tolerance for each component of the + Newton step dx, defaults to 1.0e-6. + :type tol: float, or numpy.ndarray, optional + + :param F_tol: Convergence tolerance for each component of F(x), + defaults to 1.0e-8. + :type F_tol: float, or numpy.ndarray, optional :param max_iterations: Maximum number of Newton iterations to perform, defaults to 100. @@ -105,6 +185,7 @@ def __init__( self.J = J self.guess = guess self.tol = tol + self.F_tol = F_tol self.max_iterations = max_iterations self.lambda_bounds = lambda_bounds self.linear_constraints = linear_constraints @@ -533,7 +614,7 @@ def _termination_info( return False, 3, f"The solver reached max_iterations ({max_iterations})" raise Exception("Unknown termination of solver") - def solve(self) -> SimpleNamespace: + def solve(self) -> Solution: """ Execute the damped Newton solver to find a root of F(x) = 0, optionally subject to linear inequality constraints A·x + b <= 0. @@ -564,17 +645,14 @@ def solve(self) -> SimpleNamespace: * **F** (list[np.ndarray]) -- Function evaluation for each iteration. * **lmda** (list[float]) -- Step length scaling factor for each iteration. - :rtype: SimpleNamespace + :rtype: Solution instance """ - sol = namedtuple( - "Solution", ["x", "n_it", "F", "F_norm", "J", "code", "text", "success"] - ) + sol = Solution() sol.x = self.guess sol.F = self.F(sol.x) if self.store_iterates: - sol.iterates = namedtuple("iterates", ["x", "F", "lmda"]) - sol.iterates.x, sol.iterates.F, sol.iterates.lmda = [sol.x], [sol.F], [0.0] + sol.iterates.append(sol.x, sol.F, 0.0) lmda, dxprev, dxbar = 0.0, [1.0], [1.0] sol.n_it = 0 @@ -632,6 +710,11 @@ def solve(self) -> SimpleNamespace: dxbar_j_norm = np.linalg.norm(dxbar_j, ord=2) converged = self._check_convergence(dxbar_j, dx, lmda, lmda_bounds) + + # Additional convergence check on F(x) + if converged and not all(np.abs(F_j) < self.F_tol): + converged = False + require_posteriori_loop = not converged loop_vars = self._posteriori_loop( @@ -655,9 +738,7 @@ def solve(self) -> SimpleNamespace: sol.n_it += 1 if self.store_iterates: - sol.iterates.x.append(sol.x) - sol.iterates.F.append(sol.F) - sol.iterates.lmda.append(lmda) + sol.iterates.append(sol.x, sol.F, lmda) # Final adjustment for constraints if converged and not persistent_bound_violation: @@ -669,10 +750,6 @@ def solve(self) -> SimpleNamespace: sol.F_norm = np.linalg.norm(sol.F, ord=2) sol.J = self.J(sol.x) - if self.store_iterates: - sol.iterates.x = np.array(sol.iterates.x) - sol.iterates.F = np.array(sol.iterates.F) - sol.success, sol.code, sol.text = self._termination_info( converged, minimum_lmda, @@ -689,7 +766,8 @@ def damped_newton_solve( F, J, guess, - tol: float = 1.0e-6, + tol: float | npt.NDArray = 1.0e-6, + F_tol: float | npt.NDArray = 1.0e-8, max_iterations: int = 100, lambda_bounds=lambda dx, x: (1.0e-8, 1.0), linear_constraints=(0.0, np.array([-1.0])), @@ -737,6 +815,7 @@ def damped_newton_solve( J, guess, tol, + F_tol, max_iterations, lambda_bounds, linear_constraints, diff --git a/burnman/tools/equilibration.py b/burnman/tools/equilibration.py index 60d51025..802cf0fa 100644 --- a/burnman/tools/equilibration.py +++ b/burnman/tools/equilibration.py @@ -169,6 +169,38 @@ def set_compositions_and_state_from_parameters(assemblage, parameters): return None +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 # Pa + elif type_c == "T": + F_tolerances[i] = 1.0e-6 # K + elif type_c == "S": + F_tolerances[i] = 1.0e-8 # 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 # m^3 + elif type_c == "PT_ellipse": + F_tolerances[i] = 1.0e-3 # [no units] + elif type_c == "X": + F_tolerances[i] = 1.0e-8 # [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 + F_tolerances[i + assemblage.n_reactions :] = 1.0e-8 + return F_tolerances + + def F( x, assemblage, @@ -320,7 +352,7 @@ def jacobian(x, assemblage, equality_constraints, reduced_free_composition_vecto -assemblage.n_moles * assemblage.molar_volume / assemblage.isothermal_bulk_modulus_reuss, - assemblage.n_moles * assemblage.molar_volume, + assemblage.n_moles * assemblage.molar_volume * assemblage.alpha, ] j = 2 for k, n in enumerate(assemblage.endmembers_per_phase): @@ -934,7 +966,8 @@ def equilibrate( f = 1.0 / float(n_phases) assemblage.set_fractions([f for i in range(n_phases)]) - assemblage.n_moles = sum(composition.values()) / sum(assemblage.formula.values()) + n_atoms = sum(composition.values()) + assemblage.n_moles = n_atoms / sum(assemblage.formula.values()) n_equality_constraints = len(equality_constraints) n_free_compositional_vectors = len(free_compositional_vectors) @@ -1007,6 +1040,8 @@ def equilibrate( 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) + # Set the initial fractions and compositions # of the phases in the assemblage: sol = damped_newton_solve( @@ -1029,6 +1064,7 @@ def equilibrate( guess=parameters, linear_constraints=(prm.constraint_matrix, prm.constraint_vector), tol=tol, + F_tol=F_tol, store_iterates=store_iterates, max_iterations=max_iterations, ) @@ -1036,10 +1072,6 @@ def equilibrate( if sol.success and len(assemblage.reaction_affinities) > 0.0: maxres = np.max(np.abs(assemblage.reaction_affinities)) + 1.0e-5 assemblage.equilibrium_tolerance = maxres - assert maxres < 10.0, ( - "Equilibrium not satisfied to high enough precision. " - f"Max reaction affinity is {maxres} J/mol. Try reducing tol." - ) if store_assemblage: sol.assemblage = assemblage.copy() diff --git a/burnman/tools/polytope.py b/burnman/tools/polytope.py index 328a8beb..d53c7dc6 100644 --- a/burnman/tools/polytope.py +++ b/burnman/tools/polytope.py @@ -174,6 +174,19 @@ def composite_polytope_at_constrained_composition( return MaterialPolytope(equalities, inequalities, return_fractions=return_fractions) +def reorder_to_echelon(A): + """ + Reorders the rows of a matrix A so that it is in echelon form. + :param A: The input matrix. + :type A: 2D numpy array + :returns: The reordered matrix. + :rtype: 2D numpy array + """ + A = np.array(A, dtype=float) + pivot_cols = [np.argmax(r != 0) if np.any(r != 0) else np.inf for r in A] + return A[np.argsort(pivot_cols)] + + def simplify_composite_with_composition(composite, composition): """ Takes a composite and a composition, and returns the simplest composite @@ -249,6 +262,11 @@ def simplify_composite_with_composition(composite, composition): composite_changed = True + # Try to preserve the order of endmembers + # in the original solution model + # by reordering rows of the new_basis matrix + new_basis = reorder_to_echelon(new_basis) + # If the composition is set and # consistent with the new basis, # make a new solution with the composition diff --git a/misc/benchmarks/slb_2024_benchmarks.py b/misc/benchmarks/slb_2024_benchmarks.py index e0bdc267..789545af 100644 --- a/misc/benchmarks/slb_2024_benchmarks.py +++ b/misc/benchmarks/slb_2024_benchmarks.py @@ -451,19 +451,21 @@ def check_fig_7_fO2(): ax[0].imshow(fig1, extent=[0, 100.0, -14.0, 22.0], aspect="auto") for i, assemblage in enumerate([a1, a2, a3, a4, a5, a6, a7, a8, a9]): + wu.set_composition([0.8, 0.005, 0.195]) pressures = np.linspace(P_bounds[i][0], P_bounds[i][1], 11) logfO2 = np.empty_like(pressures) * np.nan T = 1500.0 for i, P in enumerate(pressures): assemblage.set_state(P, T) - sol, _ = equilibrate( - assemblage.formula, assemblage, [["P", P], ["T", T]], tol=1.0e-10 - ) + sol, _ = equilibrate(assemblage.formula, assemblage, [["P", P], ["T", T]]) if sol.success: mu_O2 = sol.assemblage.chemical_potential([{"O": 2.0}])[0] O2_gas.set_state(1.0e5, T) O2_gibbs = O2_gas.gibbs + f logfO2[i] = (mu_O2 - O2_gibbs) / (gas_constant * T) / np.log(10.0) + else: + print(assemblage) + print(sol) ax[0].plot(pressures / 1.0e9, logfO2, linestyle=":", linewidth=3.0) diff --git a/misc/ref/example_polytopetools.py.out b/misc/ref/example_polytopetools.py.out index d75ac92d..899b9786 100644 --- a/misc/ref/example_polytopetools.py.out +++ b/misc/ref/example_polytopetools.py.out @@ -31,7 +31,7 @@ The complete set of endmembers expressed as proportions of the independent endme [[ 1. -0. -0. -0.] [ 0. 1. -0. 0.] [-0. 0. 1. -0.] - [-0. 0. 0. 1.] + [-0. -0. -0. 1.] [-1. 2. 2. -2.]] We can also create a polytope using the site-species occupancies of a set of independent endmembers of a solution. Here we do this for the garnet solid solution provided in the Stixrude and Lithgow-Bertelloni database. @@ -42,13 +42,13 @@ Ca3Al2Si3O12 Mg4Si4O12 Na2Al2Si4O12 Endmembers as independent endmember amounts: -[[-0. -0. -0. 1. 0.] +[[ 0. -0. 0. 1. -0.] [ 1. 0. 0. 0. 0.] - [ 0. 1. 0. -0. -0.] + [-0. 1. -0. 0. -0.] [-1. 1. 0. 1. -0.] - [ 0. 0. 1. -0. 0.] - [-1. -0. 1. 1. 0.] - [ 0. -0. -0. 0. 1.]] + [ 0. -0. 1. 0. 0.] + [-1. -0. 1. 1. -0.] + [ 0. -0. 0. -0. 1.]] Site occupancies for all endmembers: [Mg]3[Mg][Si] [Mg]3[Al][Al] @@ -72,7 +72,7 @@ Which led to the following extreme vertices satifying the bulk compositional con [0 0 0 5/2 0 4 1]] The new assemblage has fewer endmembers, with compositions: -['Mg4Si4O12', 'MgFe3Si4O12', 'Mg2SiO4', 'Fe2SiO4'] +['MgFe3Si4O12', 'Mg4Si4O12', 'Mg2SiO4', 'Fe2SiO4'] The same vertices now look like this: -[[11/6 2/3 5 0] - [5/2 0 4 1]] +[[2/3 11/6 5 0] + [0 5/2 4 1]] diff --git a/tests/test_equilibration.py b/tests/test_equilibration.py index c42397fc..f068ef3c 100644 --- a/tests/test_equilibration.py +++ b/tests/test_equilibration.py @@ -8,6 +8,16 @@ import numpy as np +def make_ol_wad_assemblage(): + ol = SLB_2011.mg_fe_olivine() + wad = SLB_2011.mg_fe_wadsleyite() + + assemblage = burnman.Composite([ol, wad], [0.7, 0.3]) + ol.set_composition([0.5, 0.5]) + wad.set_composition([0.6, 0.4]) + return assemblage + + class equilibration(BurnManTest): def test_univariant_line(self): andalusite = HP_2011_ds62.andalusite() @@ -89,13 +99,9 @@ def test_univariant_line_fo(self): self.assertArraysAlmostEqual(Ps, Ps_ref) def test_ol_wad_eqm(self): - ol = SLB_2011.mg_fe_olivine() - wad = SLB_2011.mg_fe_wadsleyite() - - assemblage = burnman.Composite([ol, wad], [0.7, 0.3]) - ol.set_composition([0.5, 0.5]) - wad.set_composition([0.6, 0.4]) - + assemblage = make_ol_wad_assemblage() + ol = assemblage.phases[0] + wad = assemblage.phases[1] assemblage.set_state(10.0e9, 1200.0) equality_constraints = [ ("P", 10.0e9), @@ -145,13 +151,8 @@ def test_binary_solution_convergence(self): self.assertFalse(any(not s.success for s in sol)) def test_incorrect_tol(self): - ol = SLB_2011.mg_fe_olivine() - wad = SLB_2011.mg_fe_wadsleyite() - - assemblage = burnman.Composite([ol, wad], [0.7, 0.3]) - ol.set_composition([0.5, 0.5]) - wad.set_composition([0.6, 0.4]) - + assemblage = make_ol_wad_assemblage() + ol = assemblage.phases[0] assemblage.set_state(10.0e9, 1200.0) equality_constraints = [ ("P", 10.0e9), @@ -182,6 +183,67 @@ def test_incorrect_tol(self): ) _ = equilibrate(composition, assemblage, equality_constraints) + def test_ol_wad_eqm_entropy(self): + assemblage = make_ol_wad_assemblage() + wad = assemblage.phases[1] + equality_constraints = [ + ("P", 10.0e9), + ("phase_fraction", (wad, 0.5)), + ] + composition = {"Mg": 1.0, "Fe": 1.0, "Si": 1.0, "O": 4.0} + + sol, _ = equilibrate(composition, assemblage, equality_constraints) + self.assertTrue(sol.success) + + P = assemblage.pressure + T = assemblage.temperature + S = assemblage.molar_entropy * assemblage.n_moles + + assemblage = make_ol_wad_assemblage() + equality_constraints = [("P", P), ("S", S)] + sol, _ = equilibrate(composition, assemblage, equality_constraints) + self.assertTrue(sol.success) + self.assertAlmostEqual(assemblage.temperature, T, places=5) + + def test_ol_wad_eqm_volume(self): + assemblage = make_ol_wad_assemblage() + wad = assemblage.phases[1] + equality_constraints = [ + ("P", 10.0e9), + ("phase_fraction", (wad, 0.5)), + ] + composition = {"Mg": 1.0, "Fe": 1.0, "Si": 1.0, "O": 4.0} + + sol, _ = equilibrate(composition, assemblage, equality_constraints) + self.assertTrue(sol.success) + + P = assemblage.pressure + T = assemblage.temperature + V = assemblage.molar_volume * assemblage.n_moles + + assemblage = make_ol_wad_assemblage() + equality_constraints = [("P", P), ("V", V)] + sol, _ = equilibrate(composition, assemblage, equality_constraints) + self.assertTrue(sol.success) + self.assertAlmostEqual(assemblage.temperature, T, places=5) + + def test_ol_wad_eqm_compositional_constraint(self): + assemblage = make_ol_wad_assemblage() + ol = assemblage.phases[0] + wad = assemblage.phases[1] + assemblage.set_state(10.0e9, 1200.0) + equality_constraints = [ + ("P", 10.0e9), + ("X", [[0.0, 0.0, 0.0, 2.0, 0.0, 0.0], 0.9]), + ] + composition = {"Mg": 1.0, "Fe": 1.0, "Si": 1.0, "O": 4.0} + + _, _ = equilibrate(composition, assemblage, equality_constraints) + self.assertArraysAlmostEqual( + [assemblage.temperature, ol.molar_fractions[1], wad.molar_fractions[1]], + [1620.532183457096, 0.45, 0.6791743], + ) + if __name__ == "__main__": unittest.main() diff --git a/tests/test_polytope.py b/tests/test_polytope.py index f7b1039c..91ddb4e2 100644 --- a/tests/test_polytope.py +++ b/tests/test_polytope.py @@ -6,7 +6,7 @@ import importlib from burnman import Composite -from burnman.minerals import SLB_2011, JH_2015 +from burnman.minerals import SLB_2011, SLB_2024, JH_2015 from burnman.tools.polytope import solution_polytope_from_charge_balance from burnman.tools.polytope import solution_polytope_from_endmember_occupancies from burnman.tools.polytope import simplify_composite_with_composition @@ -135,6 +135,15 @@ def test_cddlib_versions(self): except ImportError: self.assertTrue(type(poly.raw_vertices[0][0]) is float) + def test_simplified_composite_endmember_order(self): + wu = SLB_2024.ferropericlase() + wu.set_composition([0.0, 0.97, 0.005, 0.0, 0.025]) + wu = simplify_composite_with_composition( + Composite([wu]), {"Fe": 1, "O": 1.01} + ).phases[0] + self.assertTrue(len(wu.endmembers) == 3) + self.assertArraysAlmostEqual(wu.molar_fractions, [0.97, 0.005, 0.025]) + if __name__ == "__main__": unittest.main()