From 0d0f7759143618f48a67fa01f680ff642d1c9a58 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sun, 27 Sep 2026 12:22:30 -0400 Subject: [PATCH 01/17] Decide tstop state reuse by step length and respect dtmin Keep the step's own state at a tstop only when the accepted step length is exactly the distance to that stop; timestamps within one ulp say nothing about how far past the stop the step integrated (8 time units at t = 2^26 in Float32). Otherwise interpolate to the stop as before. Never clamp the next step to a stop closer than the minimum step size, and leave stops within the snap threshold of `tf` to the endpoint snap. Adjacent representable stops are then reached by the following step's tstop branch instead of forcing `dt < dtmin`. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 5.5 (1M context) Agent-Harness: Claude Code 2.1.280 Agent-Model: claude-opus-5-5[1m] Agent-Session: https://claude.ai/code/session_01SJfy3ez7u4CqvAerKrT5np --- .../integrators/integrator_utils.jl | 14 ++++--- .../perform_step/gpu_kvaerno3_perform_step.jl | 4 +- .../perform_step/gpu_kvaerno5_perform_step.jl | 4 +- .../perform_step/gpu_rodas4_perform_step.jl | 4 +- .../perform_step/gpu_rodas5P_perform_step.jl | 4 +- .../gpu_rosenbrock23_perform_step.jl | 4 +- .../perform_step/gpu_tsit5_perform_step.jl | 4 +- .../perform_step/gpu_vern7_perform_step.jl | 4 +- .../perform_step/gpu_vern9_perform_step.jl | 4 +- test/gpu_kernel_de/adaptive_endpoint.jl | 40 +++++++++++++++++++ 10 files changed, 65 insertions(+), 21 deletions(-) diff --git a/src/ensemblegpukernel/integrators/integrator_utils.jl b/src/ensemblegpukernel/integrators/integrator_utils.jl index 2f59284c..1c9442a4 100644 --- a/src/ensemblegpukernel/integrators/integrator_utils.jl +++ b/src/ensemblegpukernel/integrators/integrator_utils.jl @@ -14,21 +14,25 @@ function build_adaptive_controller_cache(::GPUTsit5IController, ::Type{T}) where return T(1 / 5), zero(T), T(5), T(1 / 5), T(9 / 10), T(1.0e-4), T(1.0e-4) end -# Whether a pending tstop before `tf` falls inside the step just accepted from `integ.t`. +# Whether a pending tstop falls inside the step just accepted from `integ.t`. Stops +# within the endpoint-snap threshold of `tf` are left to the snap, since landing on +# them would leave a remaining interval below the minimum step size. @inline function _tstop_in_step(integ, tf, ::Type{T}) where {T} tstops = integ.tstops (tstops === nothing || integ.tstops_idx > length(tstops)) && return false stop = @inbounds tstops[integ.tstops_idx] - return stop < tf && stop - integ.t - integ.dt - T(100) * eps(T) < T(0) + return tf - stop >= T(1.0e-14) && stop - integ.t - integ.dt - T(100) * eps(T) < T(0) end # The next time the adaptive integrator has to land on exactly: the next pending -# tstop after `integ.t` and before `tf`, otherwise `tf`. -@inline function _next_stop(integ, tf) +# tstop before `tf`, otherwise `tf`. A stop closer than the minimum step size is +# not a clamp target (it would force `dt < dtmin`); the following step's tstop +# branch lands on it instead. +@inline function _next_stop(integ, tf, ::Type{T}) where {T} tstops = integ.tstops if tstops !== nothing && integ.tstops_idx <= length(tstops) stop = @inbounds tstops[integ.tstops_idx] - integ.t < stop < tf && return stop + tf - stop >= T(1.0e-14) && stop - integ.t >= T(1.0e-14) && return stop end return tf end diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl index ba4eca59..fe574366 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl @@ -208,7 +208,7 @@ end if _tstop_in_step(integ, tf, T) integ.t = integ.tstops[integ.tstops_idx] - if abs(integ.t - (t + dt)) > eps(integ.t) + if integ.t - t != dt integ.u = integ(integ.t) end dt = integ.t - integ.tprev @@ -219,7 +219,7 @@ end ##Advance the integrator integ.t += dt end - integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf) - integ.t)) + integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf, T) - integ.t)) end end diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl index 56de80b6..2067e4b7 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl @@ -276,7 +276,7 @@ end if _tstop_in_step(integ, tf, T) integ.t = integ.tstops[integ.tstops_idx] - if abs(integ.t - (t + dt)) > eps(integ.t) + if integ.t - t != dt integ.u = integ(integ.t) end dt = integ.t - integ.tprev @@ -287,7 +287,7 @@ end ##Advance the integrator integ.t += dt end - integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf) - integ.t)) + integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf, T) - integ.t)) end end diff --git a/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl index adecbd81..270c9024 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl @@ -253,7 +253,7 @@ end if _tstop_in_step(integ, tf, T) integ.t = integ.tstops[integ.tstops_idx] - if abs(integ.t - (t + dt)) > eps(integ.t) + if integ.t - t != dt integ.u = integ(integ.t) end dt = integ.t - integ.tprev @@ -264,7 +264,7 @@ end ##Advance the integrator integ.t += dt end - integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf) - integ.t)) + integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf, T) - integ.t)) end end diff --git a/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl index 8638caaf..e50aaa45 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl @@ -327,7 +327,7 @@ end if _tstop_in_step(integ, tf, T) integ.t = integ.tstops[integ.tstops_idx] - if abs(integ.t - (t + dt)) > eps(integ.t) + if integ.t - t != dt integ.u = integ(integ.t) end dt = integ.t - integ.tprev @@ -338,7 +338,7 @@ end ##Advance the integrator integ.t += dt end - integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf) - integ.t)) + integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf, T) - integ.t)) end end diff --git a/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl index 752f0247..70ac28e2 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl @@ -178,7 +178,7 @@ end if _tstop_in_step(integ, tf, T) integ.t = integ.tstops[integ.tstops_idx] - if abs(integ.t - (t + dt)) > eps(integ.t) + if integ.t - t != dt integ.u = integ(integ.t) end dt = integ.t - integ.tprev @@ -189,7 +189,7 @@ end ##Advance the integrator integ.t += dt end - integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf) - integ.t)) + integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf, T) - integ.t)) end end _, saved_in_cb = handle_callbacks!(integ, ts, us) diff --git a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl index e0a67e63..c35ef02a 100644 --- a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl @@ -152,7 +152,7 @@ end if _tstop_in_step(integ, tf, T) integ.t = integ.tstops[integ.tstops_idx] - if abs(integ.t - (t + dt)) > eps(integ.t) + if integ.t - t != dt integ.u = integ(integ.t) end dt = integ.t - integ.tprev @@ -163,7 +163,7 @@ end ##Advance the integrator integ.t += dt end - integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf) - integ.t)) + integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf, T) - integ.t)) end end _, saved_in_cb = handle_callbacks!(integ, ts, us) diff --git a/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl index 28ceeece..3416514c 100644 --- a/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl @@ -188,7 +188,7 @@ end if _tstop_in_step(integ, tf, T) integ.t = integ.tstops[integ.tstops_idx] - if abs(integ.t - (t + dt)) > eps(integ.t) + if integ.t - t != dt integ.u = integ(integ.t) end dt = integ.t - integ.tprev @@ -199,7 +199,7 @@ end ##Advance the integrator integ.t += dt end - integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf) - integ.t)) + integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf, T) - integ.t)) end end _, saved_in_cb = handle_callbacks!(integ, ts, us) diff --git a/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl index f1952034..c2ca2737 100644 --- a/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl @@ -297,7 +297,7 @@ end if _tstop_in_step(integ, tf, T) integ.t = integ.tstops[integ.tstops_idx] - if abs(integ.t - (t + dt)) > eps(integ.t) + if integ.t - t != dt integ.u = integ(integ.t) end dt = integ.t - integ.tprev @@ -308,7 +308,7 @@ end ##Advance the integrator integ.t += dt end - integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf) - integ.t)) + integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf, T) - integ.t)) end end _, saved_in_cb = handle_callbacks!(integ, ts, us) diff --git a/test/gpu_kernel_de/adaptive_endpoint.jl b/test/gpu_kernel_de/adaptive_endpoint.jl index d2139453..e6984973 100644 --- a/test/gpu_kernel_de/adaptive_endpoint.jl +++ b/test/gpu_kernel_de/adaptive_endpoint.jl @@ -152,3 +152,43 @@ end @test all(s -> isapprox(s.u[end], SVector(11.25f0, 11.25f0); rtol = 1.0f-6), sol.u) end end + +# At t ≈ 2^26 a Float32 ulp is 8: the first step (32) ends one ulp past the stop in +# time but has integrated 8 time units beyond it, so the state at the stop has to +# come from the dense output. `tf` is chosen so every later step is a whole number of +# ulps; otherwise `t + dt` rounding would decouple time from the integrated length. +@testset "Adaptive tstop inside a coarse step, public solve ($(nameof(typeof(alg))))" for + alg in ADAPTIVE_ALGS + t0 = Float32(2^26) + tf = t0 + 32.0f0 + prob = ODEProblem{false}( + (u, p, t) -> SVector(1.0f0, 1.0f0), SVector(0.0f0, 0.0f0), (t0, tf) + ) + sol = solve( + EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0); + trajectories = 2, adaptive = true, dt = 32.0f0, tstops = Float32[t0 + 24.0f0], + abstol = 1.0f-6, reltol = 1.0f-3, save_everystep = false + ) + @test all(s -> s.t[end] == tf, sol.u) + @test all(s -> s.u[end] ≈ SVector(32.0f0, 32.0f0), sol.u) +end + +# Stops closer together than the minimum step size must both be visited without +# clamping the next step below `dtmin`. +@testset "Adaptive adjacent tstops with callbacks, public solve ($(nameof(typeof(alg))))" for + alg in ADAPTIVE_ALGS + stops = [0.8125, nextfloat(0.8125)] + cb = DiscreteCallback( + (u, t, integrator) -> t in stops, + integrator -> (integrator.u += SVector(10.0, 10.0)); + save_positions = (false, false) + ) + prob = ODEProblem{false}((u, p, t) -> SVector(1.0, 1.0), SVector(0.0, 0.0), (0.75, 2.0)) + sol = solve( + EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0); + trajectories = 2, adaptive = true, dt = 0.0625, tstops = stops, callback = cb, + merge_callbacks = true, abstol = 1.0e-9, reltol = 1.0e-6, save_everystep = false + ) + @test all(s -> s.t[end] == 2.0, sol.u) + @test all(s -> s.u[end] ≈ SVector(21.25, 21.25), sol.u) +end From 7070a71375418444bde670446ee8e121470ddeb3 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Thu, 1 Oct 2026 12:01:26 -0400 Subject: [PATCH 02/17] Step onto tstops exactly and honor stops next to tf Shorten each trial step so it ends exactly on the next pending tstop (when that stop is at least dtmin away) instead of overshooting and replacing the step's end state with an interpolated one. The step's dense output then stays intact for saveat points it covers, and the state at the stop is the step's own result. Honor every pending stop before `tf`, including stops within dtmin of it, and allow a step shorter than dtmin only as the forced final interval to `tf`. Previously such stops were left to the endpoint snap, which skipped their callbacks (including terminating ones). Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 5.5 (1M context) Agent-Harness: Claude Code 2.1.280 Agent-Model: claude-opus-5-5[1m] Agent-Session: https://claude.ai/code/session_01SJfy3ez7u4CqvAerKrT5np --- .../integrators/integrator_utils.jl | 19 +++-- .../perform_step/gpu_kvaerno3_perform_step.jl | 3 +- .../perform_step/gpu_kvaerno5_perform_step.jl | 3 +- .../perform_step/gpu_rodas4_perform_step.jl | 3 +- .../perform_step/gpu_rodas5P_perform_step.jl | 3 +- .../gpu_rosenbrock23_perform_step.jl | 3 +- .../perform_step/gpu_tsit5_perform_step.jl | 3 +- .../perform_step/gpu_vern7_perform_step.jl | 3 +- .../perform_step/gpu_vern9_perform_step.jl | 3 +- test/gpu_kernel_de/adaptive_endpoint.jl | 81 +++++++++++++++---- 10 files changed, 95 insertions(+), 29 deletions(-) diff --git a/src/ensemblegpukernel/integrators/integrator_utils.jl b/src/ensemblegpukernel/integrators/integrator_utils.jl index 1c9442a4..c36e0f02 100644 --- a/src/ensemblegpukernel/integrators/integrator_utils.jl +++ b/src/ensemblegpukernel/integrators/integrator_utils.jl @@ -14,29 +14,34 @@ function build_adaptive_controller_cache(::GPUTsit5IController, ::Type{T}) where return T(1 / 5), zero(T), T(5), T(1 / 5), T(9 / 10), T(1.0e-4), T(1.0e-4) end -# Whether a pending tstop falls inside the step just accepted from `integ.t`. Stops -# within the endpoint-snap threshold of `tf` are left to the snap, since landing on -# them would leave a remaining interval below the minimum step size. +# Whether a pending tstop before `tf` falls inside the step just accepted from `integ.t`. @inline function _tstop_in_step(integ, tf, ::Type{T}) where {T} tstops = integ.tstops (tstops === nothing || integ.tstops_idx > length(tstops)) && return false stop = @inbounds tstops[integ.tstops_idx] - return tf - stop >= T(1.0e-14) && stop - integ.t - integ.dt - T(100) * eps(T) < T(0) + return stop < tf && stop - integ.t - integ.dt - T(100) * eps(T) < T(0) end # The next time the adaptive integrator has to land on exactly: the next pending # tstop before `tf`, otherwise `tf`. A stop closer than the minimum step size is -# not a clamp target (it would force `dt < dtmin`); the following step's tstop -# branch lands on it instead. +# not a target (stepping to it would force `dt < dtmin`); the following step's +# tstop branch lands on it by interpolation instead. @inline function _next_stop(integ, tf, ::Type{T}) where {T} tstops = integ.tstops if tstops !== nothing && integ.tstops_idx <= length(tstops) stop = @inbounds tstops[integ.tstops_idx] - tf - stop >= T(1.0e-14) && stop - integ.t >= T(1.0e-14) && return stop + stop < tf && stop - integ.t >= T(1.0e-14) && return stop end return tf end +# Shorten a trial step so that it ends exactly on the next pending tstop. The step's +# own end state is then the state at the stop, and its dense output stays intact. +@inline function _clamp_to_next_tstop(integ, dt, tf, ::Type{T}) where {T} + stop = _next_stop(integ, tf, T) + return stop < tf ? min(dt, stop - integ.t) : dt +end + @inline function savevalues!( integrator::SciMLBase.AbstractODEIntegrator{ AlgType, diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl index fe574366..c5ba3f2c 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl @@ -99,6 +99,7 @@ end t = integ.t p = integ.p tf = integ.tf + dt = _clamp_to_next_tstop(integ, dt, tf, T) tmp = integ.tmp f = integ.f @@ -122,7 +123,7 @@ end EEst = convert(T, Inf) while EEst > convert(T, 1.0) - dt < convert(T, 1.0f-14) && error("dt= T(1.0e-14) && error("dt convert(T, 1.0) - dt < convert(T, 1.0f-14) && error("dt= T(1.0e-14) && error("dt convert(T, 1.0) - dt < convert(T, 1.0f-14) && error("dt= T(1.0e-14) && error("dt convert(T, 1.0) - dt < convert(T, 1.0f-14) && error("dt= T(1.0e-14) && error("dt convert(T, 1.0) - dt < convert(T, 1.0f-14) && error("dt= T(1.0e-14) && error("dt T(1.0) - dt < T(1.0e-14) && error("dt= T(1.0e-14) && error("dt convert(T, 1.0) - dt < convert(T, 1.0f-14) && error("dt= T(1.0e-14) && error("dt convert(T, 1.0) - dt < convert(T, 1.0f-14) && error("dt= T(1.0e-14) && error("dt s.t[end] == tf, sol.u) - # Stopping at the tstop goes through the dense output, which for Float32 Vern7 - # is ~60 ulps off even for u' = 1: https://github.com/SciML/DiffEqGPU.jl/issues/554 - if alg isa GPUVern7 - @test all(s -> all(isfinite, s.u[end]), sol.u) - else - @test all(s -> isapprox(s.u[end], SVector(1.25f0, 1.25f0); rtol = 1.0f-6), sol.u) - end + @test all(s -> isapprox(s.u[end], SVector(1.25f0, 1.25f0); rtol = 1.0f-6), sol.u) end # After landing on the first stop, the next proposed step reaches `tf`; the @@ -144,13 +138,8 @@ end abstol = 1.0f-9, reltol = 1.0f-6, save_everystep = false ) @test all(s -> s.t[end] == 1.0f0, sol.u) - # u0 + (tf - t0) + 10 from the callback; Vern7 is only checked for - # the callback having fired until issue 554 is fixed. - if alg isa GPUVern7 - @test all(s -> all(>(11.0f0), s.u[end]), sol.u) - else - @test all(s -> isapprox(s.u[end], SVector(11.25f0, 11.25f0); rtol = 1.0f-6), sol.u) - end + # u0 + (tf - t0) + 10 from the callback + @test all(s -> isapprox(s.u[end], SVector(11.25f0, 11.25f0); rtol = 1.0f-6), sol.u) end # At t ≈ 2^26 a Float32 ulp is 8: the first step (32) ends one ulp past the stop in @@ -192,3 +181,67 @@ end @test all(s -> s.t[end] == 2.0, sol.u) @test all(s -> s.u[end] ≈ SVector(21.25, 21.25), sol.u) end + +# Stops at or within the minimum step size of `tf` must still be visited and their +# callbacks run; only the forced final interval to `tf` may be shorter than `dtmin`. +@testset "Adaptive tstops at the endpoint ($(nameof(typeof(alg))), $label)" for + alg in ADAPTIVE_ALGS, + (label, stops) in ( + ("at tf", [2.0]), ("one ulp before tf", [prevfloat(2.0)]), + ("pair within dtmin of tf", [2.0 - 1.5e-14, 2.0 - 5.0e-15]), + ("one ulp before and at tf", [prevfloat(2.0), 2.0]), + ) + cb = DiscreteCallback( + (u, t, integrator) -> t in stops, + integrator -> (integrator.u += SVector(10.0, 10.0)); + save_positions = (false, false) + ) + prob = ODEProblem{false}((u, p, t) -> SVector(1.0, 1.0), SVector(0.0, 0.0), (0.75, 2.0)) + sol = solve( + EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0); + trajectories = 2, adaptive = true, dt = 0.0625, tstops = stops, callback = cb, + merge_callbacks = true, abstol = 1.0e-9, reltol = 1.0e-6, save_everystep = false + ) + expected = 1.25 + 10 * length(stops) + @test all(s -> s.t[end] == 2.0, sol.u) + @test all(s -> isapprox(s.u[end], SVector(expected, expected); rtol = 1.0e-6), sol.u) +end + +@testset "Adaptive terminating callback one ulp before tf ($(nameof(typeof(alg))))" for + alg in ADAPTIVE_ALGS + stop = prevfloat(2.0) + cb = DiscreteCallback( + (u, t, integrator) -> t == stop, + integrator -> (integrator.u += SVector(10.0, 10.0); terminate!(integrator)); + save_positions = (false, false) + ) + prob = ODEProblem{false}((u, p, t) -> SVector(1.0, 1.0), SVector(0.0, 0.0), (0.75, 2.0)) + sol = solve( + EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0); + trajectories = 2, adaptive = true, dt = 1.25 - 2.0e-14, tstops = [stop], + callback = cb, merge_callbacks = true, abstol = 1.0e-9, reltol = 1.0e-6, + save_everystep = false + ) + expected = stop - 0.75 + 10 + @test all(s -> s.t[end] == stop, sol.u) + @test all(s -> isapprox(s.u[end], SVector(expected, expected); rtol = 1.0e-6), sol.u) +end + +# The step that reaches a stop must keep its dense output intact for the saveat +# points it covers: one ulp at these times is 8, so a step that overshoots the stop +# and is then overwritten with the stop's state would corrupt the sample at t0 + 8. +@testset "Adaptive saveat inside a step that reaches a tstop ($(nameof(typeof(alg))), $T)" for + alg in ADAPTIVE_ALGS, T in (Float32, Float64) + t0 = T(T === Float32 ? 2^26 : 2^55) + @assert eps(t0) == T(8) + prob = ODEProblem{false}( + (u, p, t) -> SVector(one(T), one(T)), SVector(zero(T), zero(T)), (t0, t0 + T(32)) + ) + sol = solve( + EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0); + trajectories = 2, adaptive = true, dt = T(32), tstops = T[t0 + T(24)], + saveat = T[t0, t0 + T(8)], abstol = T(1.0e-6), reltol = T(1.0e-3), + save_everystep = false + ) + @test all(s -> s.u[2] ≈ SVector(T(8), T(8)), sol.u) +end From 6a351d7fb7fefd04286fbd24146af75d37dd8ca4 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Thu, 1 Oct 2026 12:12:23 -0400 Subject: [PATCH 03/17] Update the coarse-step tstop test comment Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 5.5 (1M context) Agent-Harness: Claude Code 2.1.280 Agent-Model: claude-opus-5-5[1m] Agent-Session: https://claude.ai/code/session_01SJfy3ez7u4CqvAerKrT5np --- test/gpu_kernel_de/adaptive_endpoint.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/test/gpu_kernel_de/adaptive_endpoint.jl b/test/gpu_kernel_de/adaptive_endpoint.jl index d5161750..b8110766 100644 --- a/test/gpu_kernel_de/adaptive_endpoint.jl +++ b/test/gpu_kernel_de/adaptive_endpoint.jl @@ -142,10 +142,9 @@ end @test all(s -> isapprox(s.u[end], SVector(11.25f0, 11.25f0); rtol = 1.0f-6), sol.u) end -# At t ≈ 2^26 a Float32 ulp is 8: the first step (32) ends one ulp past the stop in -# time but has integrated 8 time units beyond it, so the state at the stop has to -# come from the dense output. `tf` is chosen so every later step is a whole number of -# ulps; otherwise `t + dt` rounding would decouple time from the integrated length. +# At t ≈ 2^26 a Float32 ulp is 8, so the initial step (32) would end one ulp past the +# stop in time while having integrated 8 time units beyond it. `tf` keeps every step +# a whole number of ulps (see https://github.com/SciML/DiffEqGPU.jl/issues/567). @testset "Adaptive tstop inside a coarse step, public solve ($(nameof(typeof(alg))))" for alg in ADAPTIVE_ALGS t0 = Float32(2^26) From 25d5a75c41a7d283d62234f6b095c182b517339e Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Thu, 1 Oct 2026 15:58:42 -0400 Subject: [PATCH 04/17] Exempt only the whole final interval from dtmin A step shorter than dtmin was allowed whenever less than dtmin remained to `tf`. When error control rejected that final interval and shrank the step, the shorter accepted step was then snapped to `tf`, reporting success with only part of the interval integrated. Exempt only a step covering exactly the remaining interval; a rejected sub-dtmin final step fails with `dt Co-Authored-By: Claude Opus 5.5 (1M context) Agent-Harness: Claude Code 2.1.280 Agent-Model: claude-opus-5-5[1m] Agent-Session: https://claude.ai/code/session_01SJfy3ez7u4CqvAerKrT5np --- .../perform_step/gpu_kvaerno3_perform_step.jl | 2 +- .../perform_step/gpu_kvaerno5_perform_step.jl | 2 +- .../perform_step/gpu_rodas4_perform_step.jl | 2 +- .../perform_step/gpu_rodas5P_perform_step.jl | 2 +- .../gpu_rosenbrock23_perform_step.jl | 2 +- .../perform_step/gpu_tsit5_perform_step.jl | 2 +- .../perform_step/gpu_vern7_perform_step.jl | 2 +- .../perform_step/gpu_vern9_perform_step.jl | 2 +- test/gpu_kernel_de/adaptive_endpoint.jl | 44 +++++++++++++++++++ 9 files changed, 52 insertions(+), 8 deletions(-) diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl index c5ba3f2c..7302e9ad 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl @@ -123,7 +123,7 @@ end EEst = convert(T, Inf) while EEst > convert(T, 1.0) - dt < convert(T, 1.0f-14) && tf - t >= T(1.0e-14) && error("dt convert(T, 1.0) - dt < convert(T, 1.0f-14) && tf - t >= T(1.0e-14) && error("dt convert(T, 1.0) - dt < convert(T, 1.0f-14) && tf - t >= T(1.0e-14) && error("dt convert(T, 1.0) - dt < convert(T, 1.0f-14) && tf - t >= T(1.0e-14) && error("dt convert(T, 1.0) - dt < convert(T, 1.0f-14) && tf - t >= T(1.0e-14) && error("dt T(1.0) - dt < T(1.0e-14) && tf - t >= T(1.0e-14) && error("dt convert(T, 1.0) - dt < convert(T, 1.0f-14) && tf - t >= T(1.0e-14) && error("dt convert(T, 1.0) - dt < convert(T, 1.0f-14) && tf - t >= T(1.0e-14) && error("dt s.u[2] ≈ SVector(T(8), T(8)), sol.u) end + +# A final interval shorter than `dtmin` is only exempt from the `dtmin` check as a +# whole. If error control rejects it, the solve must either fail with `dt isapprox(s.u[end][end], exact; rtol = 1.0e-5, atol = 1.0e-9), sol.u) +end + +@testset "Adaptive rejected sub-dtmin final interval ($(nameof(typeof(alg))))" for + alg in ADAPTIVE_ALGS + prob = ODEProblem{false}((u, p, t) -> -1.0e15 * u, SVector(1.0, 1.0), (0.0, 4.0e-15)) + @test final_interval_outcome(exp(-4.0)) do + solve( + EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0); + trajectories = 2, adaptive = true, dt = 4.0e-15, abstol = 1.0e-9, + reltol = 1.0e-6, save_everystep = false + ) + end +end + +@testset "Adaptive stiff dynamics switched on just before tf ($(nameof(typeof(alg))))" for + alg in ADAPTIVE_ALGS + stop, tf = 1.0 - 4.0e-15, 1.0 + cb = DiscreteCallback( + (u, t, integrator) -> t == stop, + integrator -> (integrator.u = SVector(1.0e15, 1.0)); + save_positions = (false, false) + ) + prob = ODEProblem{false}( + (u, p, t) -> SVector(0.0, -u[1] * u[2]), SVector(0.0, 1.0), (0.0, tf) + ) + @test final_interval_outcome(exp(-1.0e15 * (tf - stop))) do + solve( + EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0); + trajectories = 2, adaptive = true, dt = 0.5, tstops = [stop], callback = cb, + merge_callbacks = true, abstol = 1.0e-9, reltol = 1.0e-6, save_everystep = false + ) + end +end From ceeed40d1e66e72cc09540b5e6f93be329e4e51f Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Thu, 1 Oct 2026 18:32:21 -0400 Subject: [PATCH 05/17] Complete adaptive solves only when a step reaches tf exactly Replace the absolute 1e-14 endpoint snap with exact landing. Each trial step is bounded by the next landing target (the next pending tstop at least dtmin away, otherwise `tf`) and becomes exactly the remaining distance when it would reach the target or leave less than 1% of a step. The accept branch sets `t = tf` only for a step of exactly `tf - t`, so completion always means the whole interval was integrated; only the rounding of `t + (tf - t)` is absorbed. A stop lies inside an accepted step only when `stop - t <= dt` (no absolute tolerance). The dtmin check also rejects steps too small to change `t`, and exempts only a step that is exactly the remaining interval. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 5.5 (1M context) Agent-Harness: Claude Code 2.1.280 Agent-Model: claude-opus-5-5[1m] Agent-Session: https://claude.ai/code/session_01SJfy3ez7u4CqvAerKrT5np --- .../integrators/integrator_utils.jl | 28 +++++++++++------ .../perform_step/gpu_kvaerno3_perform_step.jl | 6 ++-- .../perform_step/gpu_kvaerno5_perform_step.jl | 6 ++-- .../perform_step/gpu_rodas4_perform_step.jl | 6 ++-- .../perform_step/gpu_rodas5P_perform_step.jl | 6 ++-- .../gpu_rosenbrock23_perform_step.jl | 6 ++-- .../perform_step/gpu_tsit5_perform_step.jl | 6 ++-- .../perform_step/gpu_vern7_perform_step.jl | 6 ++-- .../perform_step/gpu_vern9_perform_step.jl | 6 ++-- test/gpu_kernel_de/adaptive_endpoint.jl | 31 +++++++++++++++++++ 10 files changed, 73 insertions(+), 34 deletions(-) diff --git a/src/ensemblegpukernel/integrators/integrator_utils.jl b/src/ensemblegpukernel/integrators/integrator_utils.jl index c36e0f02..bb0344e0 100644 --- a/src/ensemblegpukernel/integrators/integrator_utils.jl +++ b/src/ensemblegpukernel/integrators/integrator_utils.jl @@ -14,18 +14,24 @@ function build_adaptive_controller_cache(::GPUTsit5IController, ::Type{T}) where return T(1 / 5), zero(T), T(5), T(1 / 5), T(9 / 10), T(1.0e-4), T(1.0e-4) end -# Whether a pending tstop before `tf` falls inside the step just accepted from `integ.t`. +# Adaptive steps land on pending tstops and on `tf` by construction: `_bounded_step` +# makes a trial step that would reach (or nearly reach) the next landing target exactly +# the remaining distance, and the accept branch only marks `t = tf` for such a step. +# There is no tolerance-based endpoint snap, so completion always means the interval +# was integrated. + +# Whether a pending tstop before `tf` lies inside the step just accepted from `integ.t`. @inline function _tstop_in_step(integ, tf, ::Type{T}) where {T} tstops = integ.tstops (tstops === nothing || integ.tstops_idx > length(tstops)) && return false stop = @inbounds tstops[integ.tstops_idx] - return stop < tf && stop - integ.t - integ.dt - T(100) * eps(T) < T(0) + return stop < tf && stop - integ.t <= integ.dt end # The next time the adaptive integrator has to land on exactly: the next pending -# tstop before `tf`, otherwise `tf`. A stop closer than the minimum step size is -# not a target (stepping to it would force `dt < dtmin`); the following step's -# tstop branch lands on it by interpolation instead. +# tstop before `tf`, otherwise `tf`. A stop closer than the minimum step size is not +# a target (stepping to it would force `dt < dtmin`); the step that covers it lands on +# it by interpolation instead. @inline function _next_stop(integ, tf, ::Type{T}) where {T} tstops = integ.tstops if tstops !== nothing && integ.tstops_idx <= length(tstops) @@ -35,11 +41,13 @@ end return tf end -# Shorten a trial step so that it ends exactly on the next pending tstop. The step's -# own end state is then the state at the stop, and its dense output stays intact. -@inline function _clamp_to_next_tstop(integ, dt, tf, ::Type{T}) where {T} - stop = _next_stop(integ, tf, T) - return stop < tf ? min(dt, stop - integ.t) : dt +# Trial step from `integ.t`: never past the next landing target, and exactly the +# remaining distance when the proposal reaches it or would leave less than 1% of a +# step before it. +@inline function _bounded_step(integ, dt, tf, ::Type{T}) where {T} + target = _next_stop(integ, tf, T) + remaining = target - integ.t + return remaining <= dt + dt / 100 || integ.t + dt >= target ? remaining : dt end @inline function savevalues!( diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl index 7302e9ad..65583eee 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl @@ -99,7 +99,7 @@ end t = integ.t p = integ.p tf = integ.tf - dt = _clamp_to_next_tstop(integ, dt, tf, T) + dt = _bounded_step(integ, dt, tf, T) tmp = integ.tmp f = integ.f @@ -123,7 +123,7 @@ end EEst = convert(T, Inf) while EEst > convert(T, 1.0) - dt < convert(T, 1.0f-14) && dt != tf - t && error("dt convert(T, 1.0) - dt < convert(T, 1.0f-14) && dt != tf - t && error("dt convert(T, 1.0) - dt < convert(T, 1.0f-14) && dt != tf - t && error("dt convert(T, 1.0) - dt < convert(T, 1.0f-14) && dt != tf - t && error("dt convert(T, 1.0) - dt < convert(T, 1.0f-14) && dt != tf - t && error("dt T(1.0) - dt < T(1.0e-14) && dt != tf - t && error("dt convert(T, 1.0) - dt < convert(T, 1.0f-14) && dt != tf - t && error("dt convert(T, 1.0) - dt < convert(T, 1.0f-14) && dt != tf - t && error("dt Date: Thu, 1 Oct 2026 20:31:58 -0400 Subject: [PATCH 06/17] Keep tstops in the integrator's time type when bounding steps A Float64 `tstops` array for a Float32 problem made the bounded step Float64, and Kvaerno3/Kvaerno5 then failed to build their nonlinear solver. Convert the pending stop to the time type before using it. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 5.5 (1M context) Agent-Harness: Claude Code 2.1.280 Agent-Model: claude-opus-5-5[1m] Agent-Session: https://claude.ai/code/session_01SJfy3ez7u4CqvAerKrT5np --- .../integrators/integrator_utils.jl | 4 ++-- test/gpu_kernel_de/adaptive_endpoint.jl | 15 +++++++++++++++ 2 files changed, 17 insertions(+), 2 deletions(-) diff --git a/src/ensemblegpukernel/integrators/integrator_utils.jl b/src/ensemblegpukernel/integrators/integrator_utils.jl index bb0344e0..3c9ec1c3 100644 --- a/src/ensemblegpukernel/integrators/integrator_utils.jl +++ b/src/ensemblegpukernel/integrators/integrator_utils.jl @@ -24,7 +24,7 @@ end @inline function _tstop_in_step(integ, tf, ::Type{T}) where {T} tstops = integ.tstops (tstops === nothing || integ.tstops_idx > length(tstops)) && return false - stop = @inbounds tstops[integ.tstops_idx] + stop = convert(typeof(tf), @inbounds tstops[integ.tstops_idx]) return stop < tf && stop - integ.t <= integ.dt end @@ -35,7 +35,7 @@ end @inline function _next_stop(integ, tf, ::Type{T}) where {T} tstops = integ.tstops if tstops !== nothing && integ.tstops_idx <= length(tstops) - stop = @inbounds tstops[integ.tstops_idx] + stop = convert(typeof(tf), @inbounds tstops[integ.tstops_idx]) stop < tf && stop - integ.t >= T(1.0e-14) && return stop end return tf diff --git a/test/gpu_kernel_de/adaptive_endpoint.jl b/test/gpu_kernel_de/adaptive_endpoint.jl index 95993ee1..b253c9b1 100644 --- a/test/gpu_kernel_de/adaptive_endpoint.jl +++ b/test/gpu_kernel_de/adaptive_endpoint.jl @@ -319,3 +319,18 @@ end end @test isempty(failures) end + +# `tstops` given as a Float64 array for a Float32 problem: the step bound must stay in +# the integrator's time type. +@testset "Adaptive Float64 tstops on a Float32 problem ($(nameof(typeof(alg))))" for + alg in ADAPTIVE_ALGS + prob = ODEProblem{false}( + (u, p, t) -> SVector(1.0f0, 1.0f0), SVector(1.0f0, 1.0f0), (0.75f0, 1.0f0) + ) + sol = solve( + EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0); + trajectories = 2, adaptive = true, dt = 0.125f0, tstops = [0.875], + abstol = 1.0f-9, reltol = 1.0f-6, save_everystep = false + ) + @test all(s -> s.t[end] == 1.0f0 && s.u[end] ≈ SVector(1.25f0, 1.25f0), sol.u) +end From 1b9019936bba23ab4e435a85570426b22bd3e5f8 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Thu, 1 Oct 2026 23:58:22 -0400 Subject: [PATCH 07/17] Land on every tstop with an integrated step MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A stop closer than dtmin to the current time was not a landing target, so the step covering it landed there by evaluating its dense output near Θ = 1, which is inaccurate for the Float32 Verner interpolants (Vern9 was off by 4e-4 on an exactly integrable problem). Make every pending stop after `t` a landing target and exempt a step that lands exactly on its target from dtmin. A step shortened only to land on a stop keeps the controller's earlier proposal for the next step, and `dtnew` is no longer clamped after acceptance (the next trial step is bounded anyway), so a one-ulp landing step does not collapse the step size. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 5.5 (1M context) Agent-Harness: Claude Code 2.1.280 Agent-Model: claude-opus-5-5[1m] Agent-Session: https://claude.ai/code/session_01SJfy3ez7u4CqvAerKrT5np --- .../integrators/integrator_utils.jl | 13 ++++++----- .../perform_step/gpu_kvaerno3_perform_step.jl | 8 +++++-- .../perform_step/gpu_kvaerno5_perform_step.jl | 8 +++++-- .../perform_step/gpu_rodas4_perform_step.jl | 8 +++++-- .../perform_step/gpu_rodas5P_perform_step.jl | 8 +++++-- .../gpu_rosenbrock23_perform_step.jl | 8 +++++-- .../perform_step/gpu_tsit5_perform_step.jl | 8 +++++-- .../perform_step/gpu_vern7_perform_step.jl | 8 +++++-- .../perform_step/gpu_vern9_perform_step.jl | 8 +++++-- test/gpu_kernel_de/adaptive_endpoint.jl | 22 +++++++++++++++++++ 10 files changed, 77 insertions(+), 22 deletions(-) diff --git a/src/ensemblegpukernel/integrators/integrator_utils.jl b/src/ensemblegpukernel/integrators/integrator_utils.jl index 3c9ec1c3..9219c3b2 100644 --- a/src/ensemblegpukernel/integrators/integrator_utils.jl +++ b/src/ensemblegpukernel/integrators/integrator_utils.jl @@ -17,8 +17,9 @@ end # Adaptive steps land on pending tstops and on `tf` by construction: `_bounded_step` # makes a trial step that would reach (or nearly reach) the next landing target exactly # the remaining distance, and the accept branch only marks `t = tf` for such a step. -# There is no tolerance-based endpoint snap, so completion always means the interval -# was integrated. +# There is no tolerance-based endpoint snap. (At large |t| a step shortened by error +# control can still round onto `tf` in `t += dt`; see +# https://github.com/SciML/DiffEqGPU.jl/issues/567.) # Whether a pending tstop before `tf` lies inside the step just accepted from `integ.t`. @inline function _tstop_in_step(integ, tf, ::Type{T}) where {T} @@ -29,14 +30,14 @@ end end # The next time the adaptive integrator has to land on exactly: the next pending -# tstop before `tf`, otherwise `tf`. A stop closer than the minimum step size is not -# a target (stepping to it would force `dt < dtmin`); the step that covers it lands on -# it by interpolation instead. +# tstop after `integ.t` and before `tf`, otherwise `tf`. Steps end on it with an +# integrated step, never by interpolating back from a longer one; a landing step is +# exempt from `dtmin`. @inline function _next_stop(integ, tf, ::Type{T}) where {T} tstops = integ.tstops if tstops !== nothing && integ.tstops_idx <= length(tstops) stop = convert(typeof(tf), @inbounds tstops[integ.tstops_idx]) - stop < tf && stop - integ.t >= T(1.0e-14) && return stop + integ.t < stop < tf && return stop end return tf end diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl index 65583eee..33126a94 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl @@ -99,7 +99,9 @@ end t = integ.t p = integ.p tf = integ.tf + dtprop = dt dt = _bounded_step(integ, dt, tf, T) + shortened_to_land = dt < dtprop tmp = integ.tmp f = integ.f @@ -123,7 +125,7 @@ end EEst = convert(T, Inf) while EEst > convert(T, 1.0) - (dt < convert(T, 1.0f-14) || t + dt == t) && dt != tf - t && error("dt 1 dt = dt / min(inv(qmin), q11 / gamma) + shortened_to_land = false else # EEst <= 1 q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew + shortened_to_land && (dtnew = max(dtnew, dtprop)) @inbounds begin # Necessary for interpolation integ.k1 = k1 @@ -220,7 +224,7 @@ end ##Advance the integrator integ.t += dt end - integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf, T) - integ.t)) + integ.dtnew = dtnew end end diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl index a08f4775..0c3c6e6e 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl @@ -130,7 +130,9 @@ end t = integ.t p = integ.p tf = integ.tf + dtprop = dt dt = _bounded_step(integ, dt, tf, T) + shortened_to_land = dt < dtprop tmp = integ.tmp f = integ.f @@ -156,7 +158,7 @@ end EEst = convert(T, Inf) while EEst > convert(T, 1.0) - (dt < convert(T, 1.0f-14) || t + dt == t) && dt != tf - t && error("dt 1 dt = dt / min(inv(qmin), q11 / gamma) + shortened_to_land = false else # EEst <= 1 q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew + shortened_to_land && (dtnew = max(dtnew, dtprop)) @inbounds begin # Necessary for interpolation integ.k1 = k1 @@ -288,7 +292,7 @@ end ##Advance the integrator integ.t += dt end - integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf, T) - integ.t)) + integ.dtnew = dtnew end end diff --git a/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl index 3e1f57e9..a69590cc 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl @@ -126,7 +126,9 @@ end t = integ.t p = integ.p tf = integ.tf + dtprop = dt dt = _bounded_step(integ, dt, tf, T) + shortened_to_land = dt < dtprop tmp = integ.tmp f = integ.f @@ -159,7 +161,7 @@ end EEst = convert(T, Inf) while EEst > convert(T, 1.0) - (dt < convert(T, 1.0f-14) || t + dt == t) && dt != tf - t && error("dt 1 dt = dt / min(inv(qmin), q11 / gamma) + shortened_to_land = false else # EEst <= 1 q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew + shortened_to_land && (dtnew = max(dtnew, dtprop)) @inbounds begin # Necessary for interpolation @unpack h21, h22, h23, h24, h25, h31, h32, h33, h34, h35 = integ.tab @@ -265,7 +269,7 @@ end ##Advance the integrator integ.t += dt end - integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf, T) - integ.t)) + integ.dtnew = dtnew end end diff --git a/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl index f27374dc..4675e9e1 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl @@ -163,7 +163,9 @@ end t = integ.t p = integ.p tf = integ.tf + dtprop = dt dt = _bounded_step(integ, dt, tf, T) + shortened_to_land = dt < dtprop tmp = integ.tmp f = integ.f @@ -196,7 +198,7 @@ end EEst = convert(T, Inf) while EEst > convert(T, 1.0) - (dt < convert(T, 1.0f-14) || t + dt == t) && dt != tf - t && error("dt 1 dt = dt / min(inv(qmin), q11 / gamma) + shortened_to_land = false else # EEst <= 1 q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew + shortened_to_land && (dtnew = max(dtnew, dtprop)) @inbounds begin # Necessary for interpolation @unpack h21, h22, h23, h24, h25, h26, h27, h28, h31, h32, h33, h34, @@ -339,7 +343,7 @@ end ##Advance the integrator integ.t += dt end - integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf, T) - integ.t)) + integ.dtnew = dtnew end end diff --git a/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl index 93624b26..3fcfa0e9 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl @@ -81,7 +81,9 @@ end t = integ.t p = integ.p tf = integ.tf + dtprop = dt dt = _bounded_step(integ, dt, tf, T) + shortened_to_land = dt < dtprop tmp = integ.tmp f = integ.f @@ -105,7 +107,7 @@ end mass_matrix = integ.f.mass_matrix while EEst > convert(T, 1.0) - (dt < convert(T, 1.0f-14) || t + dt == t) && dt != tf - t && error("dt 1 dt = dt / min(inv(qmin), q11 / gamma) + shortened_to_land = false else # EEst <= 1 q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew + shortened_to_land && (dtnew = max(dtnew, dtprop)) @inbounds begin # Necessary for interpolation integ.k1 = k1 @@ -190,7 +194,7 @@ end ##Advance the integrator integ.t += dt end - integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf, T) - integ.t)) + integ.dtnew = dtnew end end _, saved_in_cb = handle_callbacks!(integ, ts, us) diff --git a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl index 823aae5f..c637c0c7 100644 --- a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl @@ -76,7 +76,9 @@ end t = integ.t p = integ.p tf = integ.tf + dtprop = dt dt = _bounded_step(integ, dt, tf, T) + shortened_to_land = dt < dtprop a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as btilde1, btilde2, btilde3, btilde4, btilde5, btilde6, btilde7 = integ.btildes @@ -100,7 +102,7 @@ end EEst = convert(T, Inf) while EEst > T(1.0) - (dt < T(1.0e-14) || t + dt == t) && dt != tf - t && error("dt 1 dt = dt / min(inv(qmin), q11 / gamma) + shortened_to_land = false else # EEst <= 1 q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew + shortened_to_land && (dtnew = max(dtnew, dtprop)) @inbounds begin # Necessary for interpolation integ.k1 = k1 @@ -164,7 +168,7 @@ end ##Advance the integrator integ.t += dt end - integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf, T) - integ.t)) + integ.dtnew = dtnew end end _, saved_in_cb = handle_callbacks!(integ, ts, us) diff --git a/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl index 8f7d3ff8..f4c41923 100644 --- a/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl @@ -86,7 +86,9 @@ end t = integ.t p = integ.p tf = integ.tf + dtprop = dt dt = _bounded_step(integ, dt, tf, T) + shortened_to_land = dt < dtprop @unpack c2, c3, c4, c5, c6, c7, c8, a021, a031, a032, a041, a043, a051, a053, a054, a061, a063, a064, a065, a071, a073, a074, a075, a076, a081, a083, a084, @@ -113,7 +115,7 @@ end EEst = convert(T, Inf) while EEst > convert(T, 1.0) - (dt < convert(T, 1.0f-14) || t + dt == t) && dt != tf - t && error("dt 1 dt = dt / min(inv(qmin), q11 / gamma) + shortened_to_land = false else # EEst <= 1 q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew + shortened_to_land && (dtnew = max(dtnew, dtprop)) @inbounds begin # Necessary for interpolation integ.k1 = k1 @@ -200,7 +204,7 @@ end ##Advance the integrator integ.t += dt end - integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf, T) - integ.t)) + integ.dtnew = dtnew end end _, saved_in_cb = handle_callbacks!(integ, ts, us) diff --git a/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl index 5a2188b7..07fb8669 100644 --- a/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl @@ -136,7 +136,9 @@ end t = integ.t p = integ.p tf = integ.tf + dtprop = dt dt = _bounded_step(integ, dt, tf, T) + shortened_to_land = dt < dtprop @unpack c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, a0201, a0301, a0302, a0401, a0403, a0501, a0503, a0504, a0601, a0604, a0605, a0701, a0704, a0705, a0706, @@ -167,7 +169,7 @@ end EEst = convert(T, Inf) while EEst > convert(T, 1.0) - (dt < convert(T, 1.0f-14) || t + dt == t) && dt != tf - t && error("dt 1 dt = dt / min(inv(qmin), q11 / gamma) + shortened_to_land = false else # EEst <= 1 q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew + shortened_to_land && (dtnew = max(dtnew, dtprop)) @inbounds begin # Necessary for interpolation integ.k1 = k1 @@ -309,7 +313,7 @@ end ##Advance the integrator integ.t += dt end - integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf, T) - integ.t)) + integ.dtnew = dtnew end end _, saved_in_cb = handle_callbacks!(integ, ts, us) diff --git a/test/gpu_kernel_de/adaptive_endpoint.jl b/test/gpu_kernel_de/adaptive_endpoint.jl index b253c9b1..1f56f86c 100644 --- a/test/gpu_kernel_de/adaptive_endpoint.jl +++ b/test/gpu_kernel_de/adaptive_endpoint.jl @@ -334,3 +334,25 @@ end ) @test all(s -> s.t[end] == 1.0f0 && s.u[end] ≈ SVector(1.25f0, 1.25f0), sol.u) end + +# A stop less than `dtmin` after `t` must be reached by an integrated step, not by +# evaluating the dense output of a longer step near Θ = 1 (inaccurate for Float32 +# Verner interpolants). `u' = (1e14, 0)` is integrated exactly by every method. +@testset "Adaptive tstop closer than dtmin, public solve ($(nameof(typeof(alg))))" for + alg in ADAPTIVE_ALGS + dt = 1.0f-14 + stop = prevfloat(dt) + cb = DiscreteCallback( + (u, t, integrator) -> t == stop, terminate!; save_positions = (false, false) + ) + prob = ODEProblem{false}( + (u, p, t) -> SVector(1.0f14, 0.0f0), SVector(0.0f0, 1.0f0), (0.0f0, 2dt) + ) + sol = solve( + EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0); + trajectories = 2, adaptive = true, dt, tstops = [stop], callback = cb, + merge_callbacks = true, abstol = 1.0f-9, reltol = 1.0f-6, save_everystep = false + ) + exact = Float64(BigFloat(1.0f14) * BigFloat(stop)) + @test all(s -> s.t[end] == stop && abs(s.u[end][1] - exact) <= 64eps(Float32), sol.u) +end From 1fde4876ed5684be47ee3c55cc0e9c11c2968b1e Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Fri, 2 Oct 2026 01:44:51 -0400 Subject: [PATCH 08/17] Solve badly scaled small W-systems by equilibration MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Cramer's rule in the 2x2/3x3 `linear_solve` forms products of matrix entries without scaling. A Rosenbrock W = I/(γ dt) - J with a tiny landing step has entries near 1/dt, so the determinant overflowed and Rodas4/Rodas5P silently returned a zero (or NaN) increment. When the determinant or the solution is not finite, solve the system divided by its largest entry instead; the solution is unchanged and well-scaled systems keep their exact arithmetic. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 5.5 (1M context) Agent-Harness: Claude Code 2.1.280 Agent-Model: claude-opus-5-5[1m] Agent-Session: https://claude.ai/code/session_01SJfy3ez7u4CqvAerKrT5np --- src/ensemblegpukernel/linalg/linsolve.jl | 26 ++++++++++++++++++++---- test/gpu_kernel_de/adaptive_endpoint.jl | 25 +++++++++++++++++++++++ 2 files changed, 47 insertions(+), 4 deletions(-) diff --git a/src/ensemblegpukernel/linalg/linsolve.jl b/src/ensemblegpukernel/linalg/linsolve.jl index 89b0bd64..e0abd041 100644 --- a/src/ensemblegpukernel/linalg/linsolve.jl +++ b/src/ensemblegpukernel/linalg/linsolve.jl @@ -14,9 +14,8 @@ end @inbounds return similar_type(b, typeof(a[1] \ b[1]))(a[1] \ b[1]) end -@inline function _linear_solve( +@inline function _cramer_solve( ::Size{(2, 2)}, - ::Size{(2,)}, a::StaticMatrix{<:Any, <:Any, Ta}, b::StaticVector{<:Any, Tb} ) where {Ta, Tb} @@ -28,9 +27,8 @@ end ) end -@inline function _linear_solve( +@inline function _cramer_solve( ::Size{(3, 3)}, - ::Size{(3,)}, a::StaticMatrix{<:Any, <:Any, Ta}, b::StaticVector{<:Any, Tb} ) where {Ta, Tb} @@ -56,6 +54,26 @@ end ) end +# Cramer's rule forms products of up to n matrix entries, which overflow (or underflow) +# for badly scaled systems such as a Rosenbrock W = I/(γ dt) - J with a tiny dt. When +# that happens, solve the system divided by its largest entry instead: the solution is +# unchanged, and well-scaled systems keep their exact unscaled arithmetic. +for Sa in [(2, 2), (3, 3)] + @eval @inline function _linear_solve( + ::Size{$Sa}, + ::Size{($Sa[1],)}, + a::StaticMatrix{<:Any, <:Any, Ta}, + b::StaticVector{<:Any, Tb} + ) where {Ta, Tb} + d = det(a) + x = _cramer_solve(Size($Sa), a, b) + isfinite(d) && !iszero(d) && all(isfinite, x) && return x + s = maximum(abs, a) + (isfinite(s) && !iszero(s)) || return x + return _cramer_solve(Size($Sa), a / s, b / s) + end +end + for Sa in [(2, 2), (3, 3)] # not needed for Sa = (1, 1); @eval begin @inline function _linear_solve( diff --git a/test/gpu_kernel_de/adaptive_endpoint.jl b/test/gpu_kernel_de/adaptive_endpoint.jl index 1f56f86c..86369c57 100644 --- a/test/gpu_kernel_de/adaptive_endpoint.jl +++ b/test/gpu_kernel_de/adaptive_endpoint.jl @@ -356,3 +356,28 @@ end exact = Float64(BigFloat(1.0f14) * BigFloat(stop)) @test all(s -> s.t[end] == stop && abs(s.u[end][1] - exact) <= 64eps(Float32), sol.u) end + +# A landing step far below `dtmin` must stay accurate for every method family. For the +# Rosenbrock methods W = I/(γ dt) - J then has entries near 1/dt, which overflowed the +# unscaled Cramer solve and silently zeroed (or NaN'd) the increment. +@testset "Adaptive tiny landing step ($(nameof(typeof(alg))), rate $rate)" for + alg in ADAPTIVE_ALGS, rate in (1.0f14, 1.0f20) + stop = 1.0f-20 + cb = DiscreteCallback( + (u, t, integrator) -> t == stop, terminate!; save_positions = (false, false) + ) + prob = ODEProblem{false}( + (u, p, t) -> SVector(rate, rate), SVector(0.0f0, 0.0f0), (0.0f0, 2.0f-14) + ) + sol = solve( + EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0); + trajectories = 2, adaptive = true, dt = 1.0f-14, tstops = [stop], callback = cb, + merge_callbacks = true, save_everystep = false, abstol = 1.0f-12, reltol = 1.0f-3 + ) + exact = Float64(BigFloat(rate) * BigFloat(stop)) + @test all( + s -> s.t[end] == stop && + all(x -> isapprox(Float64(x), exact; rtol = 64eps(Float32)), s.u[end]), + sol.u + ) +end From e8992bf3b8df441479b25374a1fc5de508a773c3 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Fri, 2 Oct 2026 16:25:47 -0400 Subject: [PATCH 09/17] Interpolate to stops below a landing floor MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Landing on a stop with an integrated step of length near floatmin overflows the `C/dt` terms in the Rosenbrock stage formulas, and the NaN error estimate was accepted. Below `_landing_floor = 1024 / floatmax(T)` (|C| ≤ 166 in Rodas5P), a stop is reached by interpolating the covering step as before, and a final interval shorter than the floor fails `dtmin` instead of running the step. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 5.5 (1M context) Agent-Harness: Claude Code 2.1.280 Agent-Model: claude-opus-5-5[1m] Agent-Session: https://claude.ai/code/session_01SJfy3ez7u4CqvAerKrT5np --- .../integrators/integrator_utils.jl | 11 ++++--- .../perform_step/gpu_kvaerno3_perform_step.jl | 4 ++- .../perform_step/gpu_kvaerno5_perform_step.jl | 4 ++- .../perform_step/gpu_rodas4_perform_step.jl | 4 ++- .../perform_step/gpu_rodas5P_perform_step.jl | 4 ++- .../gpu_rosenbrock23_perform_step.jl | 4 ++- .../perform_step/gpu_tsit5_perform_step.jl | 4 ++- .../perform_step/gpu_vern7_perform_step.jl | 4 ++- .../perform_step/gpu_vern9_perform_step.jl | 4 ++- test/gpu_kernel_de/adaptive_endpoint.jl | 32 +++++++++++++++++-- 10 files changed, 60 insertions(+), 15 deletions(-) diff --git a/src/ensemblegpukernel/integrators/integrator_utils.jl b/src/ensemblegpukernel/integrators/integrator_utils.jl index 9219c3b2..ca9c29cd 100644 --- a/src/ensemblegpukernel/integrators/integrator_utils.jl +++ b/src/ensemblegpukernel/integrators/integrator_utils.jl @@ -29,15 +29,18 @@ end return stop < tf && stop - integ.t <= integ.dt end +# Shortest step taken to land on a stop or `tf`. Below it, the `C/dt` terms of the +# Rosenbrock stage formulas (|C| ≤ 166 for Rodas5P) can overflow, so closer stops are +# reached by interpolating the covering step and a shorter final interval fails `dtmin`. +@inline _landing_floor(::Type{T}) where {T} = T(1024) / floatmax(T) + # The next time the adaptive integrator has to land on exactly: the next pending -# tstop after `integ.t` and before `tf`, otherwise `tf`. Steps end on it with an -# integrated step, never by interpolating back from a longer one; a landing step is -# exempt from `dtmin`. +# tstop at least `_landing_floor` after `integ.t` and before `tf`, otherwise `tf`. @inline function _next_stop(integ, tf, ::Type{T}) where {T} tstops = integ.tstops if tstops !== nothing && integ.tstops_idx <= length(tstops) stop = convert(typeof(tf), @inbounds tstops[integ.tstops_idx]) - integ.t < stop < tf && return stop + stop < tf && stop - integ.t >= _landing_floor(T) && return stop end return tf end diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl index 33126a94..cf98e361 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl @@ -125,7 +125,9 @@ end EEst = convert(T, Inf) while EEst > convert(T, 1.0) - (dt < convert(T, 1.0f-14) || t + dt == t) && dt != _next_stop(integ, tf, T) - t && error("dt convert(T, 1.0) - (dt < convert(T, 1.0f-14) || t + dt == t) && dt != _next_stop(integ, tf, T) - t && error("dt convert(T, 1.0) - (dt < convert(T, 1.0f-14) || t + dt == t) && dt != _next_stop(integ, tf, T) - t && error("dt convert(T, 1.0) - (dt < convert(T, 1.0f-14) || t + dt == t) && dt != _next_stop(integ, tf, T) - t && error("dt convert(T, 1.0) - (dt < convert(T, 1.0f-14) || t + dt == t) && dt != _next_stop(integ, tf, T) - t && error("dt T(1.0) - (dt < T(1.0e-14) || t + dt == t) && dt != _next_stop(integ, tf, T) - t && error("dt convert(T, 1.0) - (dt < convert(T, 1.0f-14) || t + dt == t) && dt != _next_stop(integ, tf, T) - t && error("dt convert(T, 1.0) - (dt < convert(T, 1.0f-14) || t + dt == t) && dt != _next_stop(integ, tf, T) - t && error("dt s.t[end] == stop && abs(s.u[end][1] - exact) <= 64eps(Float32), sol.u) end -# A landing step far below `dtmin` must stay accurate for every method family. For the -# Rosenbrock methods W = I/(γ dt) - J then has entries near 1/dt, which overflowed the -# unscaled Cramer solve and silently zeroed (or NaN'd) the increment. +# Landing on a stop far below `dtmin` keeps the state accurate for every method family, +# including Rosenbrock methods whose W = I/(γ dt) - J then has entries near 1/dt. @testset "Adaptive tiny landing step ($(nameof(typeof(alg))), rate $rate)" for alg in ADAPTIVE_ALGS, rate in (1.0f14, 1.0f20) stop = 1.0f-20 @@ -381,3 +380,30 @@ end sol.u ) end + +# Stops below `_landing_floor` (where Rosenbrock `C/dt` stage terms overflow) are reached +# by interpolating the covering step: the result is finite and accurate, never NaN. Vern7's +# dense output is checked only for finiteness: https://github.com/SciML/DiffEqGPU.jl/issues/554 +@testset "Adaptive stop below the landing floor ($(nameof(typeof(alg))), $T)" for + alg in ADAPTIVE_ALGS, T in (Float32, Float64) + stop = T === Float32 ? 5.0f-38 : 1.0e-307 + rate = T(1.0e14) + cb = DiscreteCallback( + (u, t, integrator) -> t == stop, terminate!; save_positions = (false, false) + ) + prob = ODEProblem{false}( + (u, p, t) -> SVector(rate, rate), SVector(zero(T), zero(T)), (zero(T), T(2.0e-14)) + ) + sol = solve( + EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0); + trajectories = 2, adaptive = true, dt = T(1.0e-14), tstops = [stop], callback = cb, + merge_callbacks = true, save_everystep = false, abstol = T(1.0e-12), reltol = T(1.0e-3) + ) + exact = Float64(BigFloat(rate) * BigFloat(stop)) + @test all(s -> s.t[end] == stop && all(isfinite, s.u[end]), sol.u) + if !(alg isa GPUVern7) + @test all( + s -> all(x -> isapprox(Float64(x), exact; rtol = 128eps(T)), s.u[end]), sol.u + ) + end +end From 7ef20a59c85c0e91c15cdf68bac07fb60055aa75 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Fri, 2 Oct 2026 20:19:44 -0400 Subject: [PATCH 10/17] Never accept a non-finite step; save step endpoints exactly MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A step whose state or scaled error is not finite is never accepted. A step shortened to land on a stop retries with the controller's covering step and lands on the stop by interpolation; any other step is rejected and shrunk. The check is on the norm's inputs because `ODE_DEFAULT_NORM` uses fast-math, under which `isfinite` on its result is folded to true. This replaces the fixed landing floor, which a scaled mass matrix could still overflow. `savevalues!` returns the accepted state (or the previous one) for saveat times at the ends of a step; the dense output, which is not exact at Θ = 0 or 1 for every method, is used only strictly inside the step. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 5.5 (1M context) Agent-Harness: Claude Code 2.1.280 Agent-Model: claude-opus-5-5[1m] Agent-Session: https://claude.ai/code/session_01SJfy3ez7u4CqvAerKrT5np --- .../integrators/integrator_utils.jl | 26 ++++---- .../perform_step/gpu_kvaerno3_perform_step.jl | 11 +++- .../perform_step/gpu_kvaerno5_perform_step.jl | 11 +++- .../perform_step/gpu_rodas4_perform_step.jl | 11 +++- .../perform_step/gpu_rodas5P_perform_step.jl | 11 +++- .../gpu_rosenbrock23_perform_step.jl | 11 +++- .../perform_step/gpu_tsit5_perform_step.jl | 11 +++- .../perform_step/gpu_vern7_perform_step.jl | 11 +++- .../perform_step/gpu_vern9_perform_step.jl | 11 +++- test/gpu_kernel_de/adaptive_endpoint.jl | 59 +++++++++++++++++-- 10 files changed, 148 insertions(+), 25 deletions(-) diff --git a/src/ensemblegpukernel/integrators/integrator_utils.jl b/src/ensemblegpukernel/integrators/integrator_utils.jl index 4d41ebe5..61f0bf17 100644 --- a/src/ensemblegpukernel/integrators/integrator_utils.jl +++ b/src/ensemblegpukernel/integrators/integrator_utils.jl @@ -34,18 +34,15 @@ end return stop < tf && stop - integ.t <= integ.dt end -# Shortest step taken to land on a stop or `tf`. Below it, the `C/dt` terms of the -# Rosenbrock stage formulas (|C| ≤ 166 for Rodas5P) can overflow, so closer stops are -# reached by interpolating the covering step and a shorter final interval fails `dtmin`. -@inline _landing_floor(::Type{T}) where {T} = T(1024) / floatmax(T) - # The next time the adaptive integrator has to land on exactly: the next pending -# tstop at least `_landing_floor` after `integ.t` and before `tf`, otherwise `tf`. +# tstop after `integ.t` and before `tf`, otherwise `tf`. A landing step whose arithmetic +# overflows (non-finite error estimate) is retried with the controller's unshortened +# step, which covers the stop and lands on it by interpolation. @inline function _next_stop(integ, tf, ::Type{T}) where {T} tstops = integ.tstops if tstops !== nothing && integ.tstops_idx <= length(tstops) stop = convert(typeof(tf), @inbounds tstops[integ.tstops_idx]) - stop < tf && stop - integ.t >= _landing_floor(T) && return stop + integ.t < stop < tf && return stop end return tf end @@ -85,11 +82,16 @@ end savedexactly = true while integrator.cur_t <= length(saveat) && saveat[integrator.cur_t] <= integrator.t savet = saveat[integrator.cur_t] - Θ = (savet - integrator.tprev) / integrator.dt - @inbounds us[integrator.cur_t] = _ode_interpolant( - Θ, integrator.dt, - integrator.uprev, integrator - ) + # The step's endpoints are known exactly; only interior points need the + # dense output, which is not exact at Θ = 0 or 1 for every method. + @inbounds us[integrator.cur_t] = if savet == integrator.t + integrator.u + elseif savet == integrator.tprev + integrator.uprev + else + Θ = (savet - integrator.tprev) / integrator.dt + _ode_interpolant(Θ, integrator.dt, integrator.uprev, integrator) + end @inbounds ts[integrator.cur_t] = savet integrator.cur_t += 1 end diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl index cf98e361..0780a29f 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl @@ -126,7 +126,7 @@ end while EEst > convert(T, 1.0) (dt < convert(T, 1.0f-14) || t + dt == t) && - (dt != _next_stop(integ, tf, T) - t || dt < _landing_floor(T)) && + dt != _next_stop(integ, tf, T) - t && error("dt convert(T, 1.0) (dt < convert(T, 1.0f-14) || t + dt == t) && - (dt != _next_stop(integ, tf, T) - t || dt < _landing_floor(T)) && + dt != _next_stop(integ, tf, T) - t && error("dt convert(T, 1.0) (dt < convert(T, 1.0f-14) || t + dt == t) && - (dt != _next_stop(integ, tf, T) - t || dt < _landing_floor(T)) && + dt != _next_stop(integ, tf, T) - t && error("dt convert(T, 1.0) (dt < convert(T, 1.0f-14) || t + dt == t) && - (dt != _next_stop(integ, tf, T) - t || dt < _landing_floor(T)) && + dt != _next_stop(integ, tf, T) - t && error("dt convert(T, 1.0) (dt < convert(T, 1.0f-14) || t + dt == t) && - (dt != _next_stop(integ, tf, T) - t || dt < _landing_floor(T)) && + dt != _next_stop(integ, tf, T) - t && error("dt T(1.0) (dt < T(1.0e-14) || t + dt == t) && - (dt != _next_stop(integ, tf, T) - t || dt < _landing_floor(T)) && + dt != _next_stop(integ, tf, T) - t && error("dt convert(T, 1.0) (dt < convert(T, 1.0f-14) || t + dt == t) && - (dt != _next_stop(integ, tf, T) - t || dt < _landing_floor(T)) && + dt != _next_stop(integ, tf, T) - t && error("dt convert(T, 1.0) (dt < convert(T, 1.0f-14) || t + dt == t) && - (dt != _next_stop(integ, tf, T) - t || dt < _landing_floor(T)) && + dt != _next_stop(integ, tf, T) - t && error("dt SVector(rate, rate); mass_matrix = mass * SMatrix{2, 2, T}(I) + ) + prob = ODEProblem{false}(f, zero(SVector{2, T}), (zero(T), T(2.0e-14))) + cb = DiscreteCallback( + (u, t, integrator) -> t == stop, terminate!; save_positions = (false, false) + ) + sol = solve( + EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0); + trajectories = 2, adaptive = true, dt = T(1.0e-14), tstops = [stop], callback = cb, + merge_callbacks = true, save_everystep = false, abstol = T(1.0e-12), reltol = T(1.0e-3) + ) + exact = BigFloat(rate) * BigFloat(stop) / BigFloat(mass) + @test all( + s -> s.t[end] == stop && + all(x -> isapprox(BigFloat(x), exact; rtol = 128eps(T)), s.u[end]), + sol.u + ) +end + +# `saveat` points at the ends of an accepted step are the step's own states; the dense +# output is only used strictly inside the step. +@testset "Adaptive saveat at step endpoints ($(nameof(typeof(alg))), $T)" for + alg in ADAPTIVE_ALGS, T in (Float32, Float64) + stop = T(0.125) + cb = DiscreteCallback( + (u, t, integrator) -> t == stop, terminate!; save_positions = (false, false) + ) + prob = ODEProblem{false}( + (u, p, t) -> SVector(one(T), one(T)), zero(SVector{2, T}), (zero(T), one(T)) + ) + sol = solve( + EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0); + trajectories = 2, adaptive = true, dt = T(0.5), tstops = [stop], callback = cb, + merge_callbacks = true, saveat = T[0, stop / 2, stop], abstol = T(1.0e-12), + reltol = T(1.0e-6) + ) + @test all(s -> s.t[end] == stop && s.u[1] == zero(SVector{2, T}), sol.u) + @test all(s -> all(x -> isapprox(x, stop; rtol = 128eps(T)), s.u[end]), sol.u) +end From fa8c70f0e87e5bb8a8be3a3eb69a6898b0982f22 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Fri, 2 Oct 2026 23:04:01 -0400 Subject: [PATCH 11/17] Retry an overflowing landing step without `continue` `continue` in the error-control loop produced SPIR-V that fails validation on OpenCL (a block emitted before its dominator). Fold the retry into the rejection branch instead: a non-finite step sets `EEst = Inf`, and a rejected landing step retries with the controller's covering step rather than a shrunken one. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 5.5 (1M context) Agent-Harness: Claude Code 2.1.280 Agent-Model: claude-opus-5-5[1m] Agent-Session: https://claude.ai/code/session_01SJfy3ez7u4CqvAerKrT5np --- .../perform_step/gpu_kvaerno3_perform_step.jl | 14 +++++--------- .../perform_step/gpu_kvaerno5_perform_step.jl | 14 +++++--------- .../perform_step/gpu_rodas4_perform_step.jl | 14 +++++--------- .../perform_step/gpu_rodas5P_perform_step.jl | 14 +++++--------- .../perform_step/gpu_rosenbrock23_perform_step.jl | 14 +++++--------- .../perform_step/gpu_tsit5_perform_step.jl | 14 +++++--------- .../perform_step/gpu_vern7_perform_step.jl | 14 +++++--------- .../perform_step/gpu_vern9_perform_step.jl | 14 +++++--------- 8 files changed, 40 insertions(+), 72 deletions(-) diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl index 0780a29f..1c6b868b 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl @@ -187,14 +187,8 @@ end (abstol .+ max.(abs.(uprev), abs.(u)) * reltol) EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t) # `ODE_DEFAULT_NORM` is computed with fast-math, so NaN is checked on its inputs. - if !(all(isfinite, u) && all(isfinite, tmp)) - EEst = T(Inf) - if shortened_to_land - dt = dtprop - shortened_to_land = false - continue - end - end + finite = mapreduce(isfinite, &, u) & mapreduce(isfinite, &, tmp) + finite || (EEst = T(Inf)) q11 = EEst^beta1 if iszero(EEst) @@ -204,7 +198,9 @@ end end if EEst > 1 - dt = dt / min(inv(qmin), q11 / gamma) + # A landing step that overflowed retries with the controller's covering step, + # which lands on the stop by interpolation. + dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma) shortened_to_land = false else # EEst <= 1 q = max(inv(qmax), min(inv(qmin), q / gamma)) diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl index 6ca325f7..ddeadfff 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl @@ -255,14 +255,8 @@ end (abstol .+ max.(abs.(uprev), abs.(u)) * reltol) EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t) # `ODE_DEFAULT_NORM` is computed with fast-math, so NaN is checked on its inputs. - if !(all(isfinite, u) && all(isfinite, tmp)) - EEst = T(Inf) - if shortened_to_land - dt = dtprop - shortened_to_land = false - continue - end - end + finite = mapreduce(isfinite, &, u) & mapreduce(isfinite, &, tmp) + finite || (EEst = T(Inf)) q11 = EEst^beta1 if iszero(EEst) @@ -272,7 +266,9 @@ end end if EEst > 1 - dt = dt / min(inv(qmin), q11 / gamma) + # A landing step that overflowed retries with the controller's covering step, + # which lands on the stop by interpolation. + dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma) shortened_to_land = false else # EEst <= 1 q = max(inv(qmax), min(inv(qmin), q / gamma)) diff --git a/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl index d6f88239..74409455 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl @@ -231,14 +231,8 @@ end tmp = k6 ./ (abstol .+ max.(abs.(uprev), abs.(u)) * reltol) EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t) # `ODE_DEFAULT_NORM` is computed with fast-math, so NaN is checked on its inputs. - if !(all(isfinite, u) && all(isfinite, tmp)) - EEst = T(Inf) - if shortened_to_land - dt = dtprop - shortened_to_land = false - continue - end - end + finite = mapreduce(isfinite, &, u) & mapreduce(isfinite, &, tmp) + finite || (EEst = T(Inf)) q11 = EEst^beta1 if iszero(EEst) @@ -248,7 +242,9 @@ end end if EEst > 1 - dt = dt / min(inv(qmin), q11 / gamma) + # A landing step that overflowed retries with the controller's covering step, + # which lands on the stop by interpolation. + dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma) shortened_to_land = false else # EEst <= 1 q = max(inv(qmax), min(inv(qmin), q / gamma)) diff --git a/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl index ffaf95c9..2c1fc504 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl @@ -300,14 +300,8 @@ end tmp = k8 ./ (abstol .+ max.(abs.(uprev), abs.(u)) * reltol) EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t) # `ODE_DEFAULT_NORM` is computed with fast-math, so NaN is checked on its inputs. - if !(all(isfinite, u) && all(isfinite, tmp)) - EEst = T(Inf) - if shortened_to_land - dt = dtprop - shortened_to_land = false - continue - end - end + finite = mapreduce(isfinite, &, u) & mapreduce(isfinite, &, tmp) + finite || (EEst = T(Inf)) q11 = EEst^beta1 if iszero(EEst) @@ -317,7 +311,9 @@ end end if EEst > 1 - dt = dt / min(inv(qmin), q11 / gamma) + # A landing step that overflowed retries with the controller's covering step, + # which lands on the stop by interpolation. + dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma) shortened_to_land = false else # EEst <= 1 q = max(inv(qmax), min(inv(qmin), q / gamma)) diff --git a/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl index d2a87c05..5ca1feb2 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl @@ -157,14 +157,8 @@ end tmp = tmp ./ (abstol .+ max.(abs.(uprev), abs.(u)) * reltol) EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t) # `ODE_DEFAULT_NORM` is computed with fast-math, so NaN is checked on its inputs. - if !(all(isfinite, u) && all(isfinite, tmp)) - EEst = T(Inf) - if shortened_to_land - dt = dtprop - shortened_to_land = false - continue - end - end + finite = mapreduce(isfinite, &, u) & mapreduce(isfinite, &, tmp) + finite || (EEst = T(Inf)) q11 = EEst^beta1 if iszero(EEst) @@ -174,7 +168,9 @@ end end if EEst > 1 - dt = dt / min(inv(qmin), q11 / gamma) + # A landing step that overflowed retries with the controller's covering step, + # which lands on the stop by interpolation. + dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma) shortened_to_land = false else # EEst <= 1 q = max(inv(qmax), min(inv(qmin), q / gamma)) diff --git a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl index 98d4def7..70671ecc 100644 --- a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl @@ -126,14 +126,8 @@ end tmp = tmp ./ (abstol .+ max.(abs.(uprev), abs.(u)) * reltol) EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t) # `ODE_DEFAULT_NORM` is computed with fast-math, so NaN is checked on its inputs. - if !(all(isfinite, u) && all(isfinite, tmp)) - EEst = T(Inf) - if shortened_to_land - dt = dtprop - shortened_to_land = false - continue - end - end + finite = mapreduce(isfinite, &, u) & mapreduce(isfinite, &, tmp) + finite || (EEst = T(Inf)) q11 = EEst^beta1 if iszero(EEst) @@ -143,7 +137,9 @@ end end if EEst > 1 - dt = dt / min(inv(qmin), q11 / gamma) + # A landing step that overflowed retries with the controller's covering step, + # which lands on the stop by interpolation. + dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma) shortened_to_land = false else # EEst <= 1 q = max(inv(qmax), min(inv(qmin), q / gamma)) diff --git a/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl index e02e43c5..aa8f718a 100644 --- a/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl @@ -159,14 +159,8 @@ end tmp = tmp ./ (abstol .+ max.(abs.(uprev), abs.(u)) * reltol) EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t) # `ODE_DEFAULT_NORM` is computed with fast-math, so NaN is checked on its inputs. - if !(all(isfinite, u) && all(isfinite, tmp)) - EEst = T(Inf) - if shortened_to_land - dt = dtprop - shortened_to_land = false - continue - end - end + finite = mapreduce(isfinite, &, u) & mapreduce(isfinite, &, tmp) + finite || (EEst = T(Inf)) q11 = EEst^beta1 if iszero(EEst) @@ -176,7 +170,9 @@ end end if EEst > 1 - dt = dt / min(inv(qmin), q11 / gamma) + # A landing step that overflowed retries with the controller's covering step, + # which lands on the stop by interpolation. + dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma) shortened_to_land = false else # EEst <= 1 q = max(inv(qmax), min(inv(qmin), q / gamma)) diff --git a/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl index 0647c583..f87fcd10 100644 --- a/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl @@ -268,14 +268,8 @@ end tmp = tmp ./ (abstol .+ max.(abs.(uprev), abs.(u)) * reltol) EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t) # `ODE_DEFAULT_NORM` is computed with fast-math, so NaN is checked on its inputs. - if !(all(isfinite, u) && all(isfinite, tmp)) - EEst = T(Inf) - if shortened_to_land - dt = dtprop - shortened_to_land = false - continue - end - end + finite = mapreduce(isfinite, &, u) & mapreduce(isfinite, &, tmp) + finite || (EEst = T(Inf)) q11 = EEst^beta1 if iszero(EEst) @@ -285,7 +279,9 @@ end end if EEst > 1 - dt = dt / min(inv(qmin), q11 / gamma) + # A landing step that overflowed retries with the controller's covering step, + # which lands on the stop by interpolation. + dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma) shortened_to_land = false else # EEst <= 1 q = max(inv(qmax), min(inv(qmin), q / gamma)) From 142fc7cc224eb61dfd65107bc2cef78b5199442c Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 3 Oct 2026 01:10:08 -0400 Subject: [PATCH 12/17] Land a covering retry that reaches past tf on tf The retry after an overflowing landing step uses the controller's covering step, which can reach past `tf`; the step was then accepted past `tf` with the wrong time and state. Such a step now lands on `tf` by interpolation, so endpoint callbacks still run. A state obtained by interpolating to a stop or to `tf` that is not finite fails explicitly instead of returning success. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 5.5 (1M context) Agent-Harness: Claude Code 2.1.280 Agent-Model: claude-opus-5-5[1m] Agent-Session: https://claude.ai/code/session_01SJfy3ez7u4CqvAerKrT5np --- .../perform_step/gpu_kvaerno3_perform_step.jl | 6 ++ .../perform_step/gpu_kvaerno5_perform_step.jl | 6 ++ .../perform_step/gpu_rodas4_perform_step.jl | 6 ++ .../perform_step/gpu_rodas5P_perform_step.jl | 6 ++ .../gpu_rosenbrock23_perform_step.jl | 6 ++ .../perform_step/gpu_tsit5_perform_step.jl | 6 ++ .../perform_step/gpu_vern7_perform_step.jl | 6 ++ .../perform_step/gpu_vern9_perform_step.jl | 6 ++ test/gpu_kernel_de/adaptive_endpoint.jl | 55 +++++++++++++++++++ 9 files changed, 103 insertions(+) diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl index 1c6b868b..f4f8fda5 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl @@ -222,11 +222,17 @@ end integ.t = integ.tstops[integ.tstops_idx] if integ.t - t != dt integ.u = integ(integ.t) + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") end dt = integ.t - integ.tprev integ.tstops_idx += 1 elseif dt == tf - t integ.t = tf + elseif dt > tf - t + # Only a covering retry reaches past `tf`; it lands on `tf` by interpolation. + integ.t = tf + integ.u = integ(tf) + mapreduce(isfinite, &, integ.u) || error("non-finite state at tf") else ##Advance the integrator integ.t += dt diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl index ddeadfff..da2adf16 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl @@ -290,11 +290,17 @@ end integ.t = integ.tstops[integ.tstops_idx] if integ.t - t != dt integ.u = integ(integ.t) + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") end dt = integ.t - integ.tprev integ.tstops_idx += 1 elseif dt == tf - t integ.t = tf + elseif dt > tf - t + # Only a covering retry reaches past `tf`; it lands on `tf` by interpolation. + integ.t = tf + integ.u = integ(tf) + mapreduce(isfinite, &, integ.u) || error("non-finite state at tf") else ##Advance the integrator integ.t += dt diff --git a/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl index 74409455..373d16e9 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl @@ -267,11 +267,17 @@ end integ.t = integ.tstops[integ.tstops_idx] if integ.t - t != dt integ.u = integ(integ.t) + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") end dt = integ.t - integ.tprev integ.tstops_idx += 1 elseif dt == tf - t integ.t = tf + elseif dt > tf - t + # Only a covering retry reaches past `tf`; it lands on `tf` by interpolation. + integ.t = tf + integ.u = integ(tf) + mapreduce(isfinite, &, integ.u) || error("non-finite state at tf") else ##Advance the integrator integ.t += dt diff --git a/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl index 2c1fc504..77a74d13 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl @@ -341,11 +341,17 @@ end integ.t = integ.tstops[integ.tstops_idx] if integ.t - t != dt integ.u = integ(integ.t) + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") end dt = integ.t - integ.tprev integ.tstops_idx += 1 elseif dt == tf - t integ.t = tf + elseif dt > tf - t + # Only a covering retry reaches past `tf`; it lands on `tf` by interpolation. + integ.t = tf + integ.u = integ(tf) + mapreduce(isfinite, &, integ.u) || error("non-finite state at tf") else ##Advance the integrator integ.t += dt diff --git a/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl index 5ca1feb2..8f9487af 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl @@ -192,11 +192,17 @@ end integ.t = integ.tstops[integ.tstops_idx] if integ.t - t != dt integ.u = integ(integ.t) + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") end dt = integ.t - integ.tprev integ.tstops_idx += 1 elseif dt == tf - t integ.t = tf + elseif dt > tf - t + # Only a covering retry reaches past `tf`; it lands on `tf` by interpolation. + integ.t = tf + integ.u = integ(tf) + mapreduce(isfinite, &, integ.u) || error("non-finite state at tf") else ##Advance the integrator integ.t += dt diff --git a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl index 70671ecc..7db8d658 100644 --- a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl @@ -166,11 +166,17 @@ end integ.t = integ.tstops[integ.tstops_idx] if integ.t - t != dt integ.u = integ(integ.t) + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") end dt = integ.t - integ.tprev integ.tstops_idx += 1 elseif dt == tf - t integ.t = tf + elseif dt > tf - t + # Only a covering retry reaches past `tf`; it lands on `tf` by interpolation. + integ.t = tf + integ.u = integ(tf) + mapreduce(isfinite, &, integ.u) || error("non-finite state at tf") else ##Advance the integrator integ.t += dt diff --git a/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl index aa8f718a..050b75dd 100644 --- a/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl @@ -202,11 +202,17 @@ end integ.t = integ.tstops[integ.tstops_idx] if integ.t - t != dt integ.u = integ(integ.t) + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") end dt = integ.t - integ.tprev integ.tstops_idx += 1 elseif dt == tf - t integ.t = tf + elseif dt > tf - t + # Only a covering retry reaches past `tf`; it lands on `tf` by interpolation. + integ.t = tf + integ.u = integ(tf) + mapreduce(isfinite, &, integ.u) || error("non-finite state at tf") else ##Advance the integrator integ.t += dt diff --git a/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl index f87fcd10..199cbd39 100644 --- a/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl @@ -311,11 +311,17 @@ end integ.t = integ.tstops[integ.tstops_idx] if integ.t - t != dt integ.u = integ(integ.t) + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") end dt = integ.t - integ.tprev integ.tstops_idx += 1 elseif dt == tf - t integ.t = tf + elseif dt > tf - t + # Only a covering retry reaches past `tf`; it lands on `tf` by interpolation. + integ.t = tf + integ.u = integ(tf) + mapreduce(isfinite, &, integ.u) || error("non-finite state at tf") else ##Advance the integrator integ.t += dt diff --git a/test/gpu_kernel_de/adaptive_endpoint.jl b/test/gpu_kernel_de/adaptive_endpoint.jl index 9bdffa8c..5de415b8 100644 --- a/test/gpu_kernel_de/adaptive_endpoint.jl +++ b/test/gpu_kernel_de/adaptive_endpoint.jl @@ -456,3 +456,58 @@ end @test all(s -> s.t[end] == stop && s.u[1] == zero(SVector{2, T}), sol.u) @test all(s -> all(x -> isapprox(x, stop; rtol = 128eps(T)), s.u[end]), sol.u) end + +# A final landing step that overflows retries with the controller's covering step, which +# may reach past `tf`; it lands on `tf` by interpolation and runs the endpoint callbacks. +@testset "Adaptive overflowing final landing step ($(nameof(typeof(alg))), $T, $mode)" for + alg in (GPURosenbrock23(), GPURodas4(), GPURodas5P()), T in (Float32, Float64), + mode in (:tiny, :scaled, :scaled_callback) + mass = mode == :tiny ? one(T) : T === Float32 ? T(1.0e30) : T(1.0e300) + rate = mode == :tiny ? T(1.0e14) : mass + tf = mode == :tiny ? (T === Float32 ? T(5.0e-38) : T(1.0e-307)) : T(1.0e-10) + f = ODEFunction{false}( + (u, p, t) -> SVector(rate, rate); mass_matrix = mass * SMatrix{2, 2, T}(I) + ) + prob = ODEProblem{false}(f, zero(SVector{2, T}), (zero(T), tf)) + cb = DiscreteCallback( + (u, t, integrator) -> t == tf, integrator -> (integrator.u += SVector(one(T), one(T))); + save_positions = (false, false) + ) + kw = mode == :scaled_callback ? (; tstops = T[tf], callback = cb, merge_callbacks = true) : + (;) + sol = solve( + EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0); + trajectories = 2, adaptive = true, dt = mode == :tiny ? T(1.0e-14) : T(0.001), + save_everystep = false, abstol = T(1.0e-12), reltol = T(1.0e-3), kw... + ) + exact = BigFloat(rate) * BigFloat(tf) / BigFloat(mass) + (mode == :scaled_callback) + @test all( + s -> s.t[end] == tf && + all(x -> isapprox(BigFloat(x), exact; rtol = 128eps(T)), s.u[end]), + sol.u + ) +end + +# A right-hand side that is NaN exactly at `tf` must not yield a successful solve past +# `tf` or a non-finite state: the solve either ends at `tf` with an accurate state or +# fails explicitly. +@testset "Adaptive NaN right-hand side at tf ($(nameof(typeof(alg))), $T)" for + alg in ADAPTIVE_ALGS, T in (Float32, Float64) + tf = T(0.125) + prob = ODEProblem{false}( + (u, p, t) -> t == tf ? SVector(T(NaN), T(NaN)) : SVector(one(T), one(T)), + zero(SVector{2, T}), (zero(T), tf) + ) + outcome = try + sol = solve( + EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0); + trajectories = 2, adaptive = true, dt = T(0.5), save_everystep = false, + abstol = T(1.0e-12), reltol = T(1.0e-3) + ) + all(s -> s.t[end] == tf && all(x -> isapprox(x, tf; rtol = 1.0e-3), s.u[end]), sol.u) + catch err + err isa ErrorException && + (occursin("dt Date: Sat, 3 Oct 2026 02:33:18 -0400 Subject: [PATCH 13/17] Land covered stops and tf through one interpolation site Fold the covering retry's landing on `tf` into the stop-landing branch so each adaptive stepper keeps a single dense-output evaluation in the accept branch, which keeps the kernels' code size (and Julia 1.12 compile time) where it was. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 5.5 (1M context) Agent-Harness: Claude Code 2.1.280 Agent-Model: claude-opus-5-5[1m] Agent-Session: https://claude.ai/code/session_01SJfy3ez7u4CqvAerKrT5np --- .../perform_step/gpu_kvaerno3_perform_step.jl | 16 +++++++--------- .../perform_step/gpu_kvaerno5_perform_step.jl | 16 +++++++--------- .../perform_step/gpu_rodas4_perform_step.jl | 16 +++++++--------- .../perform_step/gpu_rodas5P_perform_step.jl | 16 +++++++--------- .../gpu_rosenbrock23_perform_step.jl | 16 +++++++--------- .../perform_step/gpu_tsit5_perform_step.jl | 16 +++++++--------- .../perform_step/gpu_vern7_perform_step.jl | 16 +++++++--------- .../perform_step/gpu_vern9_perform_step.jl | 16 +++++++--------- 8 files changed, 56 insertions(+), 72 deletions(-) diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl index f4f8fda5..a28bf01a 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl @@ -218,21 +218,19 @@ end integ.tprev = t integ.u = u - if _tstop_in_step(integ, tf, T) - integ.t = integ.tstops[integ.tstops_idx] + stop_in_step = _tstop_in_step(integ, tf, T) + # Only a covering retry reaches past `tf`; like a covered stop, it lands by + # interpolation. + if stop_in_step || dt > tf - t + integ.t = stop_in_step ? integ.tstops[integ.tstops_idx] : tf if integ.t - t != dt integ.u = integ(integ.t) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop or tf") end dt = integ.t - integ.tprev - integ.tstops_idx += 1 + stop_in_step && (integ.tstops_idx += 1) elseif dt == tf - t integ.t = tf - elseif dt > tf - t - # Only a covering retry reaches past `tf`; it lands on `tf` by interpolation. - integ.t = tf - integ.u = integ(tf) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tf") else ##Advance the integrator integ.t += dt diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl index da2adf16..5db25220 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl @@ -286,21 +286,19 @@ end integ.tprev = t integ.u = u - if _tstop_in_step(integ, tf, T) - integ.t = integ.tstops[integ.tstops_idx] + stop_in_step = _tstop_in_step(integ, tf, T) + # Only a covering retry reaches past `tf`; like a covered stop, it lands by + # interpolation. + if stop_in_step || dt > tf - t + integ.t = stop_in_step ? integ.tstops[integ.tstops_idx] : tf if integ.t - t != dt integ.u = integ(integ.t) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop or tf") end dt = integ.t - integ.tprev - integ.tstops_idx += 1 + stop_in_step && (integ.tstops_idx += 1) elseif dt == tf - t integ.t = tf - elseif dt > tf - t - # Only a covering retry reaches past `tf`; it lands on `tf` by interpolation. - integ.t = tf - integ.u = integ(tf) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tf") else ##Advance the integrator integ.t += dt diff --git a/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl index 373d16e9..23aa8a46 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl @@ -263,21 +263,19 @@ end integ.tprev = t integ.u = u - if _tstop_in_step(integ, tf, T) - integ.t = integ.tstops[integ.tstops_idx] + stop_in_step = _tstop_in_step(integ, tf, T) + # Only a covering retry reaches past `tf`; like a covered stop, it lands by + # interpolation. + if stop_in_step || dt > tf - t + integ.t = stop_in_step ? integ.tstops[integ.tstops_idx] : tf if integ.t - t != dt integ.u = integ(integ.t) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop or tf") end dt = integ.t - integ.tprev - integ.tstops_idx += 1 + stop_in_step && (integ.tstops_idx += 1) elseif dt == tf - t integ.t = tf - elseif dt > tf - t - # Only a covering retry reaches past `tf`; it lands on `tf` by interpolation. - integ.t = tf - integ.u = integ(tf) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tf") else ##Advance the integrator integ.t += dt diff --git a/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl index 77a74d13..d5cecd96 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl @@ -337,21 +337,19 @@ end integ.tprev = t integ.u = u - if _tstop_in_step(integ, tf, T) - integ.t = integ.tstops[integ.tstops_idx] + stop_in_step = _tstop_in_step(integ, tf, T) + # Only a covering retry reaches past `tf`; like a covered stop, it lands by + # interpolation. + if stop_in_step || dt > tf - t + integ.t = stop_in_step ? integ.tstops[integ.tstops_idx] : tf if integ.t - t != dt integ.u = integ(integ.t) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop or tf") end dt = integ.t - integ.tprev - integ.tstops_idx += 1 + stop_in_step && (integ.tstops_idx += 1) elseif dt == tf - t integ.t = tf - elseif dt > tf - t - # Only a covering retry reaches past `tf`; it lands on `tf` by interpolation. - integ.t = tf - integ.u = integ(tf) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tf") else ##Advance the integrator integ.t += dt diff --git a/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl index 8f9487af..3f3a334d 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl @@ -188,21 +188,19 @@ end integ.tprev = t integ.u = u - if _tstop_in_step(integ, tf, T) - integ.t = integ.tstops[integ.tstops_idx] + stop_in_step = _tstop_in_step(integ, tf, T) + # Only a covering retry reaches past `tf`; like a covered stop, it lands by + # interpolation. + if stop_in_step || dt > tf - t + integ.t = stop_in_step ? integ.tstops[integ.tstops_idx] : tf if integ.t - t != dt integ.u = integ(integ.t) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop or tf") end dt = integ.t - integ.tprev - integ.tstops_idx += 1 + stop_in_step && (integ.tstops_idx += 1) elseif dt == tf - t integ.t = tf - elseif dt > tf - t - # Only a covering retry reaches past `tf`; it lands on `tf` by interpolation. - integ.t = tf - integ.u = integ(tf) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tf") else ##Advance the integrator integ.t += dt diff --git a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl index 7db8d658..5409b87e 100644 --- a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl @@ -162,21 +162,19 @@ end integ.tprev = t integ.u = u - if _tstop_in_step(integ, tf, T) - integ.t = integ.tstops[integ.tstops_idx] + stop_in_step = _tstop_in_step(integ, tf, T) + # Only a covering retry reaches past `tf`; like a covered stop, it lands by + # interpolation. + if stop_in_step || dt > tf - t + integ.t = stop_in_step ? integ.tstops[integ.tstops_idx] : tf if integ.t - t != dt integ.u = integ(integ.t) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop or tf") end dt = integ.t - integ.tprev - integ.tstops_idx += 1 + stop_in_step && (integ.tstops_idx += 1) elseif dt == tf - t integ.t = tf - elseif dt > tf - t - # Only a covering retry reaches past `tf`; it lands on `tf` by interpolation. - integ.t = tf - integ.u = integ(tf) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tf") else ##Advance the integrator integ.t += dt diff --git a/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl index 050b75dd..8468e263 100644 --- a/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl @@ -198,21 +198,19 @@ end integ.tprev = t integ.u = u - if _tstop_in_step(integ, tf, T) - integ.t = integ.tstops[integ.tstops_idx] + stop_in_step = _tstop_in_step(integ, tf, T) + # Only a covering retry reaches past `tf`; like a covered stop, it lands by + # interpolation. + if stop_in_step || dt > tf - t + integ.t = stop_in_step ? integ.tstops[integ.tstops_idx] : tf if integ.t - t != dt integ.u = integ(integ.t) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop or tf") end dt = integ.t - integ.tprev - integ.tstops_idx += 1 + stop_in_step && (integ.tstops_idx += 1) elseif dt == tf - t integ.t = tf - elseif dt > tf - t - # Only a covering retry reaches past `tf`; it lands on `tf` by interpolation. - integ.t = tf - integ.u = integ(tf) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tf") else ##Advance the integrator integ.t += dt diff --git a/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl index 199cbd39..9acb1a15 100644 --- a/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl @@ -307,21 +307,19 @@ end integ.tprev = t integ.u = u - if _tstop_in_step(integ, tf, T) - integ.t = integ.tstops[integ.tstops_idx] + stop_in_step = _tstop_in_step(integ, tf, T) + # Only a covering retry reaches past `tf`; like a covered stop, it lands by + # interpolation. + if stop_in_step || dt > tf - t + integ.t = stop_in_step ? integ.tstops[integ.tstops_idx] : tf if integ.t - t != dt integ.u = integ(integ.t) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop or tf") end dt = integ.t - integ.tprev - integ.tstops_idx += 1 + stop_in_step && (integ.tstops_idx += 1) elseif dt == tf - t integ.t = tf - elseif dt > tf - t - # Only a covering retry reaches past `tf`; it lands on `tf` by interpolation. - integ.t = tf - integ.u = integ(tf) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tf") else ##Advance the integrator integ.t += dt From a2e999d60e743a471d3ec85df887e42f2f9eb4a1 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 3 Oct 2026 05:33:22 -0400 Subject: [PATCH 14/17] Retry only stop landings with a covering step A covering retry that could reach past `tf` made the dense-output call in the accept branch reachable for every problem, including those without tstops, where it is otherwise dead code. Compiling the Verner interpolants for ForwardDiff duals inline in the step loop took 236 s (Vern7) and 1177 s (Vern9) on Julia 1.12, against ~20 s and ~50 s on master. The covering retry now applies only to a landing on a stop before `tf`, so it never reaches past `tf`. A final interval whose step is not finite fails explicitly (`dt Co-Authored-By: Claude Opus 5.5 (1M context) Agent-Harness: Claude Code 2.1.280 Agent-Model: claude-opus-5-5[1m] Agent-Session: https://claude.ai/code/session_01SJfy3ez7u4CqvAerKrT5np --- .../perform_step/gpu_kvaerno3_perform_step.jl | 17 +++---- .../perform_step/gpu_kvaerno5_perform_step.jl | 17 +++---- .../perform_step/gpu_rodas4_perform_step.jl | 17 +++---- .../perform_step/gpu_rodas5P_perform_step.jl | 17 +++---- .../gpu_rosenbrock23_perform_step.jl | 17 +++---- .../perform_step/gpu_tsit5_perform_step.jl | 17 +++---- .../perform_step/gpu_vern7_perform_step.jl | 17 +++---- .../perform_step/gpu_vern9_perform_step.jl | 17 +++---- test/gpu_kernel_de/adaptive_endpoint.jl | 44 +++++++++++-------- 9 files changed, 82 insertions(+), 98 deletions(-) diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl index a28bf01a..33b30b60 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl @@ -101,7 +101,7 @@ end tf = integ.tf dtprop = dt dt = _bounded_step(integ, dt, tf, T) - shortened_to_land = dt < dtprop + shortened_to_land = dt < dtprop && _next_stop(integ, tf, T) < tf tmp = integ.tmp f = integ.f @@ -198,8 +198,8 @@ end end if EEst > 1 - # A landing step that overflowed retries with the controller's covering step, - # which lands on the stop by interpolation. + # A landing step on a stop that overflowed retries with the controller's covering + # step, which lands on the stop by interpolation; it never reaches past `tf`. dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma) shortened_to_land = false else # EEst <= 1 @@ -218,17 +218,14 @@ end integ.tprev = t integ.u = u - stop_in_step = _tstop_in_step(integ, tf, T) - # Only a covering retry reaches past `tf`; like a covered stop, it lands by - # interpolation. - if stop_in_step || dt > tf - t - integ.t = stop_in_step ? integ.tstops[integ.tstops_idx] : tf + if _tstop_in_step(integ, tf, T) + integ.t = integ.tstops[integ.tstops_idx] if integ.t - t != dt integ.u = integ(integ.t) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop or tf") + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") end dt = integ.t - integ.tprev - stop_in_step && (integ.tstops_idx += 1) + integ.tstops_idx += 1 elseif dt == tf - t integ.t = tf else diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl index 5db25220..f9d06599 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl @@ -132,7 +132,7 @@ end tf = integ.tf dtprop = dt dt = _bounded_step(integ, dt, tf, T) - shortened_to_land = dt < dtprop + shortened_to_land = dt < dtprop && _next_stop(integ, tf, T) < tf tmp = integ.tmp f = integ.f @@ -266,8 +266,8 @@ end end if EEst > 1 - # A landing step that overflowed retries with the controller's covering step, - # which lands on the stop by interpolation. + # A landing step on a stop that overflowed retries with the controller's covering + # step, which lands on the stop by interpolation; it never reaches past `tf`. dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma) shortened_to_land = false else # EEst <= 1 @@ -286,17 +286,14 @@ end integ.tprev = t integ.u = u - stop_in_step = _tstop_in_step(integ, tf, T) - # Only a covering retry reaches past `tf`; like a covered stop, it lands by - # interpolation. - if stop_in_step || dt > tf - t - integ.t = stop_in_step ? integ.tstops[integ.tstops_idx] : tf + if _tstop_in_step(integ, tf, T) + integ.t = integ.tstops[integ.tstops_idx] if integ.t - t != dt integ.u = integ(integ.t) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop or tf") + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") end dt = integ.t - integ.tprev - stop_in_step && (integ.tstops_idx += 1) + integ.tstops_idx += 1 elseif dt == tf - t integ.t = tf else diff --git a/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl index 23aa8a46..53a012bc 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl @@ -128,7 +128,7 @@ end tf = integ.tf dtprop = dt dt = _bounded_step(integ, dt, tf, T) - shortened_to_land = dt < dtprop + shortened_to_land = dt < dtprop && _next_stop(integ, tf, T) < tf tmp = integ.tmp f = integ.f @@ -242,8 +242,8 @@ end end if EEst > 1 - # A landing step that overflowed retries with the controller's covering step, - # which lands on the stop by interpolation. + # A landing step on a stop that overflowed retries with the controller's covering + # step, which lands on the stop by interpolation; it never reaches past `tf`. dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma) shortened_to_land = false else # EEst <= 1 @@ -263,17 +263,14 @@ end integ.tprev = t integ.u = u - stop_in_step = _tstop_in_step(integ, tf, T) - # Only a covering retry reaches past `tf`; like a covered stop, it lands by - # interpolation. - if stop_in_step || dt > tf - t - integ.t = stop_in_step ? integ.tstops[integ.tstops_idx] : tf + if _tstop_in_step(integ, tf, T) + integ.t = integ.tstops[integ.tstops_idx] if integ.t - t != dt integ.u = integ(integ.t) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop or tf") + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") end dt = integ.t - integ.tprev - stop_in_step && (integ.tstops_idx += 1) + integ.tstops_idx += 1 elseif dt == tf - t integ.t = tf else diff --git a/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl index d5cecd96..c8bae065 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl @@ -165,7 +165,7 @@ end tf = integ.tf dtprop = dt dt = _bounded_step(integ, dt, tf, T) - shortened_to_land = dt < dtprop + shortened_to_land = dt < dtprop && _next_stop(integ, tf, T) < tf tmp = integ.tmp f = integ.f @@ -311,8 +311,8 @@ end end if EEst > 1 - # A landing step that overflowed retries with the controller's covering step, - # which lands on the stop by interpolation. + # A landing step on a stop that overflowed retries with the controller's covering + # step, which lands on the stop by interpolation; it never reaches past `tf`. dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma) shortened_to_land = false else # EEst <= 1 @@ -337,17 +337,14 @@ end integ.tprev = t integ.u = u - stop_in_step = _tstop_in_step(integ, tf, T) - # Only a covering retry reaches past `tf`; like a covered stop, it lands by - # interpolation. - if stop_in_step || dt > tf - t - integ.t = stop_in_step ? integ.tstops[integ.tstops_idx] : tf + if _tstop_in_step(integ, tf, T) + integ.t = integ.tstops[integ.tstops_idx] if integ.t - t != dt integ.u = integ(integ.t) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop or tf") + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") end dt = integ.t - integ.tprev - stop_in_step && (integ.tstops_idx += 1) + integ.tstops_idx += 1 elseif dt == tf - t integ.t = tf else diff --git a/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl index 3f3a334d..31480bf8 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl @@ -83,7 +83,7 @@ end tf = integ.tf dtprop = dt dt = _bounded_step(integ, dt, tf, T) - shortened_to_land = dt < dtprop + shortened_to_land = dt < dtprop && _next_stop(integ, tf, T) < tf tmp = integ.tmp f = integ.f @@ -168,8 +168,8 @@ end end if EEst > 1 - # A landing step that overflowed retries with the controller's covering step, - # which lands on the stop by interpolation. + # A landing step on a stop that overflowed retries with the controller's covering + # step, which lands on the stop by interpolation; it never reaches past `tf`. dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma) shortened_to_land = false else # EEst <= 1 @@ -188,17 +188,14 @@ end integ.tprev = t integ.u = u - stop_in_step = _tstop_in_step(integ, tf, T) - # Only a covering retry reaches past `tf`; like a covered stop, it lands by - # interpolation. - if stop_in_step || dt > tf - t - integ.t = stop_in_step ? integ.tstops[integ.tstops_idx] : tf + if _tstop_in_step(integ, tf, T) + integ.t = integ.tstops[integ.tstops_idx] if integ.t - t != dt integ.u = integ(integ.t) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop or tf") + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") end dt = integ.t - integ.tprev - stop_in_step && (integ.tstops_idx += 1) + integ.tstops_idx += 1 elseif dt == tf - t integ.t = tf else diff --git a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl index 5409b87e..4a7f05b9 100644 --- a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl @@ -78,7 +78,7 @@ end tf = integ.tf dtprop = dt dt = _bounded_step(integ, dt, tf, T) - shortened_to_land = dt < dtprop + shortened_to_land = dt < dtprop && _next_stop(integ, tf, T) < tf a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as btilde1, btilde2, btilde3, btilde4, btilde5, btilde6, btilde7 = integ.btildes @@ -137,8 +137,8 @@ end end if EEst > 1 - # A landing step that overflowed retries with the controller's covering step, - # which lands on the stop by interpolation. + # A landing step on a stop that overflowed retries with the controller's covering + # step, which lands on the stop by interpolation; it never reaches past `tf`. dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma) shortened_to_land = false else # EEst <= 1 @@ -162,17 +162,14 @@ end integ.tprev = t integ.u = u - stop_in_step = _tstop_in_step(integ, tf, T) - # Only a covering retry reaches past `tf`; like a covered stop, it lands by - # interpolation. - if stop_in_step || dt > tf - t - integ.t = stop_in_step ? integ.tstops[integ.tstops_idx] : tf + if _tstop_in_step(integ, tf, T) + integ.t = integ.tstops[integ.tstops_idx] if integ.t - t != dt integ.u = integ(integ.t) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop or tf") + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") end dt = integ.t - integ.tprev - stop_in_step && (integ.tstops_idx += 1) + integ.tstops_idx += 1 elseif dt == tf - t integ.t = tf else diff --git a/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl index 8468e263..cbd1ce08 100644 --- a/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl @@ -88,7 +88,7 @@ end tf = integ.tf dtprop = dt dt = _bounded_step(integ, dt, tf, T) - shortened_to_land = dt < dtprop + shortened_to_land = dt < dtprop && _next_stop(integ, tf, T) < tf @unpack c2, c3, c4, c5, c6, c7, c8, a021, a031, a032, a041, a043, a051, a053, a054, a061, a063, a064, a065, a071, a073, a074, a075, a076, a081, a083, a084, @@ -170,8 +170,8 @@ end end if EEst > 1 - # A landing step that overflowed retries with the controller's covering step, - # which lands on the stop by interpolation. + # A landing step on a stop that overflowed retries with the controller's covering + # step, which lands on the stop by interpolation; it never reaches past `tf`. dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma) shortened_to_land = false else # EEst <= 1 @@ -198,17 +198,14 @@ end integ.tprev = t integ.u = u - stop_in_step = _tstop_in_step(integ, tf, T) - # Only a covering retry reaches past `tf`; like a covered stop, it lands by - # interpolation. - if stop_in_step || dt > tf - t - integ.t = stop_in_step ? integ.tstops[integ.tstops_idx] : tf + if _tstop_in_step(integ, tf, T) + integ.t = integ.tstops[integ.tstops_idx] if integ.t - t != dt integ.u = integ(integ.t) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop or tf") + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") end dt = integ.t - integ.tprev - stop_in_step && (integ.tstops_idx += 1) + integ.tstops_idx += 1 elseif dt == tf - t integ.t = tf else diff --git a/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl index 9acb1a15..62ea89ae 100644 --- a/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl @@ -138,7 +138,7 @@ end tf = integ.tf dtprop = dt dt = _bounded_step(integ, dt, tf, T) - shortened_to_land = dt < dtprop + shortened_to_land = dt < dtprop && _next_stop(integ, tf, T) < tf @unpack c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, a0201, a0301, a0302, a0401, a0403, a0501, a0503, a0504, a0601, a0604, a0605, a0701, a0704, a0705, a0706, @@ -279,8 +279,8 @@ end end if EEst > 1 - # A landing step that overflowed retries with the controller's covering step, - # which lands on the stop by interpolation. + # A landing step on a stop that overflowed retries with the controller's covering + # step, which lands on the stop by interpolation; it never reaches past `tf`. dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma) shortened_to_land = false else # EEst <= 1 @@ -307,17 +307,14 @@ end integ.tprev = t integ.u = u - stop_in_step = _tstop_in_step(integ, tf, T) - # Only a covering retry reaches past `tf`; like a covered stop, it lands by - # interpolation. - if stop_in_step || dt > tf - t - integ.t = stop_in_step ? integ.tstops[integ.tstops_idx] : tf + if _tstop_in_step(integ, tf, T) + integ.t = integ.tstops[integ.tstops_idx] if integ.t - t != dt integ.u = integ(integ.t) - mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop or tf") + mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop") end dt = integ.t - integ.tprev - stop_in_step && (integ.tstops_idx += 1) + integ.tstops_idx += 1 elseif dt == tf - t integ.t = tf else diff --git a/test/gpu_kernel_de/adaptive_endpoint.jl b/test/gpu_kernel_de/adaptive_endpoint.jl index 5de415b8..273e650a 100644 --- a/test/gpu_kernel_de/adaptive_endpoint.jl +++ b/test/gpu_kernel_de/adaptive_endpoint.jl @@ -457,8 +457,24 @@ end @test all(s -> all(x -> isapprox(x, stop; rtol = 128eps(T)), s.u[end]), sol.u) end -# A final landing step that overflows retries with the controller's covering step, which -# may reach past `tf`; it lands on `tf` by interpolation and runs the endpoint callbacks. +# A covering retry never reaches past `tf`: when the final interval cannot be integrated +# the solve fails explicitly. The only acceptable outcomes are an accurate state at +# exactly `tf`, or an explicit failure; success at another time, or with a non-finite +# state, is never acceptable. +function endpoint_or_failure(solve_thunk, tf, exact; rtol) + sol = try + solve_thunk() + catch err + return err isa ErrorException && + (occursin("dt isfinite(x) && isapprox(BigFloat(x), exact; rtol), s.u[end]) + end +end + @testset "Adaptive overflowing final landing step ($(nameof(typeof(alg))), $T, $mode)" for alg in (GPURosenbrock23(), GPURodas4(), GPURodas5P()), T in (Float32, Float64), mode in (:tiny, :scaled, :scaled_callback) @@ -475,22 +491,19 @@ end ) kw = mode == :scaled_callback ? (; tstops = T[tf], callback = cb, merge_callbacks = true) : (;) - sol = solve( + solve_final() = solve( EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0); trajectories = 2, adaptive = true, dt = mode == :tiny ? T(1.0e-14) : T(0.001), save_everystep = false, abstol = T(1.0e-12), reltol = T(1.0e-3), kw... ) exact = BigFloat(rate) * BigFloat(tf) / BigFloat(mass) + (mode == :scaled_callback) - @test all( - s -> s.t[end] == tf && - all(x -> isapprox(BigFloat(x), exact; rtol = 128eps(T)), s.u[end]), - sol.u - ) + @test endpoint_or_failure(tf, exact; rtol = 128eps(T)) do + solve_final() + end end -# A right-hand side that is NaN exactly at `tf` must not yield a successful solve past -# `tf` or a non-finite state: the solve either ends at `tf` with an accurate state or -# fails explicitly. +# A right-hand side that is NaN exactly at `tf` ends either accurately at `tf` or in an +# explicit failure. @testset "Adaptive NaN right-hand side at tf ($(nameof(typeof(alg))), $T)" for alg in ADAPTIVE_ALGS, T in (Float32, Float64) tf = T(0.125) @@ -498,16 +511,11 @@ end (u, p, t) -> t == tf ? SVector(T(NaN), T(NaN)) : SVector(one(T), one(T)), zero(SVector{2, T}), (zero(T), tf) ) - outcome = try - sol = solve( + @test endpoint_or_failure(tf, BigFloat(tf); rtol = 1.0e-3) do + solve( EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0); trajectories = 2, adaptive = true, dt = T(0.5), save_everystep = false, abstol = T(1.0e-12), reltol = T(1.0e-3) ) - all(s -> s.t[end] == tf && all(x -> isapprox(x, tf; rtol = 1.0e-3), s.u[end]), sol.u) - catch err - err isa ErrorException && - (occursin("dt Date: Sat, 3 Oct 2026 06:35:42 -0400 Subject: [PATCH 15/17] State the endpoint invariant with the stop-only covering retry Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 5.5 (1M context) Agent-Harness: Claude Code 2.1.280 Agent-Model: claude-opus-5-5[1m] Agent-Session: https://claude.ai/code/session_01SJfy3ez7u4CqvAerKrT5np --- src/ensemblegpukernel/integrators/integrator_utils.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/ensemblegpukernel/integrators/integrator_utils.jl b/src/ensemblegpukernel/integrators/integrator_utils.jl index 61f0bf17..b2ac63a6 100644 --- a/src/ensemblegpukernel/integrators/integrator_utils.jl +++ b/src/ensemblegpukernel/integrators/integrator_utils.jl @@ -16,7 +16,10 @@ end # Adaptive steps land on pending tstops and on `tf` by construction: `_bounded_step` # makes a trial step that would reach (or nearly reach) the next landing target exactly -# the remaining distance, and the accept branch only marks `t = tf` for such a step. +# the remaining distance, and the accept branch marks `t = tf` only for an accepted step +# of length exactly `tf - t`. A step whose state is not finite is never accepted. Only a +# landing on a stop before `tf` retries with the controller's covering step (landing on +# the stop by interpolation), so no accepted step reaches past `tf`. # There is no tolerance-based endpoint snap. (At large |t| a step shortened by error # control can still round onto `tf` in `t += dt`; see # https://github.com/SciML/DiffEqGPU.jl/issues/567.) @@ -35,9 +38,7 @@ end end # The next time the adaptive integrator has to land on exactly: the next pending -# tstop after `integ.t` and before `tf`, otherwise `tf`. A landing step whose arithmetic -# overflows (non-finite error estimate) is retried with the controller's unshortened -# step, which covers the stop and lands on it by interpolation. +# tstop after `integ.t` and before `tf`, otherwise `tf`. @inline function _next_stop(integ, tf, ::Type{T}) where {T} tstops = integ.tstops if tstops !== nothing && integ.tstops_idx <= length(tstops) From 305f4bb895821c064c3160de3ff2dae0f69a4ca2 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 3 Oct 2026 10:43:42 -0400 Subject: [PATCH 16/17] Build large test times with float powers (32-bit Int overflow) `2^55` is an Int, which overflows to 0 on 32-bit Julia, so the saveat testsets' `eps(t0) == T(8)` assertion failed on the x86 CPU job. `T(2)^55` and `T(2)^26` give bit-identical times on 64-bit. Co-Authored-By: Chris Rackauckas Co-Authored-By: Devin <158243242+devin-ai-integration[bot]@users.noreply.github.com> Co-Authored-By: Claude Opus 5.5 (1M context) Agent-Harness: Devin CLI (edit); Claude Code 2.1.280 (verification and commit) Agent-Model: swe-2; claude-opus-5-5[1m] Agent-Session: local Devin job ~/sandbox/swe2-queue/work/difeqgpu-561-x86-int-overflow-r2 on amdci2; https://claude.ai/code/session_01SJfy3ez7u4CqvAerKrT5np Claude-Session: https://claude.ai/code/session_01SJfy3ez7u4CqvAerKrT5np --- test/gpu_kernel_de/adaptive_endpoint.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/gpu_kernel_de/adaptive_endpoint.jl b/test/gpu_kernel_de/adaptive_endpoint.jl index 273e650a..05bb71cc 100644 --- a/test/gpu_kernel_de/adaptive_endpoint.jl +++ b/test/gpu_kernel_de/adaptive_endpoint.jl @@ -147,7 +147,7 @@ end # a whole number of ulps (see https://github.com/SciML/DiffEqGPU.jl/issues/567). @testset "Adaptive tstop inside a coarse step, public solve ($(nameof(typeof(alg))))" for alg in ADAPTIVE_ALGS - t0 = Float32(2^26) + t0 = Float32(2)^26 tf = t0 + 32.0f0 prob = ODEProblem{false}( (u, p, t) -> SVector(1.0f0, 1.0f0), SVector(0.0f0, 0.0f0), (t0, tf) @@ -231,7 +231,7 @@ end # and is then overwritten with the stop's state would corrupt the sample at t0 + 8. @testset "Adaptive saveat inside a step that reaches a tstop ($(nameof(typeof(alg))), $T)" for alg in ADAPTIVE_ALGS, T in (Float32, Float64) - t0 = T(T === Float32 ? 2^26 : 2^55) + t0 = T === Float32 ? T(2)^26 : T(2)^55 @assert eps(t0) == T(8) prob = ODEProblem{false}( (u, p, t) -> SVector(one(T), one(T)), SVector(zero(T), zero(T)), (t0, t0 + T(32)) From 8adb69ae5bbb42c5ff303383ff4c260b18ed6b18 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Wed, 7 Oct 2026 10:52:27 -0400 Subject: [PATCH 17/17] Save saveat points before a continuous event moves the state back `apply_callback!` moved `integrator.u` to the event state and then saved the pending `saveat` points with the step's dense output. The Rosenbrock (Rodas4, Rodas5P) and Hermite (Kvaerno3, Kvaerno5) interpolants use `integrator.u` as the step's end state, so every `saveat` point between `tprev` and the event was interpolated against the wrong endpoint. Save those points first, while the step's state is intact. The defect predates this branch. Master's clamp of `dtnew` to the distance to `tf` usually kept steps short enough that no `saveat` point shared a step with an event; 1b90199 removed the clamp, so the step after the bounce at t = 9 in the stiff ContinuousCallback test reaches tf, contains the bounce at t = 15, and the save at t = 9.1 was off by 0.71 (GPURodas4, CUDA). Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 5.5 (1M context) Agent-Harness: Claude Code 2.1.280 Agent-Model: claude-opus-5-5 Agent-Session: https://claude.ai/code/session_01SJfy3ez7u4CqvAerKrT5np --- .../integrators/integrator_utils.jl | 40 ++++++++++++------- test/gpu_kernel_de/adaptive_endpoint.jl | 23 +++++++++++ 2 files changed, 48 insertions(+), 15 deletions(-) diff --git a/src/ensemblegpukernel/integrators/integrator_utils.jl b/src/ensemblegpukernel/integrators/integrator_utils.jl index b2ac63a6..bcd9f0fb 100644 --- a/src/ensemblegpukernel/integrators/integrator_utils.jl +++ b/src/ensemblegpukernel/integrators/integrator_utils.jl @@ -81,26 +81,35 @@ end elseif saveat !== nothing saved = true savedexactly = true - while integrator.cur_t <= length(saveat) && saveat[integrator.cur_t] <= integrator.t - savet = saveat[integrator.cur_t] - # The step's endpoints are known exactly; only interior points need the - # dense output, which is not exact at Θ = 0 or 1 for every method. - @inbounds us[integrator.cur_t] = if savet == integrator.t - integrator.u - elseif savet == integrator.tprev - integrator.uprev - else - Θ = (savet - integrator.tprev) / integrator.dt - _ode_interpolant(Θ, integrator.dt, integrator.uprev, integrator) - end - @inbounds ts[integrator.cur_t] = savet - integrator.cur_t += 1 - end + _save_saveat!(integrator, ts, us, integrator.t) end return saved, savedexactly end +# Save the pending `saveat` points up to `tsave` from the current step. The dense output +# of several methods uses `integrator.u` as the step's end state, so this must run +# before `integrator.u` is moved back to an event time. +@inline function _save_saveat!(integrator, ts, us, tsave) + saveat = integrator.saveat + while integrator.cur_t <= length(saveat) && saveat[integrator.cur_t] <= tsave + savet = saveat[integrator.cur_t] + # The step's endpoints are known exactly; only interior points need the + # dense output, which is not exact at Θ = 0 or 1 for every method. + @inbounds us[integrator.cur_t] = if savet == integrator.t + integrator.u + elseif savet == integrator.tprev + integrator.uprev + else + Θ = (savet - integrator.tprev) / integrator.dt + _ode_interpolant(Θ, integrator.dt, integrator.uprev, integrator) + end + @inbounds ts[integrator.cur_t] = savet + integrator.cur_t += 1 + end + return nothing +end + @inline function DiffEqBase.terminate!( integrator::SciMLBase.AbstractODEIntegrator{ AlgType, @@ -287,6 +296,7 @@ end cb_time, prev_sign, event_idx, ts, us ) where {AlgType <: GPUODEAlgorithm, IIP, S, T} + integrator.saveat === nothing || _save_saveat!(integrator, ts, us, cb_time) DiffEqBase.change_t_via_interpolation!(integrator, cb_time) # The new absolute-time callback handling can leave dtnew ≈ 0 when diff --git a/test/gpu_kernel_de/adaptive_endpoint.jl b/test/gpu_kernel_de/adaptive_endpoint.jl index 05bb71cc..9a8666a9 100644 --- a/test/gpu_kernel_de/adaptive_endpoint.jl +++ b/test/gpu_kernel_de/adaptive_endpoint.jl @@ -519,3 +519,26 @@ end ) end end + +# A `saveat` point inside a step that a ContinuousCallback cuts short at an event is +# saved from that step's dense output, whose end state is the step's own, not the state +# moved back to the event. +@testset "Adaptive saveat before a continuous event ($(nameof(typeof(alg))), $T)" for + alg in ADAPTIVE_ALGS, T in (Float32, Float64) + cb = ContinuousCallback( + (u, t, integrator) -> u[1], + integrator -> (integrator.u = SVector(integrator.u[1], -integrator.u[2])); + save_positions = (false, false) + ) + prob = ODEProblem{false}( + (u, p, t) -> SVector(u[2], zero(T)), SVector{2, T}(1, -1), (zero(T), T(2)) + ) + sol = solve( + EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0); + trajectories = 2, adaptive = true, dt = T(2), callback = cb, + merge_callbacks = true, saveat = T[0.5, 2] + ) + @test all(s -> s.t == T[0.5, 2], sol.u) + @test all(s -> isapprox(s.u[1], SVector{2, T}(0.5, -1); atol = 1.0e-4), sol.u) + @test all(s -> isapprox(s.u[2], SVector{2, T}(1, 1); atol = 1.0e-4), sol.u) +end