Reuse nonlinear caches for ImplicitDiscreteSolve reinitialization

[ChrisRackauckas-Claude]

Important

Ignore this draft until reviewed by @ChrisRackauckas.

Summary

• Replace ImplicitDiscreteSolve's hand-unrolled loop over private NonlinearSolve cache fields with the public, allocation-free NonlinearSolveBase.solve_cache! driver.
• Preserve the existing per-Newton residual-contraction history through a reusable step observer.
• Reuse the integrator's nonlinear cache during initial-condition solving and reinit!, including rectangular nonlinear least-squares problems.
• Pass the accepted constant extrapolant explicitly as the nonlinear initial guess so a rejected or externally modified iterate cannot leak into the next step.
• Run ImplicitDiscrete initialization before save_start during integrator reinitialization and preserve InitialFailure instead of resetting it afterward.

Dependencies

This draft depends on NonlinearSolve.jl #1100, specifically NonlinearSolveBase.solve_cache! in NonlinearSolveBase 2.40. The tested dependency head was dca80d83c2d74483d1308c70f4315b69a9cbdf85.

QA also exposed an independent SciMLTesting clean-master harness regression. It was bisected and fixed in SciMLTesting.jl #38; both QA runs below used its exact head 0afa642162a91f8e4df88803690664ded5c52daf. This is not a runtime dependency of the patch.

The related reusable homotopy-cache integration is OrdinaryDiffEq.jl #3985. The final design/benchmark follow-up is also recorded on the merged controller PR #3908.

Root cause

The ImplicitDiscrete initialization added in #2624/#2626 only ran on initial integrator construction. Generic OrdinaryDiffEq reinit! saved the raw u0 before the ImplicitDiscrete consistency solve, and its later retcode reset could erase an initialization failure.

ImplicitDiscreteSolve also rebuilt and solved a fresh NonlinearProblem during every initialization. Its time-step path reused a cache, but drove it by reading and mutating private NonlinearSolve fields and helpers. Reinitializing only the parameters left the previous nonlinear iterate in place instead of applying the documented constant extrapolant.

This patch keeps the native IDS Newton/Kantorovich controller rather than replacing it with an artificial per-step HomotopyProblem. It shares the new reusable NonlinearSolve cache-driving path with the homotopy solvers while retaining the substantially faster native recurrence solve measured below.

Correctness and allocation coverage

The new tests cover:

• fresh solve versus reinit! + reused solve, with save_start=true and false;
• nonlinear-cache identity across reinitialization;
• rectangular nonlinear least-squares cache reuse;
• initialization failure and subsequent recovery;
• stale nonlinear iterates being replaced by the accepted constant extrapolant;
• exact zero steady-state allocations for cached step! and public integrator reinit! on Julia 1.11+.

A separate 20-sample warmed probe on Julia 1.12.6 produced:

• step!: first measurement-site warmup 144 B, then 0 B for all 19 samples;
• reinit! with a preallocated u0: first measurement-site warmup 144 B, then 0 B for all 19 samples;
• solve!: 80 B in every sample because it constructs/returns the solution wrapper.

The official regression asserts the two cache operations that are expected to be allocation-free; it does not disguise the solution-construction allocation.

Performance comparison

Twenty warmed Julia 1.12.6 samples, using AutoFiniteDiff consistently and solving the same 50-step implicit recurrence:

Path	Median time	Bytes	Residual calls	Endpoint
fresh IDS	0.1925435 ms	8,896	565	(1.5531006208104776, 0.4090814367409633)
cached IDS	0.103924 ms	192	561	(1.5531006208104776, 0.4090814367409633)
one-shot Kantorovich	3.3344115 ms	901,104	2,244	(1.5531006208114124, 0.4090814367406754)
cached Kantorovich	0.280863 ms	32	2,040	(1.5531006208114124, 0.4090814367406754)
one-shot HomotopySweep	3.339966 ms	891,312	2,550	(1.5531006237647473, 0.4090814358834119)
cached HomotopySweep	0.294148 ms	32	2,346	(1.5531006237647473, 0.4090814358834119)

Cache reuse makes IDS 1.85x faster and removes 8,704 B (97.84%) while preserving the endpoint exactly and using four fewer residual evaluations. Cached IDS is 2.70x faster than cached Kantorovich and 2.83x faster than cached HomotopySweep on this recurrence. Between the two homotopy methods, cached Kantorovich is 4.52% faster and uses 13.04% fewer residual calls than Sweep. The Kantorovich and Sweep endpoint deltas from IDS are approximately (9.35e-13, -2.88e-13) and (2.95e-9, -8.58e-10), respectively.

These numbers support retaining the IDS-specific solve strategy; they are workload-specific rather than a general ranking of homotopy algorithms.

Local validation

All commands ran from the exact upstream master base fd8861edd500f8bba76a4f2a865e304d10881fe8 on Julia 1.12:

• ODEDIFFEQ_TEST_GROUP=Core ... Pkg.test() for ImplicitDiscreteSolve: 53/53 passed, including cached allocation 2/2, reinitialization 15/15, constant extrapolant 2/2, NLLS 3/3, null state 3/3, and initialization failure 3/3.
• GROUP=Integrators_I ... Pkg.test() at the repository root: all active checks passed, including reinitialization 31/31 and integrator interface 36/36. The two existing marked-broken checks remained marked broken.
• ImplicitDiscreteSolve QA against SciMLTesting update #38: JET 1/1 and Aqua 20/20 passed.
• OrdinaryDiffEqCore QA against SciMLTesting update #38: AllocCheck 1/1, JET 30 passed with one pre-existing marked-broken check, and Aqua 19/19 passed.
• Repository-wide Runic --inplace and --check: exit 0.
• git diff --check: exit 0.

No test was skipped, silenced, deleted, loosened, or converted to a warning.

Process notes

The work followed the request to evaluate whether IDS should become a homotopy adapter, then narrowed the actual issue to reusable init/cache solving. It traced the initialization history, benchmarked fresh and cached IDS against both current homotopy solvers, added the required public cache-driver API upstream instead of depending on private NonlinearSolve internals, implemented the downstream reuse path, reproduced the QA failure on clean master, and validated the independently bisected SciMLTesting fix before publishing this draft.

[ChrisRackauckas-Claude]

Release metadata follow-up: the reusable cache API was ultimately registered as NonlinearSolveBase 2.38.0 after the unreleased version sequence was normalized, while this merged PR still left ImplicitDiscreteSolve at the nonexistent 2.40 compat floor. The focused one-line correction is draft PR #4051: #4051. Local functional tests pass 53/53, QA passes JET 1/1 and Aqua 20/20, and whole-repository Runic passes against the registered 2.38.0 release.