From e617de46ae28763d6944357e90e712fe0e0a4fc6 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Wed, 23 Sep 2026 04:53:34 -0400 Subject: [PATCH 1/4] Fix Float32 endpoint stall in adaptive GPU kernel steppers The accept branch estimated the remaining interval as `tf - t - dt`, which is evaluated as `(tf - t) - dt` and can differ from the gap left after the rounded update `t + dt`. In Float32 this can leave `t` one ulp short of `tf` with a proposed step of half an ulp; `t + dt` then rounds back to `t`, the absolute `1e-14` snap never fires, and the kernel loops forever. Computing the gap as `tf - (t + dt)` makes the next step land exactly on `tf`. 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 | 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 +++++++++++++++++++ test/runtests.jl | 3 ++ 10 files changed, 51 insertions(+), 16 deletions(-) create mode 100644 test/gpu_kernel_de/adaptive_endpoint.jl diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl index 8843a99a..f1fcdcd6 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl @@ -195,7 +195,7 @@ end q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew - dtnew = min(abs(dtnew), abs(tf - t - dt)) + dtnew = min(abs(dtnew), abs(tf - (t + dt))) @inbounds begin # Necessary for interpolation integ.k1 = k1 @@ -208,7 +208,7 @@ end integ.tprev = t integ.u = u - if (tf - t - dt) < convert(T, 1.0f-14) + if tf - (t + dt) < convert(T, 1.0f-14) integ.t = tf else if integ.tstops !== nothing && integ.tstops_idx <= length(integ.tstops) && diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl index ed67189e..c70a941c 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl @@ -263,7 +263,7 @@ end q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew - dtnew = min(abs(dtnew), abs(tf - t - dt)) + dtnew = min(abs(dtnew), abs(tf - (t + dt))) @inbounds begin # Necessary for interpolation integ.k1 = k1 @@ -276,7 +276,7 @@ end integ.tprev = t integ.u = u - if (tf - t - dt) < convert(T, 1.0f-14) + if tf - (t + dt) < convert(T, 1.0f-14) integ.t = tf else if integ.tstops !== nothing && integ.tstops_idx <= length(integ.tstops) && diff --git a/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl index 819d61ad..c2d72cca 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl @@ -239,7 +239,7 @@ end q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew - dtnew = min(abs(dtnew), abs(tf - t - dt)) + dtnew = min(abs(dtnew), abs(tf - (t + dt))) @inbounds begin # Necessary for interpolation @unpack h21, h22, h23, h24, h25, h31, h32, h33, h34, h35 = integ.tab @@ -253,7 +253,7 @@ end integ.tprev = t integ.u = u - if (tf - t - dt) < convert(T, 1.0f-14) + if tf - (t + dt) < convert(T, 1.0f-14) integ.t = tf else if integ.tstops !== nothing && integ.tstops_idx <= length(integ.tstops) && diff --git a/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl index c86fb74e..01c871eb 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl @@ -308,7 +308,7 @@ end q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew - dtnew = min(abs(dtnew), abs(tf - t - dt)) + dtnew = min(abs(dtnew), abs(tf - (t + dt))) @inbounds begin # Necessary for interpolation @unpack h21, h22, h23, h24, h25, h26, h27, h28, h31, h32, h33, h34, @@ -327,7 +327,7 @@ end integ.tprev = t integ.u = u - if (tf - t - dt) < convert(T, 1.0f-14) + if tf - (t + dt) < convert(T, 1.0f-14) integ.t = tf else if integ.tstops !== nothing && integ.tstops_idx <= length(integ.tstops) && diff --git a/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl index 89187501..df65cb82 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl @@ -165,7 +165,7 @@ end q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew - dtnew = min(abs(dtnew), abs(tf - t - dt)) + dtnew = min(abs(dtnew), abs(tf - (t + dt))) @inbounds begin # Necessary for interpolation integ.k1 = k1 @@ -178,7 +178,7 @@ end integ.tprev = t integ.u = u - if (tf - t - dt) < convert(T, 1.0f-14) + if tf - (t + dt) < convert(T, 1.0f-14) integ.t = tf else if integ.tstops !== nothing && integ.tstops_idx <= length(integ.tstops) && diff --git a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl index dae609d2..a3fdecfd 100644 --- a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl @@ -134,7 +134,7 @@ end q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew - dtnew = min(abs(dtnew), abs(tf - t - dt)) + dtnew = min(abs(dtnew), abs(tf - (t + dt))) @inbounds begin # Necessary for interpolation integ.k1 = k1 @@ -152,7 +152,7 @@ end integ.tprev = t integ.u = u - if (tf - t - dt) < T(1.0e-14) + if tf - (t + dt) < T(1.0e-14) integ.t = tf else if integ.tstops !== nothing && integ.tstops_idx <= length(integ.tstops) && diff --git a/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl index 20594698..30545810 100644 --- a/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl @@ -167,7 +167,7 @@ end q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew - dtnew = min(abs(dtnew), abs(tf - t - dt)) + dtnew = min(abs(dtnew), abs(tf - (t + dt))) @inbounds begin # Necessary for interpolation integ.k1 = k1 @@ -188,7 +188,7 @@ end integ.tprev = t integ.u = u - if (tf - t - dt) < convert(T, 1.0f-14) + if tf - (t + dt) < convert(T, 1.0f-14) integ.t = tf else if integ.tstops !== nothing && integ.tstops_idx <= length(integ.tstops) && diff --git a/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl index 3a265fd0..1befac94 100644 --- a/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl @@ -276,7 +276,7 @@ end q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew - dtnew = min(abs(dtnew), abs(tf - t - dt)) + dtnew = min(abs(dtnew), abs(tf - (t + dt))) @inbounds begin # Necessary for interpolation integ.k1 = k1 @@ -297,7 +297,7 @@ end integ.tprev = t integ.u = u - if (tf - t - dt) < convert(T, 1.0e-14) + if tf - (t + dt) < convert(T, 1.0e-14) integ.t = tf else if integ.tstops !== nothing && integ.tstops_idx <= length(integ.tstops) && diff --git a/test/gpu_kernel_de/adaptive_endpoint.jl b/test/gpu_kernel_de/adaptive_endpoint.jl new file mode 100644 index 00000000..655e321d --- /dev/null +++ b/test/gpu_kernel_de/adaptive_endpoint.jl @@ -0,0 +1,32 @@ +using DiffEqGPU, SciMLBase, StaticArrays, Test + +# Drives `step!` directly with an iteration cap so that an integrator that stops +# advancing before `tf` fails the test instead of hanging the kernel loop. +function reaches_endpoint(alg, ::Type{T}, tf, λ, reltol) where {T} + f = ODEFunction{false}((u, p, t) -> SVector(-u[1], -p[1] * u[2])) + integ = DiffEqGPU.init( + alg, f, false, SVector{2, T}(1, 1), zero(T), T(tf), T(tf), SVector{1, T}(λ), + T(reltol / 1000), T(reltol), DiffEqGPU.DiffEqBase.ODE_DEFAULT_NORM, nothing, + CallbackSet(nothing), nothing + ) + ts = zeros(T, 2) + us = zeros(SVector{2, T}, 2) + for _ in 1:10_000 + integ.t < T(tf) || return integ.t == T(tf) + DiffEqGPU.step!(integ, ts, us) + end + return false +end + +@testset "Adaptive endpoint termination ($(nameof(typeof(alg))), $T)" for + alg in ( + GPUTsit5(), GPUTsit5IController(), GPUVern7(), GPUVern9(), + GPURosenbrock23(), GPURodas4(), GPURodas5P(), GPUKvaerno3(), GPUKvaerno5(), + ), + T in (Float32, Float64) + @test all( + reaches_endpoint(alg, T, tf, λ, reltol) + for tf in range(0.05, 2.0; length = 200), λ in (1, 10, 100), + reltol in (1.0e-3, 1.0e-4, 1.0e-6) + ) +end diff --git a/test/runtests.jl b/test/runtests.jl index 85eb326e..9cd579b6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -121,6 +121,9 @@ if GROUP == "CPU" @time @safetestset "Public generic interfaces" begin include("public_interface.jl") end + @time @safetestset "Adaptive endpoint termination" begin + include("gpu_kernel_de/adaptive_endpoint.jl") + end end if GROUP == "JLArrays" From 97c0e867f1c24f635fdf31a811fd1d21049f39fc Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Wed, 23 Sep 2026 04:42:26 -0400 Subject: [PATCH 2/4] Fuse adaptive GPUTsit5 stage arithmetic Apply `@muladd` to the adaptive Tsit5 stage combinations and embedded error estimate, matching the fixed-step stepper. The macro is scoped to that block so the step-size controller and endpoint logic keep their existing rounding, and the function keeps its `@inline` annotation: without it the mutable integrator is heap-allocated in device code (OpenCL `gpu_malloc` fails) and Enzyme on Julia 1.12 returns a wrong adaptive Float64 gradient. 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 --- AGENTS.md | 4 ++ .../perform_step/gpu_tsit5_perform_step.jl | 38 ++++++++++--------- 2 files changed, 25 insertions(+), 17 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index e0d82d61..c393617f 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -5,3 +5,7 @@ Preserve the existing `GPUTsit5` PI-controller behavior. Alternative controllers # Differentiation validation Gradient tests must use the package and its loaded extensions without adding differentiation rules in test code. Keep CPU Enzyme coverage alongside CUDA coverage. Exercise all active problem fields when validating transfer rules, including times and captured function values; parameter-only losses cannot detect discarded sensitivities. Keep initialization regressions enabled across backends and fix their owning dependencies rather than gating the tests. + +# Kernel stepper inlining + +Keep `@inline` on the `step!` methods in `src/ensemblegpukernel/perform_step/`. Without it the mutable integrator escapes into a call and is heap-allocated in device code: on OpenCL the resulting `gpu_malloc` fails and the kernel exits without writing output, and on Julia 1.12 Enzyme returns an incorrect adaptive Float64 gradient (the primal is unchanged). Scope `@muladd` to arithmetic blocks rather than removing the inline annotation. diff --git a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl index a3fdecfd..84106907 100644 --- a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl @@ -101,23 +101,27 @@ end while EEst > T(1.0) dt < T(1.0e-14) && error("dt Date: Sat, 26 Sep 2026 03:53:22 -0400 Subject: [PATCH 3/4] Keep the length-based endpoint snap in adaptive steppers The stall fix only needs the `dtnew` clamp `tf - (t + dt)`. Rewriting the snap test the same way made a step whose length already covers the interval skip the snap whenever `t + dt` rounds one ulp short of `tf`, so the solver took an extra ulp-sized step and shifted the final state by `eps(tf)`. Restore `(tf - t - dt) < 1e-14` in all eight steppers and add a state check for `u' = 1` over a span where `t0 + (tf - t0)` rounds below `tf`. 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 | 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 | 29 +++++++++++++++++++ 9 files changed, 37 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 f1fcdcd6..4d4e1eea 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 integ.tprev = t integ.u = u - if tf - (t + dt) < convert(T, 1.0f-14) + if (tf - t - dt) < convert(T, 1.0f-14) integ.t = tf else if integ.tstops !== nothing && integ.tstops_idx <= length(integ.tstops) && diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl index c70a941c..908b5a60 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 integ.tprev = t integ.u = u - if tf - (t + dt) < convert(T, 1.0f-14) + if (tf - t - dt) < convert(T, 1.0f-14) integ.t = tf else if integ.tstops !== nothing && integ.tstops_idx <= length(integ.tstops) && diff --git a/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl index c2d72cca..24af7449 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 integ.tprev = t integ.u = u - if tf - (t + dt) < convert(T, 1.0f-14) + if (tf - t - dt) < convert(T, 1.0f-14) integ.t = tf else if integ.tstops !== nothing && integ.tstops_idx <= length(integ.tstops) && diff --git a/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl index 01c871eb..0e267765 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 integ.tprev = t integ.u = u - if tf - (t + dt) < convert(T, 1.0f-14) + if (tf - t - dt) < convert(T, 1.0f-14) integ.t = tf else if integ.tstops !== nothing && integ.tstops_idx <= length(integ.tstops) && diff --git a/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl index df65cb82..04da2765 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 integ.tprev = t integ.u = u - if tf - (t + dt) < convert(T, 1.0f-14) + if (tf - t - dt) < convert(T, 1.0f-14) integ.t = tf else if integ.tstops !== nothing && integ.tstops_idx <= length(integ.tstops) && diff --git a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl index a3fdecfd..b211fa17 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 integ.tprev = t integ.u = u - if tf - (t + dt) < T(1.0e-14) + if (tf - t - dt) < T(1.0e-14) integ.t = tf else if integ.tstops !== nothing && integ.tstops_idx <= length(integ.tstops) && diff --git a/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl index 30545810..ac2c71c4 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 integ.tprev = t integ.u = u - if tf - (t + dt) < convert(T, 1.0f-14) + if (tf - t - dt) < convert(T, 1.0f-14) integ.t = tf else if integ.tstops !== nothing && integ.tstops_idx <= length(integ.tstops) && diff --git a/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl index 1befac94..d4f26910 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 integ.tprev = t integ.u = u - if tf - (t + dt) < convert(T, 1.0e-14) + if (tf - t - dt) < convert(T, 1.0e-14) integ.t = tf else if integ.tstops !== nothing && integ.tstops_idx <= length(integ.tstops) && diff --git a/test/gpu_kernel_de/adaptive_endpoint.jl b/test/gpu_kernel_de/adaptive_endpoint.jl index 655e321d..fc097c67 100644 --- a/test/gpu_kernel_de/adaptive_endpoint.jl +++ b/test/gpu_kernel_de/adaptive_endpoint.jl @@ -30,3 +30,32 @@ end reltol in (1.0e-3, 1.0e-4, 1.0e-6) ) end + +# `t0 + (tf - t0)` rounds one ulp below `tf` here. A step that already covers the +# interval must snap to `tf` rather than take a further ulp-sized step, which would +# leave `t == tf` but shift the state of this exactly integrable problem by `eps(tf)`. +@testset "Adaptive endpoint snap ($(nameof(typeof(alg))))" for + alg in ( + GPUTsit5(), GPUTsit5IController(), GPUVern7(), GPUVern9(), + GPURosenbrock23(), GPURodas4(), GPURodas5P(), GPUKvaerno3(), GPUKvaerno5(), + ) + t0, tf = 2.094426f7, 9.4428696f7 + @assert t0 + (tf - t0) != tf + f = ODEFunction{false}((u, p, t) -> SVector(1.0f0, 1.0f0)) + integ = DiffEqGPU.init( + alg, f, false, SVector(1.0f0, 1.0f0), t0, tf, tf - t0, SVector(0.0f0), + 1.0f-9, 1.0f-6, DiffEqGPU.DiffEqBase.ODE_DEFAULT_NORM, nothing, + CallbackSet(nothing), nothing + ) + ts, us = zeros(Float32, 2), zeros(SVector{2, Float32}, 2) + nsteps = 0 + while integ.t < tf && nsteps < 10 + DiffEqGPU.step!(integ, ts, us) + nsteps += 1 + end + @test integ.t == tf + @test nsteps == 1 + if alg isa Union{GPUTsit5, GPUTsit5IController, GPUVern7, GPUVern9, GPURosenbrock23, GPUKvaerno5} + @test integ.u == SVector(1.0f0, 1.0f0) .+ (tf - t0) + end +end From 6211289184fdaee28bc0667b7423b0f31a0510fb Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 26 Sep 2026 15:38:13 -0400 Subject: [PATCH 4/4] Clamp the next adaptive step from the accepted time Clamping `dtnew` against `tf - (t + dt)` before the tstops adjustment gave `dtnew = 0` whenever `t + dt` rounded to `tf` but the step then stopped at an earlier tstop, so the next step threw `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 | 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 | 64 ++++++++++++++++--- 9 files changed, 63 insertions(+), 25 deletions(-) diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl index 4d4e1eea..979cccf2 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl @@ -195,7 +195,6 @@ end q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew - dtnew = min(abs(dtnew), abs(tf - (t + dt))) @inbounds begin # Necessary for interpolation integ.k1 = k1 @@ -203,7 +202,6 @@ end end integ.dt = dt - integ.dtnew = dtnew integ.qold = qold integ.tprev = t integ.u = u @@ -223,6 +221,7 @@ end integ.t += dt end end + integ.dtnew = min(abs(dtnew), abs(tf - 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 908b5a60..03587507 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl @@ -263,7 +263,6 @@ end q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew - dtnew = min(abs(dtnew), abs(tf - (t + dt))) @inbounds begin # Necessary for interpolation integ.k1 = k1 @@ -271,7 +270,6 @@ end end integ.dt = dt - integ.dtnew = dtnew integ.qold = qold integ.tprev = t integ.u = u @@ -291,6 +289,7 @@ end integ.t += dt end end + integ.dtnew = min(abs(dtnew), abs(tf - 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 24af7449..7a9e253d 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl @@ -239,7 +239,6 @@ end q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew - dtnew = min(abs(dtnew), abs(tf - (t + dt))) @inbounds begin # Necessary for interpolation @unpack h21, h22, h23, h24, h25, h31, h32, h33, h34, h35 = integ.tab @@ -248,7 +247,6 @@ end end integ.dt = dt - integ.dtnew = dtnew integ.qold = qold integ.tprev = t integ.u = u @@ -268,6 +266,7 @@ end integ.t += dt end end + integ.dtnew = min(abs(dtnew), abs(tf - 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 0e267765..4ee64221 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl @@ -308,7 +308,6 @@ end q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew - dtnew = min(abs(dtnew), abs(tf - (t + dt))) @inbounds begin # Necessary for interpolation @unpack h21, h22, h23, h24, h25, h26, h27, h28, h31, h32, h33, h34, @@ -322,7 +321,6 @@ end h47 * k7 + h48 * k8 end integ.dt = dt - integ.dtnew = dtnew integ.qold = qold integ.tprev = t integ.u = u @@ -342,6 +340,7 @@ end integ.t += dt end end + integ.dtnew = min(abs(dtnew), abs(tf - 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 04da2765..1c946d6a 100644 --- a/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl @@ -165,7 +165,6 @@ end q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew - dtnew = min(abs(dtnew), abs(tf - (t + dt))) @inbounds begin # Necessary for interpolation integ.k1 = k1 @@ -173,7 +172,6 @@ end end integ.dt = dt - integ.dtnew = dtnew integ.qold = qold integ.tprev = t integ.u = u @@ -193,6 +191,7 @@ end integ.t += dt end end + integ.dtnew = min(abs(dtnew), abs(tf - 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 b211fa17..467b418c 100644 --- a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl @@ -134,7 +134,6 @@ end q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew - dtnew = min(abs(dtnew), abs(tf - (t + dt))) @inbounds begin # Necessary for interpolation integ.k1 = k1 @@ -147,7 +146,6 @@ end end integ.dt = dt - integ.dtnew = dtnew integ.qold = qold integ.tprev = t integ.u = u @@ -167,6 +165,7 @@ end integ.t += dt end end + integ.dtnew = min(abs(dtnew), abs(tf - 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 ac2c71c4..6e34f92f 100644 --- a/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl @@ -167,7 +167,6 @@ end q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew - dtnew = min(abs(dtnew), abs(tf - (t + dt))) @inbounds begin # Necessary for interpolation integ.k1 = k1 @@ -183,7 +182,6 @@ end end integ.dt = dt - integ.dtnew = dtnew integ.qold = qold integ.tprev = t integ.u = u @@ -203,6 +201,7 @@ end integ.t += dt end end + integ.dtnew = min(abs(dtnew), abs(tf - 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 d4f26910..96d0c67d 100644 --- a/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl @@ -276,7 +276,6 @@ end q = max(inv(qmax), min(inv(qmin), q / gamma)) qold = max(EEst, qoldinit) dtnew = dt / q #dtnew - dtnew = min(abs(dtnew), abs(tf - (t + dt))) @inbounds begin # Necessary for interpolation integ.k1 = k1 @@ -292,7 +291,6 @@ end end integ.dt = dt - integ.dtnew = dtnew integ.qold = qold integ.tprev = t integ.u = u @@ -312,6 +310,7 @@ end integ.t += dt end end + integ.dtnew = min(abs(dtnew), abs(tf - 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 fc097c67..f206f76e 100644 --- a/test/gpu_kernel_de/adaptive_endpoint.jl +++ b/test/gpu_kernel_de/adaptive_endpoint.jl @@ -1,4 +1,9 @@ -using DiffEqGPU, SciMLBase, StaticArrays, Test +using DiffEqGPU, KernelAbstractions, SciMLBase, StaticArrays, Test + +const ADAPTIVE_ALGS = ( + GPUTsit5(), GPUTsit5IController(), GPUVern7(), GPUVern9(), + GPURosenbrock23(), GPURodas4(), GPURodas5P(), GPUKvaerno3(), GPUKvaerno5(), +) # Drives `step!` directly with an iteration cap so that an integrator that stops # advancing before `tf` fails the test instead of hanging the kernel loop. @@ -19,10 +24,7 @@ function reaches_endpoint(alg, ::Type{T}, tf, λ, reltol) where {T} end @testset "Adaptive endpoint termination ($(nameof(typeof(alg))), $T)" for - alg in ( - GPUTsit5(), GPUTsit5IController(), GPUVern7(), GPUVern9(), - GPURosenbrock23(), GPURodas4(), GPURodas5P(), GPUKvaerno3(), GPUKvaerno5(), - ), + alg in ADAPTIVE_ALGS, T in (Float32, Float64) @test all( reaches_endpoint(alg, T, tf, λ, reltol) @@ -35,10 +37,7 @@ end # interval must snap to `tf` rather than take a further ulp-sized step, which would # leave `t == tf` but shift the state of this exactly integrable problem by `eps(tf)`. @testset "Adaptive endpoint snap ($(nameof(typeof(alg))))" for - alg in ( - GPUTsit5(), GPUTsit5IController(), GPUVern7(), GPUVern9(), - GPURosenbrock23(), GPURodas4(), GPURodas5P(), GPUKvaerno3(), GPUKvaerno5(), - ) + alg in ADAPTIVE_ALGS t0, tf = 2.094426f7, 9.4428696f7 @assert t0 + (tf - t0) != tf f = ODEFunction{false}((u, p, t) -> SVector(1.0f0, 1.0f0)) @@ -59,3 +58,50 @@ end @test integ.u == SVector(1.0f0, 1.0f0) .+ (tf - t0) end end + +# The first step covers the whole interval and `t0 + dt` rounds to `tf`, but the +# step must stop at the intermediate tstop, so the next step size has to be +# computed from the stop, not from `tf`. +@testset "Adaptive tstop before endpoint ($(nameof(typeof(alg))))" for alg in ADAPTIVE_ALGS + t0, tf, stop = 0.027657658f0, 1.0f0, 0.5f0 + dt0 = prevfloat(tf - t0) + @assert t0 + dt0 == tf + f = ODEFunction{false}((u, p, t) -> SVector(1.0f0, 1.0f0)) + integ = DiffEqGPU.init( + alg, f, false, SVector(1.0f0, 1.0f0), t0, tf, dt0, SVector(0.0f0), + 1.0f-9, 1.0f-6, DiffEqGPU.DiffEqBase.ODE_DEFAULT_NORM, Float32[stop], + CallbackSet(nothing), nothing + ) + ts, us = zeros(Float32, 2), zeros(SVector{2, Float32}, 2) + visited_stop = false + positive_dtnew = true + nsteps = 0 + while integ.t < tf && nsteps < 100 + DiffEqGPU.step!(integ, ts, us) + nsteps += 1 + visited_stop |= integ.t == stop + integ.t < tf && (positive_dtnew &= integ.dtnew > 0) + end + @test visited_stop + @test positive_dtnew + @test integ.t == tf +end + +@testset "Adaptive tstop before endpoint, public solve ($(nameof(typeof(alg))))" for + alg in ADAPTIVE_ALGS + t0, tf = 0.75f0, 1.0f0 + prob = ODEProblem{false}((u, p, t) -> SVector(1.0f0, 1.0f0), SVector(1.0f0, 1.0f0), (t0, tf)) + sol = solve( + EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0); + trajectories = 2, adaptive = true, dt = prevfloat(tf - t0), + tstops = Float32[0.875], abstol = 1.0f-9, reltol = 1.0f-6, save_everystep = false + ) + @test all(s -> s.t[end] == tf, sol.u) + # Stopping at the tstop goes through the dense-output interpolant, and the + # Float32 Verner interpolants carry O(1e-5) rounding error there even for u' = 1. + if alg isa Union{GPUVern7, GPUVern9} + @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 +end