From 16764b1d9c5f411be11bd37d611fa588d89af631 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Mon, 13 Nov 2023 17:35:56 +0100 Subject: [PATCH] Fix Rodas5 and Rodas5P --- src/dense/stiff_addsteps.jl | 6 ++++++ test/integrators/event_detection_tests.jl | 10 ++++++++++ 2 files changed, 16 insertions(+) diff --git a/src/dense/stiff_addsteps.jl b/src/dense/stiff_addsteps.jl index 87d1da4f35..2dfbe9ba2f 100644 --- a/src/dense/stiff_addsteps.jl +++ b/src/dense/stiff_addsteps.jl @@ -838,6 +838,9 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::Rosenbrock5Cache, veck8 = _vec(k8) @.. broadcast=false veck8=-vecu + # https://github.com/SciML/OrdinaryDiffEq.jl/issues/2055 + tmp = linsolve_tmp + @unpack h21, h22, h23, h24, h25, h26, h27, h28, h31, h32, h33, h34, h35, h36, h37, h38, h41, h42, h43, h44, h45, h46, h47, h48 = cache.tab @.. broadcast=false tmp=h21 * k1 + h22 * k2 + h23 * k3 + h24 * k4 + h25 * k5 + h26 * k6 + h27 * k7 + h28 * k8 @@ -1092,6 +1095,9 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::Rosenbrock5Cache{<:Arra @unpack h21, h22, h23, h24, h25, h26, h27, h28, h31, h32, h33, h34, h35, h36, h37, h38, h41, h42, h43, h44, h45, h46, h47, h48 = cache.tab + # https://github.com/SciML/OrdinaryDiffEq.jl/issues/2055 + tmp = linsolve_tmp + @inbounds @simd ivdep for i in eachindex(u) tmp[i] = h21 * k1[i] + h22 * k2[i] + h23 * k3[i] + h24 * k4[i] + h25 * k5[i] + h26 * k6[i] + h27 * k7[i] + h28 * k8[i] diff --git a/test/integrators/event_detection_tests.jl b/test/integrators/event_detection_tests.jl index 9caddeacda..307baf7759 100644 --- a/test/integrators/event_detection_tests.jl +++ b/test/integrators/event_detection_tests.jl @@ -1,5 +1,6 @@ using StaticArrays using OrdinaryDiffEq +using DiffEqDevTools using Test @inbounds @inline function ż(z, p, t) @@ -60,6 +61,15 @@ prob = ODEProblem(f, u0, tspan, p) sol = solve(prob, Tsit5(), callback = cb2) @test minimum(Array(sol)) > -40 +# https://github.com/SciML/OrdinaryDiffEq.jl/issues/2055 +for alg in (Rodas5(), Rodas5P()) + sol2 = solve(prob, alg; callback = cb2) + sol3 = appxtrue(sol, sol2) + @test sol3.errors[:L2] < 1e-5 + @test sol3.errors[:L∞] < 5e-5 + @test sol3.errors[:final] < 1e-5 +end + function fball(du, u, p, t) du[1] = u[2] du[2] = -p