On clean master f1d1f81, GPUKvaerno3 and GPUKvaerno5 produce incorrect endpoints for the simple nonsingular mass-matrix equation 2I * du/dt = -u. Both fixed and adaptive stepping fail. GPURosenbrock23, GPURodas4, and GPURodas5P pass the same check.
Reproducer (CPU backend, no GPU required):
using DiffEqGPU, StaticArrays, LinearAlgebra, SciMLBase, KernelAbstractions, Test
f(u, p, t) = -u
jac(u, p, t) = -one(SMatrix{2,2,Float64})
@testset "Kvaerno mass matrix" begin
for alg in (GPUKvaerno3(), GPUKvaerno5()), adaptive in (false, true)
prob = ODEProblem(ODEFunction{false}(f; jac, mass_matrix=2.0I),
SVector(1.0, 2.0), (0.0, 1.0))
sol = solve(EnsembleProblem(prob; safetycopy=false), alg,
EnsembleGPUKernel(CPU(), 0.0); trajectories=2, adaptive,
dt=0.01, abstol=1e-8, reltol=1e-8, save_everystep=false)
expected = exp(-0.5) * prob.u0
@show alg adaptive sol.u[1].u[end]
@test sol.u[1].u[end] ≈ expected rtol=1e-5 atol=1e-8
end
endExpected first component: 0.6065306597126334. Observed:
| Algorithm |
adaptive |
First component |
| GPUKvaerno3 |
false |
0.5195192295470478 |
| GPUKvaerno3 |
true |
0.519751226885421 |
| GPUKvaerno5 |
false |
0.5664710897058612 |
| GPUKvaerno5 |
true |
0.5664883425771347 |
The reproduced test failure is Test Failed ... Expression: ≈((sol.u[1]).u[end], expected, rtol = 1.0e-5, atol = 1.0e-8). A larger clean-checkout test covering these four failures, the three Rosenbrock methods, and consistent stationary 0I problems reported 16 Pass, 4 Fail, 20 Total.
The explicit first stage and fixed-step FSAL cache in both Kvaerno step implementations use f(u,p,t) as the derivative, although the implicit stages solve the mass-matrix equation. As a diagnostic only, replacing those evaluations with f.mass_matrix \ f(u,p,t) made all four nonsingular cases pass (endpoint vector error norms between 1.2e-15 and 2.8e-9). This is not a complete DAE fix: directly inverting a singular mass matrix is invalid, and singular initialization/derivative handling needs a considered implementation. No patch is proposed here.
History investigation: 6101d69 introduced Kvaerno mass-matrix support in the Newton residual and Jacobian, leaving the explicit stage/FSAL evaluations unchanged. Before that commit these implicit stages ignored the mass matrix entirely. This identifies the incomplete support change by source history; it is not a completed passing-to-failing runtime bisect, and no passing historical nonidentity-mass Kvaerno baseline has been established.
Environment: DiffEqGPU 3.20.2 at the clean master commit above, SciMLBase 3.53.1, KernelAbstractions 0.9.42, StaticArrays 1.9.20. Clean-master reproduction and diagnostic patch ran in a separate worktree with a frozen manifest. GPU hardware and singular nonstationary DAEs were not tested.
Julia Version 1.12.7
Commit 6d172b025e4 (2026-08-15 08:05 UTC)
Build Info:
Official https://julialang.org release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 128 × AMD EPYC 7502 32-Core Processor
WORD_SIZE: 64
LLVM: libLLVM-18.1.7 (ORCJIT, znver2)
GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 128 virtual cores)
AI investigation: Codex (harness version unknown; model exact ID unknown). Local session ID: 01a07fcc-1c4f-7ee3-9a1e-51eaab7df293. No shareable conversation URL was exposed.
On clean master f1d1f81, GPUKvaerno3 and GPUKvaerno5 produce incorrect endpoints for the simple nonsingular mass-matrix equation
2I * du/dt = -u. Both fixed and adaptive stepping fail. GPURosenbrock23, GPURodas4, and GPURodas5P pass the same check.Reproducer (CPU backend, no GPU required):
Expected first component:
0.6065306597126334. Observed:The reproduced test failure is
Test Failed ... Expression: ≈((sol.u[1]).u[end], expected, rtol = 1.0e-5, atol = 1.0e-8). A larger clean-checkout test covering these four failures, the three Rosenbrock methods, and consistent stationary0Iproblems reported16 Pass, 4 Fail, 20 Total.The explicit first stage and fixed-step FSAL cache in both Kvaerno step implementations use
f(u,p,t)as the derivative, although the implicit stages solve the mass-matrix equation. As a diagnostic only, replacing those evaluations withf.mass_matrix \ f(u,p,t)made all four nonsingular cases pass (endpoint vector error norms between 1.2e-15 and 2.8e-9). This is not a complete DAE fix: directly inverting a singular mass matrix is invalid, and singular initialization/derivative handling needs a considered implementation. No patch is proposed here.History investigation: 6101d69 introduced Kvaerno mass-matrix support in the Newton residual and Jacobian, leaving the explicit stage/FSAL evaluations unchanged. Before that commit these implicit stages ignored the mass matrix entirely. This identifies the incomplete support change by source history; it is not a completed passing-to-failing runtime bisect, and no passing historical nonidentity-mass Kvaerno baseline has been established.
Environment: DiffEqGPU 3.20.2 at the clean master commit above, SciMLBase 3.53.1, KernelAbstractions 0.9.42, StaticArrays 1.9.20. Clean-master reproduction and diagnostic patch ran in a separate worktree with a frozen manifest. GPU hardware and singular nonstationary DAEs were not tested.
AI investigation: Codex (harness version unknown; model exact ID unknown). Local session ID: 01a07fcc-1c4f-7ee3-9a1e-51eaab7df293. No shareable conversation URL was exposed.