quit solve on invalid step when J=singular

Author: bobmyhill

burnman/optimize/nonlinear_solvers.py

@@ -708,10 +708,16 @@ def solve(self) -> Solution:
             and not converged
         ):
             sol.J = self.J(sol.x)
+            condition_number = np.linalg.cond(sol.J)
             luJ = lu_factor(sol.J)
             dx = lu_solve(luJ, -sol.F)
             dx_norm = np.linalg.norm(dx, ord=2)
 
+            is_singular = condition_number > self.max_condition_number
+            if is_singular and any(np.isnan(dx)):
+                lmda_bounds = [0.0, 0.0]
+                break
+
             lmda_bounds = self.lambda_bounds(dx, sol.x)
             h = (
                 lmda
@@ -742,7 +748,6 @@ def solve(self) -> Solution:
             F_j = self.F(x_j)
 
             # Regularise ill-conditioned Jacobian
-            condition_number = np.linalg.cond(sol.J)
             if condition_number < self.cond_lu_thresh:
                 dxbar_j = lu_solve(luJ, -F_j)
             else:

burnman/tools/equilibration.py

@@ -498,10 +498,12 @@ def lambda_bounds(dx, x, endmembers_per_phase):
         ]
     )
     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}.")
+        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)