Change Incompatible constraints to warn, add error if too large of a difference

desc/objectives/utils.py

@@ -149,15 +149,42 @@
         y1 = con.compute_unscaled(*xpi)
         y2 = con.target
         y1, y2 = np.broadcast_arrays(y1, y2)
+
+        # If the error is very large, likely want to error out as
+        # it probably is due to a real mistake instead of just numerical
+        # roundoff errors.
         np.testing.assert_allclose(
             y1,
             y2,
-            atol=3e-14,
-            rtol=3e-14,
+            atol=1e-6,
+            rtol=1e-1,
             err_msg="Incompatible constraints detected, cannot satisfy constraint "
             + f"{con}.",
         )
 
+        # else check with tighter tols and throw an error, these tolerances
+        # could be tripped due to just numerical roundoff or poor scaling between
+        # constraints, so don't want to error out but we do want to warn the user.
+        atol = 3e-14
+        rtol = 3e-14
+
+        try:
+            np.testing.assert_allclose(
+                y1,
+                y2,
+                atol=atol,
+                rtol=rtol,
+                err_msg="Incompatible constraints detected, cannot satisfy constraint "
+                + f"{con}.",
+            )
+        except AssertionError as e:
+            warnif(
+                True,
+                UserWarning,
+                str(e) + "\n This may indicate incompatible constraints, "
+                "or be due to floating point error.",
+            )
+
     return xp, A, b, Z, unfixed_idx, project, recover
 
 

tests/test_linear_objectives.py

@@ -133,25 +133,37 @@ def test_constrain_bdry_with_only_one_mode():
 def test_constrain_asserts():
     """Test error checking for incompatible constraints."""
     eqi = Equilibrium(iota=PowerSeriesProfile(0, 0), pressure=PowerSeriesProfile(0, 0))
-    eqc = Equilibrium(current=PowerSeriesProfile(0))
+    eqc = Equilibrium(current=PowerSeriesProfile(1))
+    obj_i = get_equilibrium_objective(eqi, "force")
+    obj_c = get_equilibrium_objective(eqc, "force")
+    obj_i.build()
+    obj_c.build()
     # nonexistent toroidal current can't be constrained
     with pytest.raises(RuntimeError):
-        eqi.solve(constraints=FixCurrent(eq=eqi))
+        con = FixCurrent(eq=eqi)
+        con.build()
     # nonexistent rotational transform can't be constrained
     with pytest.raises(RuntimeError):
-        eqc.solve(constraints=FixIota(eq=eqc))
+        con = FixIota(eq=eqc)
+        con.build()
     # toroidal current and rotational transform can't be constrained simultaneously
     with pytest.raises(ValueError):
-        eqi.solve(constraints=(FixCurrent(eq=eqi), FixIota(eq=eqi)))
-    with pytest.raises(AssertionError):
-        eqi.solve(
-            constraints=(FixPressure(eq=eqi, target=2), FixPressure(eq=eqi, target=1))
-        )
+        con = (FixCurrent(eq=eqi), FixIota(eq=eqi))
+        eqi.solve(constraints=con)
     # cannot use two incompatible constraints
     with pytest.raises(AssertionError):
         con1 = FixCurrent(target=eqc.c_l, eq=eqc)
         con2 = FixCurrent(target=eqc.c_l + 1, eq=eqc)
-        eqc.solve(constraints=(con1, con2))
+        con = ObjectiveFunction((con1, con2))
+        con.build()
+        _ = factorize_linear_constraints(obj_c, con)
+    # if only slightly off, should raise only a warning
+    with pytest.warns(UserWarning):
+        con1 = FixCurrent(target=eqc.c_l, eq=eqc)
+        con2 = FixCurrent(target=eqc.c_l * (1 + 1e-9), eq=eqc)
+        con = ObjectiveFunction((con1, con2))
+        con.build()
+        _ = factorize_linear_constraints(obj_c, con)
 
 
 @pytest.mark.regression