Add SUNNonlinearSolver_Auto module that automatically switches between newton and fixed-point based on stiffness

[balos1]

This implements the algorithm from https://doi.org/10.1007/BF01933714 and enables its use from CVODE/CVODES. Usage from ARKODE is a bit more complicated if we want to be true to what the paper does (use the max convergence rate estimate across all stages).

The results wrt the benefits of using this module are mixed and we may still need to explore it more. I am opening the PR now so the implementation can be reviewed and we can decide on if it should be merged.

[balos1]

Some plots demo'ing that the algorithm does work. This is using the KPR problem:

[u]' = [ a  b ] [ (-1 + u^2 - r(t)) / (2u) ] + [ r'(t) / (2u) ]
[v]    [ b  a ] [ (-2 + v^2 - s(t)) / (2v) ]   [ s'(t) / (2v) ]

with b = 0 and varying a (stiffness).

With a=-1e-2 we see that fixed point is mostly used as one might expect:

With a=-1e2 we see a nice mix of fixed point and newton, where newton is used during the periods of rapid change:

Finally, with a=-1e6 we see newton used entirely as one might expect:

Side note: codex wrote all of the code to make these plots. I told it to use the logs from SUNLogger and the tools in suntools to make them, and it did in one shot. Also, thanks to @Steven-Roberts for the discussions.

[balos1]

Here are some stats with a=-1e2, the case where the most switching takes place

-------------------
AUTO with a=-1e2

Final Solver Statistics:
   Internal solver steps = 1208
   Total RHS evals = 1858
   Total number of nonlinear solver iterations = 4184
         newton = 630, fixed-point = 1225
   Total number of nonlinear solver convergence failures = 41
   Total number of Jacobian evaluations = 50
   Total linear solver setups = 181
   Total RHS evals for setting up the linear system = 0

-------------------
Newton with a=-1e2

Final Solver Statistics:
   Internal solver steps = 484
   Total RHS evals = 749
   Total number of nonlinear solver iterations = 746
   Total number of nonlinear solver convergence failures = 0
   Total number of Jacobian evaluations = 10
   Total linear solver setups = 132
   Total RHS evals for setting up the linear system = 0

-------------------
Fixed point with a=-1e2

Final Solver Statistics:
   Internal solver steps = 1782
   Total RHS evals = 1922
   Total number of nonlinear solver iterations = 1919
   Total number of nonlinear solver convergence failures = 0
   Total number of Jacobian evaluations = 0
   Total linear solver setups = 0
   Total RHS evals for setting up the linear system = 0

We can see that the fastest option is using the Newton nonlinear solver, second fastest using the automatic solver, and third using fixed-point. I think a good takeaway is that the automatic solver may be worth trying when you are unsure if your problem is stiff (or how stiff it is). However, its not clear how likely it is that the automatic solver will be the best option w.r.t. time to solution. As in this case, I suspect it might often be the second best.

[balos1]

Another thing to note. The performance of the switching solver does depend (quite a bit) on the values used for the parameters in the algorithm.

[balos1]

I did go ahead and try this with arkode, just letting each stage decide on switching or not, and the results are not great.

--------------------
AUTO, a=-1e2

Final Solver Statistics:
   Internal solver steps = 2482 (attempted = 3833)
   Total RHS evals (Fi) = 41024
   Total linear solver setups = 187
   Total RHS evals for setting up the linear system = 0
   Total number of Jacobian evaluations = 185
   Total number of nonlinear iterations = 46041
         newton = 2859, fixed-point = 24886
   Total number of nonlinear solver convergence failures = 982
   Total number of error test failures = 0
   Total number of failed steps from solver failure = 1351

--------------------
Newton, a=-1e2

Final Solver Statistics:
   Internal solver steps = 211 (attempted = 278)
   Total RHS evals (Fi) = 4676
   Total linear solver setups = 264
   Total RHS evals for setting up the linear system = 0
   Total number of Jacobian evaluations = 109
   Total number of nonlinear iterations = 3401
   Total number of nonlinear solver convergence failures = 108
   Total number of error test failures = 29
   Total number of failed steps from solver failure = 38

--------------------
Fixed point, a=-1e2

Final Solver Statistics:
   Internal solver steps = 272 (attempted = 398)
   Total RHS evals (Fi) = 3983
   Total linear solver setups = 0
   Total RHS evals for setting up the linear system = 0
   Total number of Jacobian evaluations = 0
   Total number of nonlinear iterations = 1990
         newton = 0, fixed-point = 0
   Total number of nonlinear solver convergence failures = 0
   Total number of error test failures = 126
   Total number of failed steps from solver failure = 0

[balos1]

Something else to think about. This moves a norm from within the integrator into the solver modules. In standalone uses of the SUNNonlinearSolver_Newton and SUNNonlinearSolver_FixedPoint modules, its possible this is an extra cost. Although, I imagine most convergence tests would need it, so it may be only in some edge cases.

[Copilot]

Pull request overview

This PR introduces a new SUNNonlinearSolver_Auto implementation that can switch between Newton and fixed-point nonlinear iterations based on stiffness-related runtime metrics, and wires that capability into CVODE/CVODES (and partially ARKODE).

Changes:

• Added SUNNonlinearSolver_Auto module (C + Fortran interface, SWIG, and sundials4py bindings) and integrated it into build/test targets.
• Extended the nonlinear solver API with an optional getdelnrm op (SUNNonlinSolGetDelNrm) and propagated SUN_NLS_SWITCH to integrators to support switching.
• Added docs and examples demonstrating the new solver and switching behavior.

Reviewed changes

Copilot reviewed 71 out of 71 changed files in this pull request and generated 3 comments.

Show a summary per file

File	Description
src/sunnonlinsol/auto/sunnonlinsol_auto.c	Implements the new auto-switching nonlinear solver and switching logic.
include/sunnonlinsol/sunnonlinsol_auto.h	Public API for constructing/configuring the auto solver.
src/sundials/sundials_nonlinearsolver.c	Adds SUNNonlinSolGetDelNrm API entry point and ops initialization.
include/sundials/sundials_nonlinearsolver.h	Extends ops struct, adds SUNNONLINEARSOLVER_HYBRID, and SUN_NLS_SWITCH.
src/sunnonlinsol/newton/sunnonlinsol_newton.c	Stores update norm and stiffness metric; supports switch return code.
src/sunnonlinsol/fixedpoint/sunnonlinsol_fixedpoint.c	Adds convergence-rate estimate, update norm storage, and switch return code.
src/cvode/*, src/cvodes/*, src/arkode/*	Wires hybrid solver type, uses SUNNonlinSolGetDelNrm, handles SUN_NLS_SWITCH.
doc/shared/sunnonlinsol/*	Documents new API and the Auto solver; adds references.
examples/*	Adds new examples demonstrating Auto solver usage.
bindings/sundials4py/*	Adds Auto solver bindings and exposes new enums/functions.

Comments suppressed due to low confidence (7)

src/sunnonlinsol/newton/sunnonlinsol_newton.c:1

• This overwrites content->delnrm (documented/used as the update norm) with a residual norm, and computes stiffr as resnorm / delnrm without guarding against delnrm == 0. This can (a) break SUNNonlinSolGetDelNrm* semantics and integrator usage, and (b) trigger division-by-zero when the correction/update norm is zero. Consider keeping content->delnrm as the update norm and storing the residual norm in a separate local variable (or dedicated field), and guard the ratio computation when delnrm is 0.
  src/sunnonlinsol/fixedpoint/sunnonlinsol_fixedpoint.c:1
• delnrmp may be 0 (or uninitialized if FP_CONTENT(NLS)->delnrm is never initialized before the first solve/iteration), making delnrm/delnrmp a division-by-zero hazard and potentially producing inf/NaN in crate. Additionally, the newly-added delnrm/crate fields should be explicitly initialized in the solver constructor to avoid reading uninitialized memory. Suggested fix: initialize content->delnrm and content->crate during creation, and guard the ratio when delnrmp <= 0 (or when it is not finite).
  src/sunnonlinsol/auto/sunnonlinsol_auto.c:1
• These setters only forward linear-solver callbacks when Newton is currently active. If the Auto solver starts in fixed-point and later switches to Newton, newton_solver may never receive LSetupFn/LSolveFn, leading to a broken Newton path after switching. For an auto-switching solver, these callbacks should be set on the Newton sub-solver regardless of the currently-active type (or stored and applied when switching). Similar concerns likely apply to other setters that currently forward only to the active sub-solver (e.g., max iters / system function), since switching changes which sub-solver needs the configuration.
  src/sunnonlinsol/auto/sunnonlinsol_auto.c:1
• SUNNonlinSolGetNumIters returns the total iterations across all solves for that solver; adding that total into *_niters_total accumulators will overcount (the total is repeatedly added after each solve). This will make SUNNonlinSolGetNumItersByType_Auto report incorrect values. A concrete fix is to remove the manual accumulation and have SUNNonlinSolGetNumItersByType_Auto query SUNNonlinSolGetNumIters(fp_solver, ...) and SUNNonlinSolGetNumIters(newton_solver, ...) directly, or track the previous totals and only add per-solve deltas.
  src/sunnonlinsol/auto/sunnonlinsol_auto.c:1
• nsolves_since_switch is incremented even when the underlying solve returns SUN_NLS_SWITCH. Since the switching logic resets nsolves_since_switch = 0 inside the convergence test right before returning SUN_NLS_SWITCH, this increment makes the post-switch state 1 (off-by-one), which undermines delay gating. Consider incrementing only when retval != SUN_NLS_SWITCH (and possibly only when the solve completes successfully).
  src/sunnonlinsol/auto/sunnonlinsol_auto.c:1
• The constructor does not check whether SUNNonlinSol_FixedPoint or SUNNonlinSol_Newton returned NULL (allocation failure, etc.). Returning an Auto solver with a NULL sub-solver will lead to null dereferences later (and also leaks NLS/content when creation fails partway). Suggested fix: validate both sub-solver pointers, and on failure free any partially-created resources and return NULL (consistent with other SUNDIALS constructors).
  src/sunnonlinsol/auto/sunnonlinsol_auto.c:1
• The switching behavior (threshold/delay gating + SUN_NLS_SWITCH propagation) is new, user-visible behavior that can materially affect solver outcomes. Given the repo already has unit test infrastructure for ARKODE/CVODE/CVODES, please add unit tests covering the Auto solver's switching decision logic (including delay gating and threshold boundaries) and that integrators correctly recover/retry after SUN_NLS_SWITCH.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

You can also share your feedback on Copilot code review. Take the survey.

[Steven-Roberts]

I took a pass through the docs

I think the auto solver will need to be added to https://sundials.readthedocs.io/en/latest/sundials/Install_link.html#nonlinear-solvers