Apply the integrator tolerance to the Krylov linear solve - #4192
Conversation
|
Confirming the FIRK Identical counts and identical single failure to the branch run. |
|
Root cause of the pre-existing FIRK QA failure above, so it does not get mistaken for something this PR should carry: it is a SciMLTesting regression, not an OrdinaryDiffEq one, and it affects the whole monorepo rather than FIRK.
CI brackets the registration exactly: a run at 2026-08-05T03:11 on v2.6.1 was green, a run at 2026-08-05T12:39 on v2.6.2 was red, and the General registration of v2.6.2 merged at 04:01 between them. It is already fixed upstream and merged — SciMLTesting PR #45, commit Deliberately not adding a local Links:
|
|
Re-ran the Core group on the exact pushed commit (the Verification section's numbers were taken before the |
48fea8d to
42e2090
Compare
`dolinsolve` forwarded `reltol` as a `solve!` keyword, but `SciMLBase.solve!(::LinearCache, alg; kwargs...)` reads the tolerances off the cache and drops the keywords, so the Newton path's Krylov solves ran at LinearSolve's `sqrt(eps)` default however tight or loose the integrator was. Rosenbrock and RadauIIA3 already set both tolerances at `LinearSolve.init`; `build_nlsolver` never receives `reltol`, so `dolinsolve` is the only place the Newton path can set it. Apply it through `LinearSolve.update_tolerances!`, factored into `set_linear_reltol!` so the `dolinsolve` and GaussLegendre call sites cannot drift. Solvers needing a concrete `A` and non-scalar tolerances are skipped. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01F75XsVeZ94QH3PmCuq6QUF
42e2090 to
4f6dedc
Compare
…erance `build_nlsolver` zeroes the inner nonlinear tolerances so the integrator owns convergence, and restores a linear tolerance through `linsolve_kwargs` as a fixed `eps^(4//5)`. That value is independent of `reltol`, so a Krylov descent solve runs to ~3e-13 however loose the integrator is: at `reltol = 1e-3` it is ten orders tighter than the step can use, paid in Jacobian-vector products. Set it per step from the integrator instead, the way SciML#4192 made `dolinsolve` do for `NLNewton`, reusing that PR's `set_linear_reltol!` so a solver that needs a concrete `A` is left alone. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
…erance (#4268) `build_nlsolver` zeroes the inner nonlinear tolerances so the integrator owns convergence, and restores a linear tolerance through `linsolve_kwargs` as a fixed `eps^(4//5)`. That value is independent of `reltol`, so a Krylov descent solve runs to ~3e-13 however loose the integrator is: at `reltol = 1e-3` it is ten orders tighter than the step can use, paid in Jacobian-vector products. Set it per step from the integrator instead, the way #4192 made `dolinsolve` do for `NLNewton`, reusing that PR's `set_linear_reltol!` so a solver that needs a concrete `A` is left alone. Co-authored-by: Chris Rackauckas <accounts@chrisrackauckas.com>
2 participants
Closes #4161.
Please ignore until reviewed by @ChrisRackauckas.
What changed and why
dolinsolveforwarded the integrator's tolerance assolve!(linsolve; reltol = integrator.opts.reltol), butSciMLBase.solve!(::LinearCache, ::KrylovJL; kwargs...)takes its tolerances off the cache and ignores the keywords. The argument was dead, so every Krylov solve on the Newton path ran at LinearSolve'ssqrt(eps) = 1.49e-8default regardless of what the user asked for. Same pattern at the one other site,firk_perform_step.jl(GaussLegendre's stage solve).The fix applies the tolerance through
LinearSolve.update_tolerances!, which is the entry point that takes effect, and stops passing the dead keyword. It is factored intoset_linear_reltol!so the two call sites cannot drift.I made the tolerance live rather than deleting the argument. Two reasons, both checked:
rosenbrock_caches.jl:257) and RadauIIA3 (firk_caches.jl:116) passabstol = reltol, reltol = reltoltoLinearSolve.init, so their Krylov solves have always run at the integrator tolerance. The Newton path is the outlier:build_nlsolvernever receivesreltol, so it cannot set it atinit, anddolinsolvewas the only place left that could. This is not a new regime, it is the shipped behaviour of the most-used stiff family.set_linear_reltol!skips solvers that need a concreteA(a factorization has no tolerance andupdate_tolerances!throws for one) and skips non-scalarreltol, since aLinearCacheholds one scalar tolerance. Rosenbrock and RadauIIA3 caches are already at the right value, so the guard makes them a no-op and their behaviour is bit-identical.Failing before / passing after
All evidence below was re-run on the rebased commit (
4f6dedcf5, on top of masterb9889d3ae), not carried over from the pre-rebase branch.The new test asserts the integrator's tolerance reaches the linear cache. To show it discriminates, the two behavioural call sites were reverted to
solve!(linsolve; reltol)/LinearSolve.solve!(linsolve; reltol = integrator.opts.reltol)— i.e. exactly master's behaviour — withset_linear_reltol!and the committed test file left untouched:The one pass is the
Rodas5Pcase, which confirms Rosenbrock was already correct. With the call sites restored, same command, same environment:The original report also reproduced directly on the pre-rebase base:
FBDF(linsolve = KrylovJL_GMRES())on a 2D Brusselator (N=32, 2048 ODEs) atreltol = abstol = 1e-6left the cache at1.4901161193847656e-8before step 1 and after 30 steps; with the change it reads1.0e-6after the first solve. That reproduction was not re-run after the rebase — the testset above covers the same assertion.Effect on results
Measured against the pre-rebase base (
9d23061e) and not re-run after the rebase over #4193; #4193 changes hownfis counted on this path, so the absolutenfnumbers below no longer match what master reports today. The relative comparison and the error/step columns are unaffected by that. Three problems (2D Brusselator N=32/2048 ODEs, 1D Allen-Cahn 256, 2D nonlinear heat 576), four Newton-based integrators,KrylovJL_GMRES()throughout, error against anFBDFreference atreltol = abstol = 1e-11,BLAS.set_num_threads(1), Julia 1.12.6, LinearSolve 5.6.0.bruss2d wall time (2048 ODEs, the only size where timing is above noise):
At
reltol = 1e-6the linear solve loosens from1.49e-8and costs 11–23% fewer RHS evaluations, with relative error unchanged to three significant figures and step counts unchanged in 22 of 24 cases (Brusselator FBDF 552→553 accepted steps, QNDF 815→821 at 1e-8). Atreltol = 1e-8the integrator tolerance is tighter thansqrt(eps), so the linear solve correctly tightens and costs 1–6% more. No configuration failed and no relative error got worse.Other silently-dropped-kwarg sites
solve!(linsolve; ...)appears at exactly two places in the repo and both are in this diff. No caller ofdolinsolveever passes an explicitreltol— every one takes theintegrator.opts.reltoldefault — so the keyword's only job was the one it was not doing.Worth raising separately, not changed here:
KrylovJL'ssolve!swallowing unrecognised keywords is what let this sit since the LinearSolve v2 migration in 2023. Erroring on them instead would have surfaced it at the first call. That is a LinearSolve change.Also for LinearSolve, and again not changed here:
LinearSolve.update_tolerances!is documented in the LinearSolve manual (docs/src/advanced/algorithm_interface.md) but not declaredpublic, so ExplicitImports flags the qualified access. I added it to this sublib's existing ignore list next toneeds_concrete_A, which is in the same situation; the real fix is apublicdeclaration upstream. Separately,SimpleGMRESsnapshots its tolerances intocachevaland defines noupdate_tolerances_internal!, so it hits LinearSolve'sAbstractKrylovSubspaceMethodfallback and emits one@warn(maxlog = 1) rather than an error; itscachevalis rebuilt wheneverAchanges, so the new tolerance does take effect on the nextWrefresh, and the warning is a false alarm from the conservative default. A one-lineupdate_tolerances_internal!(cache, ::SimpleGMRES, …)upstream would remove it.Rebase over #4193
#4193 (Krylov
nfaccounting) merged first and touched the samedolinsolve. This branch is rebased onto masterb9889d3ae. Four conflicts, all "both sides added at the same place":src/OrdinaryDiffEqDifferentiation.jlpublicblock — kept all three names::drain_jvp_count!, :jvp_counter, :set_linear_reltol!.src/linsolve_utils.jldolinsolvedocstring — merged into one sentence covering both: it now says the tolerance is applied viaset_linear_reltol!and thatstats.nfis charged viadrain_jvp_count!.test/runtests.jl— both includes kept.lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl— kept Count Krylov Jacobian-vector products on the operator that performs them #4193'sdrain_jvp_count!/jvp_counterimports and addedset_linear_reltol!. Count Krylov Jacobian-vector products on the operator that performs them #4193 dropped the now-unusedimport ADTypes; that removal is kept (this branch does not need it, and ExplicitImports would flag it).The
dolinsolvebody merged without a conflict, which is where a silent bug could have hidden, so it was checked rather than trusted. The merged order isset_linear_reltol!(linsolve, reltol)→solve!(linsolve)→drain_jvp_count!(integrator, _alg, linsolve.A). That is correct for both changes and needed no rework: the tolerance must be set before the solve, and the JVP tally must be drained after it. The two do not interact —LinearSolve.update_tolerances!writescache.abstol/cache.reltoland then callsupdate_tolerances_internal!, which forKrylovJLis= nothing(LinearSolvesrc/iterative_wrappers.jl:573); it does not touchcacheval, does not setisfresh, and applies the operator zero times, so it cannot add to or resetJVPCache.njvps. Thenfaccounting test from #4193 was re-run on the rebased commit and passes (below).lib/OrdinaryDiffEqDifferentiationis bumped 3.8.0 → 3.9.0: #4193 took 3.8.0 fordrain_jvp_count!/jvp_counter, and this PR adds a further public name (set_linear_reltol!), so it is another minor bump rather than a reuse of the same version.lib/OrdinaryDiffEqFIRKgoes 2.6.2 → 2.6.3 with itsOrdinaryDiffEqDifferentiationcompat raised"3.8"→"3.9", since it now importsset_linear_reltol!.Verification
All of the following was run on the rebased commit
4f6dedcf5. Julia 1.12.6.GROUP=Core julia --project -e 'using Pkg; Pkg.test()'inlib/OrdinaryDiffEqDifferentiation— note Krylov nf accounting (#4193's test) is included and green, which is the specific regression risk of this merge:GROUP=QAinlib/OrdinaryDiffEqDifferentiation(JET + Aqua; this is where ExplicitImports runs and so where the newLinearSolve.update_tolerances!ignore entry is exercised):GROUP=Coreinlib/OrdinaryDiffEqFIRK(the other call site this touches):GROUP=QAinlib/OrdinaryDiffEqFIRK— this failed on the pre-rebase branch withisempty([:SciMLBase]), and it failed the same way on master. SciMLTesting 2.6.3 is now registered and the group passes:Runic (
Runic.main(["--check","--diff", <changed .jl files>])) exits 0, andtyposover the changed files is clean.activate_group_envrewriteslib/OrdinaryDiffEqDifferentiation/test/qa/Project.tomlduring a QA run; that churn was reverted and is not in this diff.Not verified
docs/and no rendered API entry (set_linear_reltol!is in the sublib's internalpublicblock, likedolinsolveanddrain_jvp_count!beside it, and is not rendered in the OrdinaryDiffEq manual).## Effect on resultsbenchmark tables were not re-run after the rebase; see the note there.dolinsolvecall sites: they go through the same helper and take the same default, but the SDE test groups were not run.abstol/reltolwith a Krylov solver: the new guard skips them, which is the same behaviour as master. Noterosenbrock_caches.jlandfirk_caches.jlpass such an array straight intoLinearSolve.initwithout a guard, which looks like a latent bug on those paths; not touched here.Worth pushing back on
integrator.opts.reltolat all. The adjacent NonlinearSolve inner-algorithm branch inbuild_nlsolverdeliberately uses a fixed tight_inner_lintol = eps^(4/5) ≈ 3e-13for its descent solve instead, so this PR makes the two branches of the same function disagree about the right linear tolerance. I matched Rosenbrock/RadauIIA3 because that is whatdolinsolve's argument said it was doing, but the reverse choice is defensible and a forcing-term (Eisenstat–Walker) scheme would be better than either.abstolis left at LinearSolve's default on the Newton path, so the stopping bound is1.49e-8 + reltol·‖b‖. Rosenbrock and RadauIIA3 set both.dolinsolvehas noabstolargument and adding one is beyond this fix.