From a5eedd862f169b88d6241618c22ede1d75f9f659 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Wed, 30 Sep 2026 07:58:14 -0400 Subject: [PATCH 1/3] Fix GPUKvaerno3/5 explicit-stage derivatives for nonsingular mass matrices. ESDIRK first-stage and FSAL caches used f(u,p,t) as du/dt; solve M*du=f instead and reject singular mass matrices (DAEs) with a clear ArgumentError. Co-Authored-By: Chris Rackauckas Co-Authored-By: Cursor Agent Agent-Harness: Cursor Agent CLI 2026.09.28-64d2043 Agent-Model: auto Agent-Session: local session, transcript at /home/crackauc/sandbox/goals/issue-backlog/ibd/DiffEqGPU.jl/jobs/530-cursor/log.txt on amdci2.julia.csail.mit.edu --- src/ensemblegpukernel/gpukernel_algorithms.jl | 12 ++++++-- src/ensemblegpukernel/linalg/linsolve.jl | 28 +++++++++++++++++++ src/ensemblegpukernel/lowerlevel_solve.jl | 4 +++ .../perform_step/gpu_kvaerno3_perform_step.jl | 8 +++--- .../perform_step/gpu_kvaerno5_perform_step.jl | 8 +++--- 5 files changed, 50 insertions(+), 10 deletions(-) diff --git a/src/ensemblegpukernel/gpukernel_algorithms.jl b/src/ensemblegpukernel/gpukernel_algorithms.jl index 6f177384..7b013bc7 100644 --- a/src/ensemblegpukernel/gpukernel_algorithms.jl +++ b/src/ensemblegpukernel/gpukernel_algorithms.jl @@ -185,6 +185,9 @@ struct GPURodas5P{AD} <: GPUODEImplicitAlgorithm{AD} end Third-order ESDIRK method specialized for stiff `EnsembleGPUKernel` ODE solves. +Supports nonsingular mass matrices (`M * u' = f(u,p,t)`). Singular mass matrices +(DAEs) are rejected; use `GPURosenbrock23`, `GPURodas4`, or `GPURodas5P` instead. + # Keyword Arguments - `autodiff`: whether automatic differentiation is used for derivative generation. @@ -201,7 +204,8 @@ A `GPUKvaerno3` algorithm selector. # Throws -`solve` throws when the problem does not provide GPU-compatible derivative information. +`solve` throws when the problem does not provide GPU-compatible derivative information, +or when the mass matrix is singular. """ struct GPUKvaerno3{AD} <: GPUODEImplicitAlgorithm{AD} end @@ -210,6 +214,9 @@ struct GPUKvaerno3{AD} <: GPUODEImplicitAlgorithm{AD} end Fifth-order ESDIRK method specialized for stiff `EnsembleGPUKernel` ODE solves. +Supports nonsingular mass matrices (`M * u' = f(u,p,t)`). Singular mass matrices +(DAEs) are rejected; use `GPURosenbrock23`, `GPURodas4`, or `GPURodas5P` instead. + # Keyword Arguments - `autodiff`: whether automatic differentiation is used for derivative generation. @@ -226,7 +233,8 @@ A `GPUKvaerno5` algorithm selector. # Throws -`solve` throws when the problem does not provide GPU-compatible derivative information. +`solve` throws when the problem does not provide GPU-compatible derivative information, +or when the mass matrix is singular. """ struct GPUKvaerno5{AD} <: GPUODEImplicitAlgorithm{AD} end diff --git a/src/ensemblegpukernel/linalg/linsolve.jl b/src/ensemblegpukernel/linalg/linsolve.jl index 89b0bd64..f890a0ac 100644 --- a/src/ensemblegpukernel/linalg/linsolve.jl +++ b/src/ensemblegpukernel/linalg/linsolve.jl @@ -106,3 +106,31 @@ end LUp.U \ (LUp.L \ $(length(Sb) > 1 ? :(b[LUp.p, :]) : :(b[LUp.p]))) end end + +const _KVAERNO_SINGULAR_MASS_MSG = "GPUKvaerno3 and GPUKvaerno5 do not support DAEs (singular mass matrices). Use GPURosenbrock23, GPURodas4, or GPURodas5P instead." + +@inline function _check_kvaerno_mass_matrix(alg, mass_matrix) + if alg isa Union{GPUKvaerno3, GPUKvaerno5} && _is_singular_mass_matrix(mass_matrix) + throw(ArgumentError(_KVAERNO_SINGULAR_MASS_MSG)) + end + return nothing +end + +@inline _is_singular_mass_matrix(::LinearAlgebra.UniformScaling{Bool}) = false +@inline _is_singular_mass_matrix(M::LinearAlgebra.UniformScaling) = iszero(M.λ) +@inline _is_singular_mass_matrix(M) = iszero(det(M)) + +# Explicit first-stage / FSAL derivative for ESDIRK: solve M * du = f(u,p,t). +@inline function _kvaerno_explicit_du(f, u, p, t) + return _mass_matrix_solve_du(f.mass_matrix, f(u, p, t)) +end + +@inline _mass_matrix_solve_du(::LinearAlgebra.UniformScaling{Bool}, rhs) = rhs +@inline function _mass_matrix_solve_du(M::LinearAlgebra.UniformScaling, rhs) + iszero(M.λ) && throw(ArgumentError(_KVAERNO_SINGULAR_MASS_MSG)) + return rhs / M.λ +end +@inline function _mass_matrix_solve_du(M::StaticMatrix, rhs) + iszero(det(M)) && throw(ArgumentError(_KVAERNO_SINGULAR_MASS_MSG)) + return linear_solve(M, rhs) +end diff --git a/src/ensemblegpukernel/lowerlevel_solve.jl b/src/ensemblegpukernel/lowerlevel_solve.jl index b2812956..ad1e85db 100644 --- a/src/ensemblegpukernel/lowerlevel_solve.jl +++ b/src/ensemblegpukernel/lowerlevel_solve.jl @@ -88,6 +88,8 @@ function vectorized_solve( backend = get_backend(probs) backend = maybe_prefer_blocks(backend) + _check_kvaerno_mass_matrix(alg, prob.f.mass_matrix) + prob = convert(ImmutableODEProblem, prob) dt = convert(eltype(prob.tspan), dt) saveat_converted = nothing @@ -277,6 +279,8 @@ function vectorized_asolve( backend = get_backend(probs) backend = maybe_prefer_blocks(backend) + _check_kvaerno_mass_matrix(alg, prob.f.mass_matrix) + # Get the time type from the problem Tt = eltype(prob.tspan) diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl index ba4eca59..c7a2dd64 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl @@ -25,7 +25,7 @@ end if integ.u_modified - k1 = f(uprev, p, t) + k1 = _kvaerno_explicit_du(f, uprev, p, t) integ.u_modified = false else @inbounds k1 = integ.k1 @@ -79,7 +79,7 @@ k2 = z₄ ./ dt @inbounds begin # Necessary for interpolation - integ.k1 = f(integ.u, p, t) + integ.k1 = _kvaerno_explicit_du(f, integ.u, p, t) integ.k2 = k2 end @@ -113,7 +113,7 @@ end α32 = integ.tab if integ.u_modified - k1 = f(uprev, p, t) + k1 = _kvaerno_explicit_du(f, uprev, p, t) integ.u_modified = false else @inbounds k1 = integ.k1 @@ -134,7 +134,7 @@ end # FSAL Step 1 - k1 = f(uprev, p, t) + k1 = _kvaerno_explicit_du(f, uprev, p, t) z₁ = dt * k1 diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl index 56de80b6..9c885839 100644 --- a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl +++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl @@ -27,7 +27,7 @@ end if integ.u_modified - k1 = f(uprev, p, t) + k1 = _kvaerno_explicit_du(f, uprev, p, t) integ.u_modified = false else @inbounds k1 = integ.k1 @@ -110,7 +110,7 @@ k2 = z₇ ./ dt @inbounds begin # Necessary for interpolation - integ.k1 = f(integ.u, p, t) + integ.k1 = _kvaerno_explicit_du(f, integ.u, p, t) integ.k2 = k2 end @@ -146,7 +146,7 @@ end @unpack α31, α32, α41, α42, α43, α51, α52, α53, α61, α62, α63 = integ.tab if integ.u_modified - k1 = f(uprev, p, t) + k1 = _kvaerno_explicit_du(f, uprev, p, t) integ.u_modified = false else @inbounds k1 = integ.k1 @@ -167,7 +167,7 @@ end # FSAL Step 1 - k1 = f(uprev, p, t) + k1 = _kvaerno_explicit_du(f, uprev, p, t) z₁ = dt * k1 From db010e726f28d172c70ef70022cfdf49f4cddd77 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Wed, 30 Sep 2026 07:58:14 -0400 Subject: [PATCH 2/3] Add Kvaerno nonsingular mass-matrix regression and DAE rejection tests. Covers the issue #530 analytic check (M=2I) for fixed/adaptive GPUKvaerno3/5 and updates singular-mass DAE tests to expect ArgumentError. Co-Authored-By: Chris Rackauckas Co-Authored-By: Cursor Agent Agent-Harness: Cursor Agent CLI 2026.09.28-64d2043 Agent-Model: auto Agent-Session: local session, transcript at /home/crackauc/sandbox/goals/issue-backlog/ibd/DiffEqGPU.jl/jobs/530-cursor/log.txt on amdci2.julia.csail.mit.edu --- .../stiff_ode/gpu_ode_mass_matrix.jl | 32 +++++++++++++++++++ .../stiff_ode/gpu_ode_modelingtoolkit_dae.jl | 10 ++---- 2 files changed, 34 insertions(+), 8 deletions(-) diff --git a/test/gpu_kernel_de/stiff_ode/gpu_ode_mass_matrix.jl b/test/gpu_kernel_de/stiff_ode/gpu_ode_mass_matrix.jl index 3e089675..8935f327 100644 --- a/test/gpu_kernel_de/stiff_ode/gpu_ode_mass_matrix.jl +++ b/test/gpu_kernel_de/stiff_ode/gpu_ode_mass_matrix.jl @@ -67,3 +67,35 @@ sol = solve( expected = iszero(initial) ? initial : exp(-1.0f0 / mass.λ) @test all(s -> isapprox(s.u[end][1], expected; atol = 2.0f-6), sol.u) end + +# Analytic check for M * u' = -u with nonsingular M = 2I (issue #530). +# Expected endpoint is u0 * exp(-t/2); first component is exp(-1/2). +@testset "Kvaerno nonsingular mass matrix" begin + f(u, p, t) = -u + jac(u, p, t) = -one(SMatrix{2, 2, Float64}) + expected = exp(-0.5) * SVector(1.0, 2.0) + 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(backend, 0.0); trajectories = 2, adaptive, + dt = 0.01, abstol = 1.0e-8, reltol = 1.0e-8, save_everystep = false + ) + @test sol.u[1].u[end] ≈ expected rtol = 1.0e-5 atol = 1.0e-8 + end +end + +@testset "Kvaerno rejects singular mass matrix" begin + f = ODEFunction{false}((u, p, t) -> -u; mass_matrix = 0.0I) + prob = ODEProblem(f, SVector(0.0), (0.0, 1.0)) + for alg in (GPUKvaerno3(), GPUKvaerno5()) + @test_throws ArgumentError solve( + EnsembleProblem(prob; safetycopy = false), alg, + EnsembleGPUKernel(backend, 0.0); trajectories = 2, adaptive = false, + dt = 0.01, save_everystep = false + ) + end +end diff --git a/test/gpu_kernel_de/stiff_ode/gpu_ode_modelingtoolkit_dae.jl b/test/gpu_kernel_de/stiff_ode/gpu_ode_modelingtoolkit_dae.jl index 208df41e..09f1305c 100644 --- a/test/gpu_kernel_de/stiff_ode/gpu_ode_modelingtoolkit_dae.jl +++ b/test/gpu_kernel_de/stiff_ode/gpu_ode_modelingtoolkit_dae.jl @@ -92,27 +92,21 @@ end end @testset "GPUKvaerno3 DAE" begin - sol = solve( + @test_throws ArgumentError solve( monteprob, GPUKvaerno3(), EnsembleGPUKernel(backend), trajectories = 2, dt = 0.001f0, adaptive = false ) - @test length(sol.u) == 2 - @test !any(isnan, sol.u[1].u[end]) - @test abs(sol.u[1].u[end][1] + sol.u[1].u[end][2] - 1.0f0) < 0.01f0 end @testset "GPUKvaerno5 DAE" begin - sol = solve( + @test_throws ArgumentError solve( monteprob, GPUKvaerno5(), EnsembleGPUKernel(backend), trajectories = 2, dt = 0.001f0, adaptive = false ) - @test length(sol.u) == 2 - @test !any(isnan, sol.u[1].u[end]) - @test abs(sol.u[1].u[end][1] + sol.u[1].u[end][2] - 1.0f0) < 0.01f0 end # ============================================================================ From 04395938a4941d30ed472911c5903b6c4e624649 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Wed, 30 Sep 2026 11:45:30 -0400 Subject: [PATCH 3/3] Address review: Float32 mass-matrix tests, host-side DAE check, drop in-kernel guard. Rewrite Kvaerno mass-matrix tests in Float32 for Metal/oneAPI; reject singular M before device transfer; remove dead in-kernel det/throw; treat false*I as singular. Co-Authored-By: Chris Rackauckas Co-Authored-By: Cursor Agent Agent-Harness: Cursor Agent CLI 2026.09.28-64d2043 Agent-Model: auto Agent-Session: local session, transcript at /home/crackauc/sandbox/goals/issue-backlog/ibd/DiffEqGPU.jl/jobs/530-cursor/log.txt on amdci2.julia.csail.mit.edu --- src/ensemblegpukernel/linalg/linsolve.jl | 6 ++---- src/solve.jl | 3 +++ .../stiff_ode/gpu_ode_mass_matrix.jl | 20 +++++++++---------- 3 files changed, 15 insertions(+), 14 deletions(-) diff --git a/src/ensemblegpukernel/linalg/linsolve.jl b/src/ensemblegpukernel/linalg/linsolve.jl index f890a0ac..d8f934db 100644 --- a/src/ensemblegpukernel/linalg/linsolve.jl +++ b/src/ensemblegpukernel/linalg/linsolve.jl @@ -116,21 +116,19 @@ const _KVAERNO_SINGULAR_MASS_MSG = "GPUKvaerno3 and GPUKvaerno5 do not support D return nothing end -@inline _is_singular_mass_matrix(::LinearAlgebra.UniformScaling{Bool}) = false +# Includes UniformScaling{Bool}: `I` (λ=true) is nonsingular; `false*I` (λ=false) is singular. @inline _is_singular_mass_matrix(M::LinearAlgebra.UniformScaling) = iszero(M.λ) @inline _is_singular_mass_matrix(M) = iszero(det(M)) # Explicit first-stage / FSAL derivative for ESDIRK: solve M * du = f(u,p,t). +# Singular M is rejected on the host by `_check_kvaerno_mass_matrix` before launch. @inline function _kvaerno_explicit_du(f, u, p, t) return _mass_matrix_solve_du(f.mass_matrix, f(u, p, t)) end -@inline _mass_matrix_solve_du(::LinearAlgebra.UniformScaling{Bool}, rhs) = rhs @inline function _mass_matrix_solve_du(M::LinearAlgebra.UniformScaling, rhs) - iszero(M.λ) && throw(ArgumentError(_KVAERNO_SINGULAR_MASS_MSG)) return rhs / M.λ end @inline function _mass_matrix_solve_du(M::StaticMatrix, rhs) - iszero(det(M)) && throw(ArgumentError(_KVAERNO_SINGULAR_MASS_MSG)) return linear_solve(M, rhs) end diff --git a/src/solve.jl b/src/solve.jl index a700d788..9f15063e 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -406,6 +406,8 @@ function batch_solve_up_kernel( )... ) + _check_kvaerno_mass_matrix(alg, ensembleprob.prob.f.mass_matrix) + dev = ensemblealg.dev probs = _kernel_transfer(dev, adapted_probs) @@ -422,6 +424,7 @@ function batch_solve_up_kernel( kwargs..., callback = _callback ) end + solus = _kernel_transfer(CPU(), us) solts = _kernel_transfer(CPU(), ts) return (solts, solus) diff --git a/test/gpu_kernel_de/stiff_ode/gpu_ode_mass_matrix.jl b/test/gpu_kernel_de/stiff_ode/gpu_ode_mass_matrix.jl index 8935f327..c5ccdc99 100644 --- a/test/gpu_kernel_de/stiff_ode/gpu_ode_mass_matrix.jl +++ b/test/gpu_kernel_de/stiff_ode/gpu_ode_mass_matrix.jl @@ -69,33 +69,33 @@ sol = solve( end # Analytic check for M * u' = -u with nonsingular M = 2I (issue #530). -# Expected endpoint is u0 * exp(-t/2); first component is exp(-1/2). +# Expected endpoint is u0 * exp(-t/2). Float32 so Metal/oneAPI backends can run it. @testset "Kvaerno nonsingular mass matrix" begin f(u, p, t) = -u - jac(u, p, t) = -one(SMatrix{2, 2, Float64}) - expected = exp(-0.5) * SVector(1.0, 2.0) + jac(u, p, t) = -one(SMatrix{2, 2, Float32}) + expected = exp(-0.5f0) * SVector(1.0f0, 2.0f0) 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) + ODEFunction{false}(f; jac, mass_matrix = 2.0f0 * I), + SVector(1.0f0, 2.0f0), (0.0f0, 1.0f0) ) sol = solve( EnsembleProblem(prob; safetycopy = false), alg, EnsembleGPUKernel(backend, 0.0); trajectories = 2, adaptive, - dt = 0.01, abstol = 1.0e-8, reltol = 1.0e-8, save_everystep = false + dt = 0.01f0, abstol = 1.0f-8, reltol = 1.0f-8, save_everystep = false ) - @test sol.u[1].u[end] ≈ expected rtol = 1.0e-5 atol = 1.0e-8 + @test sol.u[1].u[end] ≈ expected rtol = 1.0f-5 atol = 1.0f-5 end end @testset "Kvaerno rejects singular mass matrix" begin - f = ODEFunction{false}((u, p, t) -> -u; mass_matrix = 0.0I) - prob = ODEProblem(f, SVector(0.0), (0.0, 1.0)) + f = ODEFunction{false}((u, p, t) -> -u; mass_matrix = 0.0f0 * I) + prob = ODEProblem(f, SVector(0.0f0), (0.0f0, 1.0f0)) for alg in (GPUKvaerno3(), GPUKvaerno5()) @test_throws ArgumentError solve( EnsembleProblem(prob; safetycopy = false), alg, EnsembleGPUKernel(backend, 0.0); trajectories = 2, adaptive = false, - dt = 0.01, save_everystep = false + dt = 0.01f0, save_everystep = false ) end end