Interpolate to tstops by step length and handle sub-dtmin stop spacing

[ChrisRackauckas-Claude]

Ignore this draft until reviewed by @ChrisRackauckas.

What changed and why

#556 (merged as d6be47f) made adaptive steps land on pending tstops. Its review and six rounds of review of this PR found silent partial integrations, each coming from absolute thresholds: the 1e-14 endpoint snap, the 1e-14 dtmin, and the 100eps tstop tolerance, interacting with step control and tstops. This PR replaces those thresholds with one invariant, in all nine adaptive EnsembleGPUKernel steppers:

An adaptive solve is marked complete (t = tf) only by an accepted step whose length is exactly tf - t. A step whose state is not finite is never accepted, and only a landing on a stop before tf may retry with the controller's covering step, so no accepted step reaches past tf. Only the rounding of the time variable in t + (tf - t) is absorbed; no remainder of the dynamics ever is.

Implementation, through shared @inline helpers in integrator_utils.jl:

1. _bounded_step chooses every trial step. It never goes past the next landing target, which is the next pending tstop after t, otherwise tf. It becomes exactly the remaining distance when the proposal reaches the target or would leave less than 1% of a step. Steps therefore end on stops and on tf by construction. The state at a stop is the step's own result, and its dense output stays intact for saveat.
2. The accept branch has no threshold snap: a pending stop inside the step (stop - t <= dt, no absolute tolerance) is landed on; otherwise t = tf if dt == tf - t, else t += dt.
3. The dtmin check errors when the step is below 1e-14 or too small to change t. The only exemption is a step of exactly the distance to its landing target (a stop or tf). So a rejected or otherwise shortened step below dtmin fails with dt<dtmin instead of being snapped over.
4. Every pending stop after t is a landing target. A step that lands exactly on its target (a stop or tf) is exempt from dtmin, so close stops, including adjacent stops one ulp apart, are reached by integrated steps. A step shortened only to land on a stop keeps the controller's earlier proposal for the next step, and dtnew is not clamped after acceptance (the next trial step is bounded anyway), so a one-ulp landing step does not collapse the step size. Every stop before tf is honored and its callbacks run, including terminating ones.

The Float32 endpoint stall from #544 cannot recur: the final step is exactly tf - t, and a step that does not advance t errors instead of looping.

Bounding the step by tf also stops a first step from overshooting tf. The reproducer of #535 now gives [0.99004984, 0.3680994] (exact [0.99004984, 0.36787945], test rtol 1e-2), where master gives [0.98403025, 0.20528531]. That overlaps #537.

5. Converts each pending stop to the integrator's time type before using it. A Float64 tstops array on a Float32 problem otherwise made the bounded step Float64, and Kvaerno3/Kvaerno5 threw MethodError building their nonlinear solver (round-7 review; commit 51c325d).

6. Solves badly scaled small linear systems by equilibration. The 2×2/3×3 linear_solve uses Cramer's rule without scaling. For a Rosenbrock W = I/(γ·dt) − J with a tiny landing step (e.g. dt = 1e-20), the determinant overflowed, and Rodas4/Rodas5P silently returned a zero (or NaN) increment (round-9 review). When the determinant or the solution is not finite, linear_solve now solves the system divided by its largest entry. The solution is unchanged in exact arithmetic, and well-scaled systems keep their exact unscaled arithmetic (results are bitwise identical to before). This also covers Kvaerno3/5's Newton solves, which use the same linear_solve.

7. Never accepts a non-finite step. If a trial step's state or scaled error is not finite, it is not accepted. A step shortened to land on a stop retries with the controller's covering step and lands on the stop by interpolation, as master does. For tiny landing steps, the Rosenbrock C/dt terms and W = J - M/(γ dt) can overflow at any fixed threshold (a scaled mass matrix defeats one), so this guard replaces a fixed landing floor. Any other non-finite step is rejected and shrunk, and ends in dt<dtmin if it stays non-finite. The check is on the norm's inputs, through a non-short-circuit mapreduce(isfinite, &, …): ODE_DEFAULT_NORM uses fast-math, under which isfinite on its result folds to true and NaN > 1 would accept the step. The retry lives in the rejection branch rather than as a continue, because continue in the error-control loop produced SPIR-V that fails OpenCL validation.

8. Saves step endpoints exactly. savevalues! returns the accepted state (or the previous one) for saveat times equal to the step's end (or start). The dense output, which is not exact at Θ = 0 or 1 for every method (Vern9 Float32: 4.7e-4 at Θ = 1), is used only strictly inside a step.

9. Never accepts a step past tf. The covering retry applies only to a landing on a stop before tf, so it never reaches past tf. A final interval whose step is not finite is rejected and shrunk until it fails with an explicit dt<dtmin, rather than being interpolated. A state obtained by interpolating to a covered stop that is not finite fails with "non-finite state at tstop".

Regressions on master, public reproducers

case	master 5e8faee	this branch
Float32 t0 = 2^26, tf = t0 + 64, dt = 32, tstops = [t0 + 24], Tsit5 (expect 64)	72.0	64.0
Float64 stops [0.8125, nextfloat(0.8125)], callback +10 at each (expect 21.25)	dt<dtmin	21.25
stop prevfloat(2.0), callback +10 then terminate! (expect 11.25 at t = prevfloat(2))	dt<dtmin	11.249999999999991 at t = 1.9999999999999998
stops [prevfloat(2)], [2 - 1.5e-14, 2 - 5e-15], [prevfloat(2), 2] with callbacks, all 9 algorithms	dt<dtmin for all 27	all pass
t0 = 2^26 (Float32) / 2^55 (Float64), tf = t0 + 32, tstops = [t0 + 24], saveat = [t0, t0 + 8] (expect 8)	8 for all four stiff methods	8 for all four stiff methods

The fourth and fifth rows are the round-4 review's reproducers (endpoint_stops.jl, coarse_saveat_types.jl). On this PR's previous head 0d0f775 they gave 45 pass / 27 fail (near-tf callbacks silently dropped, terminating callback skipped: final 1.25 at t = 2.0) and 6.0–6.75 instead of 8 for Rodas4, Rodas5P, Kvaerno3 and Kvaerno5.

Tests

test/gpu_kernel_de/adaptive_endpoint.jl (CPU group), all nine adaptive algorithms, 451 assertions:

• Termination: step! with an iteration cap, Float32 and Float64, over 200 tf values × λ ∈ {1, 10, 100} × reltol ∈ {1e-3, 1e-4, 1e-6} for u' = (-u₁, -λu₂). Every configuration must reach t == tf.

• Snap: u' = (1, 1) from t0 = 2.094426f7 to tf = 9.4428696f7 with dt0 = tf - t0, where t0 + (tf - t0) rounds one ulp (8) below tf. Each algorithm must reach tf in one step. For the six methods that integrate u' = 1 exactly, the final state must equal u0 + (tf - t0).

• One tstop before the endpoint (step!): t0 = 0.027657658f0, tf = 1, dt0 = prevfloat(tf - t0), tstops = [0.5]. The stop is visited, dtnew > 0 whenever t < tf, and t == tf at the end.

• One tstop, public solve: tspan = (0.75f0, 1f0), dt = prevfloat(0.25f0), tstops = [0.875]. Final t == 1, final state ≈ 1.25 (rtol 1e-6).

• Multiple tstops (step!): same span, tstops = [0.8125, 0.875, 0.9375]. Each stop is visited exactly once and in order, times are sorted, and t == tf at the end.

• Callback at an intermediate tstop, public solve: the reproducer from the round-10 review, a DiscreteCallback adding 10 at t == 0.875 with merge_callbacks = true. Final t == 1 and final state ≈ 11.25.

• Tstop inside a coarse step, public solve: Float32, u' = (1, 1), u0 = 0, t0 = 2^26, tf = t0 + 32, dt = 32, tstops = [t0 + 24]. One ulp is 8 here, and the initial step of 32 would cross the stop. The final state must be ≈ 32. (tf = t0 + 32 keeps every step a whole number of ulps, so the test does not depend on Adaptive kernel steppers: t += dt rounding decouples time from integrated length at large t #567.)

• Adjacent tstops with callbacks, public solve: Float64, tspan = (0.75, 2), dt = 0.0625, tstops = [0.8125, nextfloat(0.8125)], and a DiscreteCallback adding 10 at each stop. The final state must be ≈ 21.25, i.e. both callbacks fire.

• Tstops at the endpoint (callbacks, public solve): stops [2], [prevfloat(2)], [2 - 1.5e-14, 2 - 5e-15] and [prevfloat(2), 2] on tspan = (0.75, 2), with a callback adding 10 at each stop. The test checks that every callback fires and t == 2 at the end.

• Terminating callback one ulp before tf: the solve must stop at prevfloat(2) with the callback's state.

• saveat inside a step that reaches a tstop: Float32 t0 = 2^26 and Float64 t0 = 2^55 (one ulp = 8), tstops = [t0 + 24], saveat = [t0, t0 + 8]. The sample at t0 + 8 must be 8.

• Rejected sub-dtmin final interval: u' = -1e15 u, tspan = (0, 4e-15), dt = 4e-15. Error control rejects the single full-interval step. The solve must either fail with dt<dtmin or return exp(-4) (rtol 1e-5), never a different value with success.

• Stiff dynamics switched on just before tf: a callback at the stop 1 - 4e-15 sets u₁ = 1e15 in u₂' = -u₁u₂, with the same pass condition against exp(-1e15 (tf - stop)).

• Property grid: for each algorithm, every combination of reltol ∈ {1e-6, 1e-9}, interval length L ∈ {1e-13, 1e-3, 1}, stiffness λL ∈ {0.1, 1, 4} for u' = -λu, initial dt ∈ {0.5, 0.999, 1}·L, and stops ∈ {none, L/4, 0.99L, prevfloat(L), L − 5e-15}. Each solve must either fail with dt<dtmin or end within 2000·reltol (relative) of exp(-λL). On this head every method's error stays below 900·reltol (Rosenbrock23, second order, is the largest); a skipped remainder produces errors of 1e-4 to 0.3. The round-6 reproducer (u' = -1e12u, L = 1e-13, dt = 5e-14, stop 2.5e-14, Rodas5P/Kvaerno5) is one of these cases.

• Float64 tstops on a Float32 problem: tspan = (0.75f0, 1f0), dt = 0.125f0, tstops = [0.875]; final t == 1 and state ≈ 1.25. It errors for Kvaerno3/Kvaerno5 on ceeed40 (310 pass / 2 errors) and passes on 51c325d (312/312).

• Stop closer than dtmin, public solve: Float32 dt = 1f-14, stop prevfloat(dt), u' = (1f14, 0), terminating callback at the stop. The state must be within 64 eps(Float32) of the BigFloat value 0.9999999015. On 51c325d Vern9 returns 1.0004045 (fails); this head passes for all nine.

• Tiny landing step: Float32 u' = (r, r) for r ∈ {1f14, 1f20}, dt = 1f-14, stop 1f-20, terminating callback; the state must be within 64 eps(Float32) of the BigFloat value. On 1b90199 Rodas4/Rodas5P return 0 (1f14) or fail (1f20); this head passes for all nine.

• Stop near floatmin: Float32 stop 5f-38 and Float64 stop 1e-307, u' = (1e14, 1e14), dt = 1e-14, terminating callback. The time must equal the stop, the state must be finite, and it must be within 128 eps of the BigFloat value (Vern7 is checked for finiteness only, because its dense output is GPUVern7 Float32 dense output is ~60 ulps off at mid-step for a constant RHS (visible at tstops) #554).

• Overflowing landing step with a mass matrix: M = 1024 I, stop 1024/floatmax(T), Rosenbrock23/Rodas4/Rodas5P, Float32 and Float64. The state must be within 128 eps of the BigFloat value. On 0d2241f Rodas4/Rodas5P return Success with NaN.

• saveat at step endpoints: all nine methods, Float32 and Float64, u' = (1, 1), dt = 0.5, terminating stop 0.125, saveat = [0, 0.0625, 0.125]. The sample at 0 must be exactly 0 and the sample at the stop within 128 eps of 0.125. On 0d2241f Vern9 saves 0.12505849 (Float32).

• Overflowing final landing step: the round-12 final_retry.jl cases (Rosenbrock23/Rodas4/Rodas5P; a tiny interval with an identity mass matrix, and a 1e30/1e300 scaled mass matrix with and without an endpoint callback). NaN right-hand side at tf: nan_at_tf.jl, all nine methods, both precisions. Each case accepts exactly two outcomes: an accurate, finite state at exactly tf, or a non-Success return or explicit error. A successful solve at another time, or with a non-finite state, fails. On fa8c70f, 30 of these checks fail, because the solve returned Success past tf.

The single-tstop and intermediate-callback tests now check the exact state for all nine algorithms. Earlier versions exempted Vern7 because stopping went through its dense output (#554). Steps now end on stops, so Vern7 and Vern9 are exact there (0 ulps; previously 62 and −3).

Before/after, Julia 1.12.7, include-ing the test file:

source	result
master 5e8faee	218 pass / 13 fail / 36 errors; property grid fails for Tsit5, Vern7 and Vern9 (relative errors up to 0.29) and raises dt<dtmin in 63–81 of 135 cases per method
0d0f775	the four round-4 regressions (near-tf stops, terminating callback, stiff saveat)
6a351d7	the round-5 regression (rejected sub-dtmin final step snapped: 0.759 instead of 0.0183)
25d5a75	297 pass / 6 fail: property grid fails for 6 methods, e.g. Rodas5P 0.90714 and Kvaerno5 0.91383 instead of 0.90484 (round-6 reproducer), Vern9 0.0237 instead of 0.0183
ceeed40	303 / 303 of its own tests; mixed-type stops error for Kvaerno3/5
51c325d	320 pass / 1 fail of the current 321: Vern9 sub-dtmin stop
1b90199	335 pass / 4 fail of the current 339: Rodas4/Rodas5P tiny landing step
1fde487	365 pass / 8 fail of the current 373: Rodas4/Rodas5P below-floor stop (NaN)
0d2241f (e8992bf + merge of master 4f726d5)	409 pass / 6 fail of the current 415: Rodas4/Rodas5P mass-matrix landing (NaN), Vern9 saveat endpoint
fa8c70f	421 pass / 30 fail of the current 451: overflowing final landing (Rodas4/Rodas5P, 12) and NaN at tf (all nine methods, 18)
2ead3c4	451 / 451 pass
a2e999d / 8da09d3 (this head; 8da09d3 changes only a comment)	451 / 451 pass

Verification

# a2e999d (8da09d3 changes only a comment); JULIA_PKG_PRECOMPILE_AUTO=0, clean worktree per version, Enzyme v0.13.209 (unpinned; #568 is on master)
GROUP=CPU    julia +lts  -> Tsit5 accuracy | 172 172; Adaptive endpoint termination | 451 451; Testing DiffEqGPU tests passed
GROUP=CPU    julia +1.12 -> Tsit5 accuracy | 172 172; Adaptive endpoint termination | 451 451; Testing DiffEqGPU tests passed
GROUP=OpenCL julia +lts  (pocl) -> Tsit5 accuracy | 112 112; Testing DiffEqGPU tests passed
GROUP=OpenCL julia +1.12 (pocl) -> Tsit5 accuracy | 112 112; Testing DiffEqGPU tests passed
GROUP=Enzyme julia +lts  -> Enzyme ensemble gradients (Float64, adaptive=true) | 1 1; Lorenz parameter gradients | 2 2; Testing DiffEqGPU tests passed
GROUP=Enzyme julia +1.12 -> Enzyme ensemble gradients (Float64, adaptive=true) | 1 1; Lorenz parameter gradients | 2 2; Testing DiffEqGPU tests passed

Review probes (Julia 1.12; final_retry.jl and nan_at_tf.jl on a2e999d, the others on 2ead3c4, which differs only in the past-tf landing):

probe	result
round-10 extreme.jl (252 cases)	252 / 252 (master: 189 / 252)
round-11 mass_floor.jl	6 / 6
round-12 final_retry.jl	6 end at tf with the expected value (Rosenbrock23); 12 fail explicitly with dt<dtmin (Rodas4/Rodas5P, whose final step overflows)
round-12 nan_at_tf.jl	all 18 fail explicitly with dt<dtmin; none returns Success
round-12 retry_callback.jl	18 / 18
round-12 probes.jl (every round-10/11 probe)	871 ok / 43 not ok, identical to the round-12 review's run on fa8c70f; the 43 are the interior dense-output samples discussed below

The interior-sample failures are the dense-output accuracy of these methods, not the landing logic. On master with no tstops at all, the same interpolants give the same errors at the same Θ: Vern9 Float32 2.1e-5 at Θ = 0.5 vs 2.4e-7 at Θ = 0.125, Rodas5P Float64 3.3e-14 vs 1.8e-14, and Vern7 Float32 6e-5 at either Θ. Landing on a stop moves an interior sample to a different Θ of a shorter step, so the observed error changes. See #554; Vern9 and Rodas5P belong in its scope too.

Runic --check, typos, and git diff --check exit 0 on the changed files.

Comparison with master on the round-12 and round-11 reproducers (Julia 1.12)

final_retry.jl (no stop before tf; the final interval itself overflows for Rodas):

case	master 4f726d5	this head
Rosenbrock23 tiny, Float32 / Float64	Success, u = 1.0 (expected 5e-24 / 1e-293)	Success, exact value
Rosenbrock23 scaled, scaled_callback (both precisions)	Success, NaN	Success, exact value (callback applied)
Rodas4 / Rodas5P tiny (both precisions)	Success, u ≈ 1.0 (expected 5e-24 / 1e-293)	explicit dt<dtmin
Rodas4 / Rodas5P scaled, scaled_callback (both precisions)	Success, NaN	explicit dt<dtmin

Master returns Success with a wrong state in all 18 cases. This head is exact in 6 and fails explicitly in 12, and never returns Success with a wrong state.

rodas_extreme.jl (stop at 5f-38 / 1e-307 before tf, terminating callback): identical on both. Rodas4 gives 5.0000008f-24 / 9.999999999999994e-294, and Rodas5P gives 5.0000035f-24 / 1.0000000000000074e-293; all 4 pass. A landing on a stop before tf keeps the covering retry, so it lands on the stop by interpolation exactly as master does.

Compile time

An earlier revision (2ead3c4) let a covering retry reach past tf and land there by interpolation. That made the dense-output call in the accept branch reachable for every problem, including those without tstops, where it is otherwise dead code. On Julia 1.12, compiling the Verner interpolants for ForwardDiff duals inline in the step loop made the CPU group's ForwardDiff test set take 33 to 36 minutes. This head restricts the covering retry to stops, which restores the compile time; the cost is the explicit failure for an overflowing final interval described above.

First adaptive solve, CPU backend, Julia 1.12, Lorenz with 6-partial ForwardDiff.Dual state and parameters (as in test/gpu_kernel_de/forward_diff.jl):

	master 4f726d5	fa8c70f	2ead3c4	a2e999d / 8da09d3 (this head)
GPUVern7	23.4 s	21.4 s	236 s	23.2 s
GPUVern9	50.5 s	53.3 s	1177 s	53.3 s

On 2ead3c4, the remaining methods took 5 to 16 s each. Making the non-finite error path @noinline did not help (221 s for Vern7). A @noinline interpolation would pass the mutable integrator to a non-inlined call, which heap-allocates it in device code, and that fails on OpenCL.

Full GROUP=CPU run, ForwardDiff test set on Julia 1.12: 33 min 29 s on 2ead3c4; 3 min 06 s on a2e999d (8da09d3 changes only a comment).

Risk assessment

• Risk: medium
• Blast radius: step selection, completion, the dtmin check, tstop landing and non-finite-step rejection in all nine adaptive EnsembleGPUKernel steppers; saveat at step endpoints; the overflow fallback in the internal 2×2/3×3 linear_solve. No public API change. Results change for adaptive solves with tstops (steps land on stops) and at rounding level for the final step of every adaptive solve.
• Evidence: 451/451 endpoint/tstop assertions on LTS and 1.12, with failing-before evidence per review round in the tables above; the CPU, OpenCL and Enzyme groups pass on both versions; extreme matrix 252/252. Remaining review-probe failures are interior dense-output samples that master shows at the same Θ (GPUVern7 Float32 dense output is ~60 ulps off at mid-step for a constant RHS (visible at tstops) #554).
• Independent review: GPT-6 Astra (Codex CLI) has reviewed every revision. Its round-13 review (medium risk) found one P1, the Julia 1.12 compile-time regression on 2ead3c4, fixed by this head (see Compile time); round 14 is pending.
• Merge: needs human review — numerics in all adaptive steppers.

Related work and known CI noise

• The one remaining case from the round-7 review is pre-existing, identical on master and before Fix Float32 endpoint stall in adaptive GPU kernel steppers #544, and is Adaptive kernel steppers: t += dt rounding decouples time from integrated length at large t #567, being fixed separately. At large |t| (Float32 t0 = 2^26, Float64 t0 = 2^55), a rejected final step whose t + dt rounds to tf can still complete without integrating the whole interval.
• Adaptive kernel steppers: t += dt rounding decouples time from integrated length at large t #567 (filed from the round-4 review): at large t, t += dt rounding decouples time from the integrated length, with or without tstops (Tsit5 dt = 12 at t = 2^26 gives 60 instead of 64, before Fix Float32 endpoint stall in adaptive GPU kernel steppers #544 and on master). It is not addressed here. In the review's t0 + 64 survey with a stop at t0 + 24, this branch gives 64 for six algorithms. Vern9, Rodas4 and Rodas5P land at 60.3, 65.6 and 63.8 because of Adaptive kernel steppers: t += dt rounding decouples time from integrated length at large t #567; master gives 72 for most algorithms, and gave 61.7/64.0/66.0 for those three before Respect pending tstops in the adaptive next-step clamp #556.
• Stiff dense output across sub-dtmin stop clusters, a pre-existing issue reproduced before Fix Float32 endpoint stall in adaptive GPU kernel steppers #544, is resolved by landing on every stop with an integrated step: the round-9 review measured its tiny/clustered stiff saveat matrix at 27/27 on 1b90199 (not re-measured on this head, whose only further change is the overflow fallback in linear_solve).
• Adaptive kernel: ContinuousCallback just before a tstop suppresses the tstop's DiscreteCallback #558: a ContinuousCallback just before a tstop suppresses that stop's DiscreteCallback; pre-existing.
• GPUVern7 Float32 dense output is ~60 ulps off at mid-step for a constant RHS (visible at tstops) #554: Vern7 Float32 dense-output accuracy.
• Adaptive GPUTsit5 accepts an initial step past tf and reports an incorrect final state #535 / Fix adaptive last-step overshoot past tf #537: see above. This PR's step bound resolves Adaptive GPUTsit5 accepts an initial step past tf and reports an incorrect final state #535's reproducer as a consequence of the invariant, and conflicts textually with Fix adaptive last-step overshoot past tf #537.
• The CUDA test/enzyme/cuda_records.jl:9 failure and the downgrade-lane failure are pre-existing on master and unrelated; they are tracked in CUDA CI: test/enzyme/cuda_records.jl fails on released CUDACore (Enzyme undef-shadow fill! for non-numeric eltypes) #546 and Raise ModelingToolkitBase floor to 1.72 (fixes downgrade CPU lane) #548.

Not verified

• CUDA/AMDGPU/Metal/oneAPI execution. This host has no GPU.
• GROUP=QA was not rerun for this revision. It passed here for the previous head (LTS 25/25, 1.12 27/27). The round-4 reviewer could not reproduce that: LTS failed Aqua's persistent-tasks load deadline (as on the parent), and 1.12 hit a host EMFILE in JET's FolderMonitor. So QA should be treated as unverified for this PR.

🤖 Generated with Claude Code (model: claude-opus-5-5[1m])

https://claude.ai/code/session_01SJfy3ez7u4CqvAerKrT5np