NonlinearSolveAlg: apply smooth_est error smoothing under W-reuse (fixes silent isnewton gating)#3823
Conversation
The SDIRK/ESDIRK smoothed error estimate (est = W \ tmp, smooth_est = true by default) was gated on isnewton(nlsolver), which is false for NonlinearSolveAlg — NSA silently error-controlled on the raw embedded estimate regardless of the smooth_est setting. The same TRBDF2() call could take a 40x different number of steps and produce orders-of-magnitude different accuracy depending on the inner nonlinear solver choice (SciML#3821). Give NonlinearSolveCache a dz buffer and an estimator-only linsolve on W (built only under W-reuse; the inner NonlinearSolve keeps owning the Newton factorization). Its A is re-set whenever _update_nlsolvealg_W! rewrites W in place, so the estimator refactorizes exactly once per W refresh and otherwise reuses the factorization. Replace the two isnewton gates in the generic SDIRK/ESDIRK error estimation with can_smooth_est(nlsolver), true for Newton caches and for NSA when the estimator linsolve exists. The OOP path (no W-reuse on master) keeps the raw estimate, as before. Verified on VdP mu=1e5 at rtol 1e-8 (naccept / final rel err): TRBDF2+NLNewton smooth=true 1133 / 4.3e-2 TRBDF2+NSA smooth=true 1341 / 3.2e-2 (before: 44949 / 2.2e-6) TRBDF2+NLNewton smooth=false 45261 / 9.6e-7 TRBDF2+NSA smooth=false 45237 / 6.8e-7 KenCarp4+NLNewton smooth=true 498 / 3.2e-2 KenCarp4+NSA smooth=true 570 / 1.5e-2 (before: 6386 / 1.5e-7) HIRES rtol 1e-8: TRBDF2+NSA 715 steps / 9.4 ms vs NLNewton 690 / 5.4 ms (before: 2840 steps / 21.1 ms). ROBER unchanged-sane (906 vs 998 steps). The smoothed estimator's own accuracy behavior on relaxation oscillators is a separate design question tracked in SciML#3821; this change only makes the two nonlinear solvers behave identically under it. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
| # caller must fall back to the raw embedded estimate. | ||
| can_smooth_est(nlsolver::AbstractNLSolver) = can_smooth_est(nlsolver.cache) | ||
| can_smooth_est(cache::AbstractNLSolverCache) = isnewton(cache) | ||
| can_smooth_est(cache::NonlinearSolveCache) = cache.linsolve !== nothing |
There was a problem hiding this comment.
what happens on polyalgs? Since they don't always have a linsolve in all algs
There was a problem hiding this comment.
Polyalgs don't interact with this: the estimator linsolve is built from alg.linsolve — the ODE algorithm's linsolve — never from the inner NonlinearSolve algorithm, so heterogeneous/absent linsolves inside a polyalg are irrelevant to it. With TRBDF2() (no linsolve given) it falls through to LinearSolve's default algorithm selection on W's type. Verified on this branch: NonlinearSolveAlg(RobustMultiNewton()), NonlinearSolveAlg(FastShortcutNonlinearPolyalg()), and NonlinearSolveAlg(TrustRegion()) all build the estimator linsolve and report can_smooth_est == true at init.
Separate pre-existing problem this check surfaced: polyalg inner solvers crash the NSA path outright, on master too — initialize! reads cache.cache.u for the resize length check and NonlinearSolvePolyAlgorithmCache has no u field (FieldError, identical on master and this branch; TrustRegion works fine). Filing as its own issue — the NSA internals need NonlinearSolveBase.get_u-style accessors rather than raw field access, which the residual-termination redesign discussed in #3817 would want anyway.
2 participants
Fixes the
NonlinearSolveAlghalf of #3821.Bug
The SDIRK/ESDIRK smoothed error estimate (
est = W \ tmp,smooth_est = truedefault) was gated onisnewton(nlsolver)— false forNonlinearSolveAlg— so NSA silently error-controlled on the raw embedded estimate no matter whatsmooth_estwas set to. SameTRBDF2()call, 40× different step counts and orders-different accuracy depending on the inner nonlinear solver choice.Fix
NonlinearSolveCachegains adzbuffer and an estimator-only linsolve on W, built only under W-reuse. The inner NonlinearSolve keeps owning the Newton factorization; the estimator linsolve refactorizes exactly once per W refresh (_update_nlsolvealg_W!re-sets itsAafter rewriting W in place) and reuses the factorization otherwise — one triangular solve per step.generic_imex_perform_step.jlbecomecan_smooth_est(nlsolver): true for Newton caches, true for NSA when the estimator linsolve exists. OOP (no W-reuse on master) keeps the raw estimate as before — zero behavior change there.Verified locally (VdP μ=1e5, rtol 1e-8; naccept / final rel err)
The
smooth_estknob now acts identically for both solvers. HIRES rtol 1e-8: TRBDF2+NSA 715 steps / 9.4 ms vs NLNewton 690 / 5.4 ms (before: 2840 steps / 21.1 ms) — this closes the TRBDF2/KenCarp4 step-inflation gap from the #3819/#3820 benchmark series. ROBER sanity: 906 vs 998 steps, errors normal.Note this deliberately makes NSA match NLNewton's estimator behavior including its VDP accuracy trap (median true local error 3.3e3× tol on accepted steps — the audit in #3821). Whether the smoothed estimator itself should change is the remaining, separate design question in that issue.
Tests
New
nsa_smooth_est_tests.jl: assertssmooth_estchanges NSA behavior (pre-fix: identical step counts) and that NSA tracks NLNewton under both settings on VdP. Passes locally along with Newton and Linear-Nonlinear suites (Julia 1.12).🤖 Generated with Claude Code