Add F(x) termination criterion to DampedNewtonSolver, fix Jacobian for volume constraints

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:
@@ -47,6 +121,7 @@ def __init__(
         J,
         guess,
         tol: float = 1.0e-6,
+        F_tol: float = 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,
@@ -690,6 +767,7 @@ def damped_newton_solve(
     J,
     guess,
     tol: float = 1.0e-6,
+    F_tol: float = 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,

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,
@@ -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,17 +1064,14 @@ 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,
         )
 
         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()

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

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)
 

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]]

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