Address #85 review: fix convergence off-by-one + guard scalar-path Newton locals

Author: petlenz

include/numsim_codegen/code_emit/linear_algebra_emitter.h

@@ -59,12 +59,16 @@ class LinearAlgebraEmitter {
   // step sets `true` when `‖R‖∞ < tol`, so the render frame can detect a
   // max-iter/singular failure after the loop. The step also `break`s without
   // setting it when the dense solve is non-finite (singular Jacobian).
+  // `iter_var` / `max_iter` (issue #85): the loop runs `iter_var <= max_iter`,
+  // so the residual of the final update is checked; the step caps the UPDATE at
+  // `max_iter` by `break`ing (post-convergence-check) when `iter_var == max_iter`.
   virtual void emit_newton_step(
       std::ostream &os, std::string const &prefix,
       std::vector<std::string> const &unknowns,
       std::vector<std::string> const &residual_rhs,
       std::vector<std::vector<std::string>> const &jacobian_rhs,
-      double tol, std::string const &converged_flag) const = 0;
+      double tol, std::string const &converged_flag,
+      std::string const &iter_var, int max_iter) const = 0;
 };
 
 // Eigen implementation — fixed-size `Eigen::Matrix<double,N,·>` + partial-pivot
@@ -86,7 +90,8 @@ class EigenLinearAlgebraEmitter final : public LinearAlgebraEmitter {
       std::vector<std::string> const &unknowns,
       std::vector<std::string> const &residual_rhs,
       std::vector<std::vector<std::string>> const &jacobian_rhs,
-      double tol, std::string const &converged_flag) const override {
+      double tol, std::string const &converged_flag,
+      std::string const &iter_var, int max_iter) const override {
     auto const n = unknowns.size();
     auto L = [&](char const *s) { return prefix + "_" + s; };
     // Residual vector R (size N).
@@ -108,10 +113,16 @@ class EigenLinearAlgebraEmitter final : public LinearAlgebraEmitter {
     // Convergence on the residual ∞-norm (issue #85: record it).
     os << "    if (" << L("r") << ".cwiseAbs().maxCoeff() < " << tol
        << ") {\n      " << converged_flag << " = true;\n      break;\n    }\n";
+    // issue #85: the loop runs one pass past max_iter to check the final
+    // update's residual (above); cap the UPDATE itself at max_iter.
+    os << "    if (" << iter_var << " == " << max_iter
+       << ") {\n      break;\n    }\n";
     // Solve J·Δx = R. issue #85: partialPivLu does NO rank check — on a
     // singular J it returns garbage (often non-finite). Guard the solution and
-    // break rather than propagate it; a near-singular-but-finite step is caught
-    // by the outer convergence flag (the loop exhausts max_iter).
+    // break rather than propagate it. (A near-singular-but-finite solution is
+    // NOT caught here — it either fails to converge within max_iter, or, like any
+    // Newton-without-globalization, may reach a spurious root; inherent, not a
+    // guarantee this guard provides.)
     os << "    Eigen::Matrix<double, " << n << ", 1> " << L("dx") << " = "
        << L("J") << ".partialPivLu().solve(" << L("r") << ");\n";
     os << "    if (!" << L("dx") << ".allFinite()) {\n      break;\n    }\n";
@@ -128,6 +139,13 @@ class EigenLinearAlgebraEmitter final : public LinearAlgebraEmitter {
 // `cwiseAbs().maxCoeff()`, `arma::solve` vs `partialPivLu().solve`), yet it
 // slots behind the same four methods. NOTE: `arma::solve` links LAPACK/BLAS, a
 // heavier downstream dependency than header-only Eigen.
+//
+// STATUS: opt-in DEMONSTRATOR, kept for a possible future non-Eigen backend.
+// Nothing in the library selects it (Eigen is the default and the only backend
+// any target uses), so its emitted code is exercised by STRING assertions only
+// (LinearAlgebraEmitterTest) — it is never compiled or run. Treat the emitted
+// Armadillo (e.g. the bool-overload `arma::solve(dx, J, r)` below) as
+// plausible-but-unverified until a gated compile test is added.
 class ArmadilloLinearAlgebraEmitter final : public LinearAlgebraEmitter {
 public:
   [[nodiscard]] auto includes() const -> std::vector<std::string> override {
@@ -145,7 +163,8 @@ class ArmadilloLinearAlgebraEmitter final : public LinearAlgebraEmitter {
       std::vector<std::string> const &unknowns,
       std::vector<std::string> const &residual_rhs,
       std::vector<std::vector<std::string>> const &jacobian_rhs,
-      double tol, std::string const &converged_flag) const override {
+      double tol, std::string const &converged_flag,
+      std::string const &iter_var, int max_iter) const override {
     auto const n = unknowns.size();
     auto L = [&](char const *s) { return prefix + "_" + s; };
     // Residual vector R (size N) — element assignment (Armadillo has no
@@ -165,6 +184,9 @@ class ArmadilloLinearAlgebraEmitter final : public LinearAlgebraEmitter {
     // Convergence on the residual ∞-norm (issue #85: record it).
     os << "    if (arma::norm(" << L("r") << ", \"inf\") < " << tol
        << ") {\n      " << converged_flag << " = true;\n      break;\n    }\n";
+    // issue #85: cap the update at max_iter (see the Eigen sibling).
+    os << "    if (" << iter_var << " == " << max_iter
+       << ") {\n      break;\n    }\n";
     // Solve J·Δx = R via arma::solve's bool overload, which returns false on a
     // singular / rank-deficient system (issue #85) — break rather than update
     // with garbage. A near-singular-but-finite step is caught by the outer

include/numsim_codegen/recipe.h

@@ -255,23 +255,43 @@ struct UpdateEquation {
   std::string doc;
 };
 
-// Phase 3a-2 (issue #75): tuning for the in-function Newton solve emitted
-// by LocalNewtonLoweringPass. `tol` is the absolute residual threshold for
-// convergence; `max_iter` caps the iteration count. Recipe-level for now
+// Phase 3a-2 (issue #75): tuning for the in-function Newton solve emitted by
+// LocalNewtonLoweringPass — i.e. the SELF-CONTAINED (standalone / MOOSE-local)
+// target path. `tol` is the ABSOLUTE residual ∞-norm threshold; `max_iter` is
+// the maximum number of Newton UPDATES (the emitted loop additionally checks the
+// residual of that final update before deciding failure). Recipe-level for now
 // (uniform across all evolution equations); per-equation tuning can be added
 // when a real consumer needs it.
+//
+// SCOPE (issue #85): these options drive ONLY the self-contained emitted loop.
+// The graph-coupled NumSimMaterialTarget (Mode B) delegates to numsim-materials'
+// `backward_euler`, which has its own convergence handling and ignores
+// `on_failure` — a residual-equation recipe emitted there is unaffected by this
+// struct.
+//
+// `tol` is ABSOLUTE: a well-scaled solve whose residual carries large physical
+// units may never reach a tiny absolute `tol` and will then be reported as
+// FAILED (below). Set `tol` to the residual's scale; a relative test is a future
+// enhancement.
 struct NewtonOptions {
-  // Failure policy (issue #85): what the generated solve does when it fails to
-  // converge — the loop exhausts `max_iter` without the residual reaching `tol`,
-  // OR the Newton step is non-finite (singular Jacobian). Without this the loop
+  // Failure policy (issue #85): what the generated solve does when it fails —
+  // the loop exhausts `max_iter` updates without the residual reaching `tol`, OR
+  // the Newton step is non-finite (singular Jacobian). Without this the loop
   // silently wrote back the last (un-converged) iterate, corrupting the host FE
   // solve with no signal.
-  //   NaN   — write a quiet NaN to the failed state variable(s); it propagates
-  //           into every output derived from them, so the host detects a
-  //           non-finite result (MOOSE cutback-friendly — no throw in a QP loop).
+  //   NaN   — write a quiet NaN to the failed state variable(s). It propagates
+  //           through ordinary arithmetic into outputs derived from the state,
+  //           so a host detects a non-finite result (MOOSE cutback-friendly — no
+  //           throw in a QP loop). CONTRACT: the host must check the STATE
+  //           out-param(s), which are always NaN'd. NaN is NOT guaranteed to
+  //           survive a derived output that passes it through `std::fmax`/`fmin`
+  //           (IEEE min/max ignore a NaN operand — e.g. a Macaulay bracket /
+  //           `max(·,0)` in a return map), so an output alone can read finite.
+  //           Use `Throw` when the failure must be unmissable.
   //   Throw — throw std::runtime_error from the generated compute function
   //           (clean for standalone drivers; unsafe inside a MOOSE QP loop).
-  // NaN is the default: universally detectable and safe on every backend.
+  // NaN is the default: universally detectable on the state and safe on every
+  // backend.
   enum class OnFailure { NaN, Throw };
 
   double tol = 1e-10;
@@ -1542,6 +1562,42 @@ inline auto tensor_arg_count(RecipeView model) -> int {
   return n;
 }
 
+// Identifiers already live in the generated function's scope: every declared
+// symbol (inputs, params, each state var's current + `_old`) and every output
+// (`<out>` and its `<out>_out` param). A Newton solve's scaffolding locals must
+// avoid all of these. (issue #85 / PR #83.)
+inline auto newton_reserved_ids(RecipeView model) -> std::set<std::string> {
+  std::set<std::string> reserved;
+  for (auto const &[nm, _] : model.scalar_symbol_map())
+    reserved.insert(nm);
+  for (auto const &[nm, _] : model.tensor_symbol_map())
+    reserved.insert(nm);
+  for (auto const &o : model.outputs()) {
+    reserved.insert(o.name);
+    reserved.insert(o.name + "_out");
+  }
+  return reserved;
+}
+
+// Mangle `seed` (by appending `_`) until none of `seed + "_" + suffix` collides
+// with a reserved identifier — so the solve's scaffolding locals (`<p>_R`,
+// `<p>_iter`, `<p>_converged`, …) are guaranteed unique in the emitted scope.
+// Shared by the scalar-segment and coupled-system loops (issue #85 closes the
+// scalar-path gap: PR #83 had mangled only the coupled path).
+inline auto collision_free_prefix(std::string seed,
+                                  std::set<std::string> const &reserved,
+                                  std::vector<std::string> const &suffixes)
+    -> std::string {
+  auto collides = [&](std::string const &pre) {
+    for (auto const &s : suffixes)
+      if (reserved.contains(pre + "_" + s)) return true;
+    return false;
+  };
+  while (collides(seed))
+    seed += "_";
+  return seed;
+}
+
 // issue #85: emit the post-loop failure action for a local Newton solve. Runs
 // after the loop, keyed on the `converged_flag` the loop sets true only on a
 // genuine `‖R‖ < tol` exit. On failure (max-iter exhausted or a singular
@@ -1790,28 +1846,47 @@ inline auto render_compute_function(
 
   // Newton iteration segments first — they declare + solve the local
   // state-variable iterates that downstream output expressions reference.
+  //
+  // The iterate keeps the state-variable name `<sv>` (output expressions
+  // reference it), but the solve's scaffolding locals use a collision-free
+  // prefix `<p>` (issue #85: PR #83 mangled only the coupled path, leaving
+  // `<sv>_R`/`_J`/`_iter`/`_converged` able to redeclare a like-named symbol).
+  auto const scalar_reserved = detail::newton_reserved_ids(model);
+  std::vector<std::string> const scalar_suffixes{"R", "J", "d", "iter",
+                                                 "converged"};
   for (auto const &seg : newton) {
     auto const &sv = seg.state_var_name;
+    auto const p =
+        detail::collision_free_prefix(sv, scalar_reserved, scalar_suffixes);
     os << "  double " << sv << " = " << sv << "_old;\n";
     // issue #85: track convergence so a max-iter exhaustion or a singular
     // Jacobian is detected rather than silently writing the last iterate.
-    os << "  bool " << sv << "_converged = false;\n";
-    os << "  for (int " << sv << "_iter = 0; " << sv << "_iter < "
-       << seg.max_iter << "; ++" << sv << "_iter) {\n";
+    os << "  bool " << p << "_converged = false;\n";
+    // `<= max_iter` so the residual of the FINAL (max_iter-th) update is checked
+    // before the loop exits — a solve that converges on its last allowed update
+    // must be recognised, not falsely NaN'd. The update itself is still capped
+    // at max_iter (the `== max_iter` guard below), so `max_iter` is the update
+    // budget.
+    os << "  for (int " << p << "_iter = 0; " << p << "_iter <= "
+       << seg.max_iter << "; ++" << p << "_iter) {\n";
     os << seg.loop_local_decls; // already 4-indented, trailing newline
-    os << "    double " << sv << "_R = " << seg.residual_rhs << ";\n";
-    os << "    double " << sv << "_J = " << seg.jacobian_rhs << ";\n";
-    os << "    if (std::abs(" << sv << "_R) < " << seg.tol << ") {\n      "
-       << sv << "_converged = true;\n      break;\n    }\n";
+    os << "    double " << p << "_R = " << seg.residual_rhs << ";\n";
+    os << "    double " << p << "_J = " << seg.jacobian_rhs << ";\n";
+    os << "    if (std::abs(" << p << "_R) < " << seg.tol << ") {\n      " << p
+       << "_converged = true;\n      break;\n    }\n";
+    os << "    if (" << p << "_iter == " << seg.max_iter
+       << ") {\n      break;\n    }\n";
     // A singular / ill-scaled Jacobian makes the update non-finite; stop rather
-    // than propagate inf/nan through further iterations. Near-singular-but-finite
-    // steps are caught by the convergence flag (the loop exhausts max_iter).
-    os << "    double " << sv << "_d = " << sv << "_R / " << sv << "_J;\n";
-    os << "    if (!std::isfinite(" << sv << "_d)) {\n      break;\n    }\n";
-    os << "    " << sv << " -= " << sv << "_d;\n";
+    // than propagate inf/nan. (A near-singular-but-finite step is not caught
+    // here — it either fails to converge within max_iter, or, like any
+    // Newton-without-globalization, may reach a spurious root; that is inherent,
+    // not introduced by this guard.)
+    os << "    double " << p << "_d = " << p << "_R / " << p << "_J;\n";
+    os << "    if (!std::isfinite(" << p << "_d)) {\n      break;\n    }\n";
+    os << "    " << sv << " -= " << p << "_d;\n";
     os << "  }\n";
     detail::emit_newton_failure_policy(
-        os, sv + "_converged", {sv}, seg.throw_on_failure, model.name(),
+        os, p + "_converged", {sv}, seg.throw_on_failure, model.name(),
         "state variable '" + sv + "'");
   }
 
@@ -1831,35 +1906,22 @@ inline auto render_compute_function(
     auto suffixes = la.local_suffixes();
     suffixes.emplace_back("iter");      // render's loop counter
     suffixes.emplace_back("converged"); // issue #85 convergence flag
-    std::set<std::string> reserved_ids;
-    for (auto const &[nm, _] : model.scalar_symbol_map())
-      reserved_ids.insert(nm);
-    for (auto const &[nm, _] : model.tensor_symbol_map())
-      reserved_ids.insert(nm);
-    for (auto const &o : model.outputs()) {
-      reserved_ids.insert(o.name);
-      reserved_ids.insert(o.name + "_out");
-    }
+    auto const reserved_ids = detail::newton_reserved_ids(model);
     for (auto const &sys : newton_systems) {
-      std::string p = sys.unknowns[0]; // seed; mangled below until collision-free
-      auto collides = [&](std::string const &pre) {
-        for (auto const &s : suffixes) {
-          if (reserved_ids.contains(pre + "_" + s)) return true;
-        }
-        return false;
-      };
-      while (collides(p)) {
-        p += "_";
-      }
+      auto const p =
+          detail::collision_free_prefix(sys.unknowns[0], reserved_ids, suffixes);
       for (auto const &u : sys.unknowns) {
         os << "  double " << u << " = " << u << "_old;\n";
       }
       os << "  bool " << p << "_converged = false;\n"; // issue #85
-      os << "  for (int " << p << "_iter = 0; " << p << "_iter < "
+      // `<= max_iter`: the residual of the final update is checked before exit
+      // (see the scalar path); emit_newton_step caps the update at max_iter.
+      os << "  for (int " << p << "_iter = 0; " << p << "_iter <= "
          << sys.max_iter << "; ++" << p << "_iter) {\n";
       os << sys.loop_local_decls; // shared CSE temps (4-indented)
       la.emit_newton_step(os, p, sys.unknowns, sys.residual_rhs,
-                          sys.jacobian_rhs, sys.tol, p + "_converged");
+                          sys.jacobian_rhs, sys.tol, p + "_converged",
+                          p + "_iter", sys.max_iter);
       os << "  }\n";
       std::string desc = "coupled system {";
       for (std::size_t i = 0; i < sys.unknowns.size(); ++i)

tests/LinearAlgebraEmitterTest.cpp

@@ -38,7 +38,8 @@ TEST(LinearAlgebraEmitter, EigenEmitsDenseSolveBody) {
   EigenLinearAlgebraEmitter const e;
   std::ostringstream os;
   e.emit_newton_step(os, "sys", {"a", "b"}, {"Ra", "Rb"},
-                     {{"J00", "J01"}, {"J10", "J11"}}, 1e-10, "sys_converged");
+                     {{"J00", "J01"}, {"J10", "J11"}}, 1e-10, "sys_converged",
+                     "sys_iter", 25);
   auto const src = os.str();
   // Fixed-size matrices, row-major fill, ∞-norm convergence, partial-pivot LU.
   EXPECT_NE(src.find("Eigen::Matrix<double, 2, 1> sys_r"), std::string::npos)
@@ -51,8 +52,9 @@ TEST(LinearAlgebraEmitter, EigenEmitsDenseSolveBody) {
       << src;
   EXPECT_NE(src.find("sys_J.partialPivLu().solve(sys_r)"), std::string::npos)
       << src;
-  // issue #85: convergence recorded, and the singular solve is guarded.
+  // issue #85: convergence recorded, update capped at max_iter, singular guard.
   EXPECT_NE(src.find("sys_converged = true"), std::string::npos) << src;
+  EXPECT_NE(src.find("sys_iter == 25"), std::string::npos) << src;
   EXPECT_NE(src.find("sys_dx.allFinite()"), std::string::npos) << src;
   // Each unknown updated from the solution vector.
   EXPECT_NE(src.find("a -= sys_dx(0)"), std::string::npos) << src;
@@ -76,7 +78,8 @@ TEST(LinearAlgebraEmitter, ArmadilloEmitsDenseSolveBody) {
   ArmadilloLinearAlgebraEmitter const a;
   std::ostringstream os;
   a.emit_newton_step(os, "sys", {"a", "b"}, {"Ra", "Rb"},
-                     {{"J00", "J01"}, {"J10", "J11"}}, 1e-10, "sys_converged");
+                     {{"J00", "J01"}, {"J10", "J11"}}, 1e-10, "sys_converged",
+                     "sys_iter", 25);
   auto const src = os.str();
   // Fixed-size types, element assignment, ∞-norm via arma::norm, arma::solve.
   EXPECT_NE(src.find("arma::vec::fixed<2> sys_r"), std::string::npos) << src;
@@ -90,6 +93,7 @@ TEST(LinearAlgebraEmitter, ArmadilloEmitsDenseSolveBody) {
   EXPECT_NE(src.find("arma::solve(sys_dx, sys_J, sys_r)"), std::string::npos)
       << src;
   EXPECT_NE(src.find("sys_converged = true"), std::string::npos) << src;
+  EXPECT_NE(src.find("sys_iter == 25"), std::string::npos) << src;
   EXPECT_NE(src.find("a -= sys_dx(0)"), std::string::npos) << src;
   EXPECT_NE(src.find("b -= sys_dx(1)"), std::string::npos) << src;
   // No Eigen leakage.
@@ -106,7 +110,8 @@ TEST(LinearAlgebraEmitter, EveryBackendEmitsItsUsageMarker) {
   for (auto const *la : backends) {
     std::ostringstream os;
     la->emit_newton_step(os, "p", {"a", "b"}, {"Ra", "Rb"},
-                         {{"J00", "J01"}, {"J10", "J11"}}, 1e-10, "p_converged");
+                         {{"J00", "J01"}, {"J10", "J11"}}, 1e-10, "p_converged",
+                         "p_iter", 10);
     auto const src = os.str();
     EXPECT_NE(src.find(la->usage_marker()), std::string::npos)
         << "marker '" << la->usage_marker() << "' absent from its emission";