Local-Newton robustness: detect non-convergence + guard singular Jacobian (scalar + coupled)

[petlenz]

Summary

The in-function local Newton solvers emit only a max-iteration guard — there is no detection of non-convergence, a singular/ill-conditioned Jacobian, or divergence. On any of those, the generated code silently writes a wrong value to the state-variable out-params with no diagnostic. This affects BOTH the single-scalar path (Phase 3a-2 NewtonSegment) and the coupled N×N path (Phase 3b-2b NewtonSystem, PR #83).

Current behavior (the gap)

The emitted loop is for (iter = 0; iter < max_iter; ++iter) { assemble R, J; if (‖R‖∞ < tol) break; <solve + update>; } followed by an unconditional writeback <sv>_out = <sv>;. So:

1. Non-convergence is silent. If the loop exhausts max_iter without the residual reaching tol, it writes the last (un-converged) iterate — no flag, no error. The caller cannot tell a converged solve from a failed one.
2. Singular Jacobian → silent garbage (coupled path). Eigen::Matrix::partialPivLu().solve(R) does no rank check — on a singular/near-singular J it returns NaN/inf or large garbage, silently (no throw). (Armadillo's arma::solve returns a bool overload that the emitter does not use, so same silent-garbage outcome.)
3. Scalar path divide-by-singular. The single-scalar update sv -= R / J with J ≈ 0 yields inf/nan — at least non-finite (detectable downstream) but still unsignalled.
4. No divergence detection. Residual growth / oscillation isn't caught; the loop runs to max_iter then silently emits whatever it has.

This is consistent with the documented scope: NewtonOptions (recipe.h ~224) carries only tol + max_iter, with the comment "no line search / divergence detection — a plain max-iter guard", and the 3b-2b scope explicitly deferred divergence detection / line search. This issue tracks closing that gap across the whole local-Newton subsystem.

Why it matters

Real plasticity solves DO fail to converge (poor initial guess, load step too large, near-yield-surface corners in multi-surface). A consistent-tangent / stress update that silently returns a bogus value corrupts the host FE solve with no signal — exactly the failure mode hardest to debug downstream. The coupled multi-surface case (the motivation for 3b-2b) is where singular/ill-conditioned local Jacobians are most likely.

Design decision needed: the failure policy

This is the crux — pick a contract (likely a NewtonOptions field):

• Throw / abort — simplest, but generated code in a MOOSE QP loop can't always throw cleanly; MOOSE prefers signalling a cut-back.
• Non-finite sentinel — write NaN to the out-params on failure so the host detects it (matches MOOSE's _console/cutback conventions loosely).
• Status out-param / return code — add a bool& <model>_converged (or per-system flag) the caller checks; most explicit, changes the signature.
• Diagonal regularization / line search — actually improve robustness rather than just report (bigger; the plan flags line search as a later item).

Recommended: start with detection + a signalling policy (status flag or NaN-on-failure, configurable), and treat line search / regularization as a separate enhancement.

Scope

• Emit a post-loop convergence check (did ‖R‖∞ < tol actually hold at exit?).
• Guard the solve: scalar path check |J| (or finiteness of the update); coupled path use a rank-revealing/isfinite check on the Eigen/Armadillo solution (e.g. solve's bool overload for Armadillo; (J*dx - r).norm() residual check or fullPivLu rank for Eigen).
• A NewtonOptions field selecting the failure policy; thread it through NewtonSegment / NewtonSystem (the IR already carries per-system tol/max_iter).
• Tests: a recipe that does not converge in max_iter, and a singular-Jacobian recipe, each asserting the chosen signal fires (not silent).

Affected code

• include/numsim_codegen/recipe.h — render_compute_function scalar-segment loop (sv -= R/J, no guard) + coupled-system loop frame; NewtonOptions (~224).
• include/numsim_codegen/code_emit/linear_algebra_emitter.h — emit_newton_step (the partialPivLu().solve / arma::solve with no failure check).

Related

• PR Phase 3b-2b: coupled vector Newton (#35) #83 (Phase 3b-2b coupled vector Newton) — surfaced this in round-5 review; the gap is pre-existing (inherited from 3a-2) and was explicitly out of Phase 3b-2b: coupled vector Newton (#35) #83's scope, hence this follow-up.
• Phase 3a-2 (Phase 3a-2 — LocalNewtonLoweringPass: in-function scalar Newton iteration #75) single-scalar Newton — same gap.