From 4cc617d3eac4034843931fefec097ea217ad1933 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Thu, 2 Nov 2023 13:10:32 -0400 Subject: [PATCH 01/17] warn for Rosenbrock methods with no differential variables --- src/solve.jl | 9 +++++++++ test/integrators/check_error.jl | 8 ++++++++ 2 files changed, 17 insertions(+) diff --git a/src/solve.jl b/src/solve.jl index ba1c2264c2..c79c60d3e4 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -96,6 +96,15 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem, error("This solver is not able to use mass matrices.") end + if alg isa OrdinaryDiffEqRosenbrockAdaptiveAlgorithm && + !(prob.f.mass_matrix isa SciMLOperators.AbstractSciMLOperator) && + all(isequal(0), prob.f.mass_matrix) + alg isa Union{Rosenbrock23, Rosenbrock32} && error("Rosenbrock23 and Rosenbrock32 require at least one differential variable to produce valid solutions") + if dense && verbose + @warn("Rosenbrock methods have questionable interpolations when applied to equations without differential states.") + end + end + if !isempty(saveat) && dense @warn("Dense output is incompatible with saveat. Please use the SavingCallback from the Callback Library to mix the two behaviors.") end diff --git a/test/integrators/check_error.jl b/test/integrators/check_error.jl index 07abab0482..9287e35dad 100644 --- a/test/integrators/check_error.jl +++ b/test/integrators/check_error.jl @@ -40,3 +40,11 @@ for i in 1:(integrator.opts.maxiters) end end @test ok + +let + function f!(out, u, _, t) + out[1] = u[1] - sin(1e8*t) + end + mprob = ODEProblem(ODEFunction(f!, mass_matrix=[0.0;;]), [0.0], (0, 2.0)) + @test_throws ErrorException solve(mprob, Rosenbrock23()) +end From 3a76bb0853159d98e219f6740505310fabe2dd0a Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Thu, 2 Nov 2023 15:25:07 -0400 Subject: [PATCH 02/17] Update solve.jl --- src/solve.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/solve.jl b/src/solve.jl index c79c60d3e4..13932fb9a3 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -97,8 +97,9 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem, end if alg isa OrdinaryDiffEqRosenbrockAdaptiveAlgorithm && - !(prob.f.mass_matrix isa SciMLOperators.AbstractSciMLOperator) && + prob.f.mass_matrix isa AbstractMatrix && all(isequal(0), prob.f.mass_matrix) + # technically this should also warn for zero operators but those are hard to check for alg isa Union{Rosenbrock23, Rosenbrock32} && error("Rosenbrock23 and Rosenbrock32 require at least one differential variable to produce valid solutions") if dense && verbose @warn("Rosenbrock methods have questionable interpolations when applied to equations without differential states.") From a8c5f45a0d3b4d902a41a87d7f4f01fec0a44cbd Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 5 Nov 2023 11:28:33 -0500 Subject: [PATCH 03/17] Remove early defaulting and fix factorization algs No longer default to GMRES on SplitODEProblem when it's not an operator equation. LUFactorization works fine. GenericFactorization has stronger assumptions than is required. This removes the early defaulting since `linsolve=nothing` is type-inferrable for defaultalg since it's a single algorithm, and thus doing it early has no benefit but significant drawbacks. One drawback case is Radau since it picks the default based on real numbers but requires complex numbers. This is also the reason for the aforementioned factorization issue on SplitODEProblem, it was simply choosing wrong. --- src/alg_utils.jl | 47 +---------- src/algorithms.jl | 4 +- .../interface/linear_solver_split_ode_test.jl | 82 +++++-------------- 3 files changed, 24 insertions(+), 109 deletions(-) diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 2bd131675f..33f7f0b02b 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -234,33 +234,6 @@ function DiffEqBase.prepare_alg(alg::Union{ OrdinaryDiffEqExponentialAlgorithm{0, AD, FDT}}, u0::AbstractArray{T}, p, prob) where {AD, FDT, T} - if alg isa OrdinaryDiffEqExponentialAlgorithm - linsolve = nothing - elseif alg.linsolve === nothing - if (prob.f isa ODEFunction && prob.f.f isa AbstractSciMLOperator) - linsolve = LinearSolve.defaultalg(prob.f.f, u0) - elseif (prob.f isa SplitFunction && - prob.f.f1.f isa AbstractSciMLOperator) - linsolve = LinearSolve.defaultalg(prob.f.f1.f, u0) - if (linsolve === nothing) || (linsolve isa LinearSolve.DefaultLinearSolver && - linsolve.alg !== LinearSolve.DefaultAlgorithmChoice.KrylovJL_GMRES) - msg = "Split ODE problem do not work with factorization linear solvers. Bug detailed in https://github.com/SciML/OrdinaryDiffEq.jl/pull/1643. Defaulting to linsolve=KrylovJL()" - @warn msg - linsolve = KrylovJL_GMRES() - end - elseif (prob isa ODEProblem || prob isa DDEProblem) && - (prob.f.mass_matrix === nothing || - (prob.f.mass_matrix !== nothing && - !(prob.f.jac_prototype isa AbstractSciMLOperator))) - linsolve = LinearSolve.defaultalg(prob.f.jac_prototype, u0) - else - # If mm is a sparse matrix and A is a MatrixOperator, then let linear - # solver choose things later - linsolve = nothing - end - else - linsolve = alg.linsolve - end # If not using autodiff or norecompile mode or very large bitsize (like a dual number u0 already) # don't use a large chunksize as it will either error or not be beneficial @@ -268,16 +241,11 @@ function DiffEqBase.prepare_alg(alg::Union{ (isbitstype(T) && sizeof(T) > 24) || (prob.f isa ODEFunction && prob.f.f isa FunctionWrappersWrappers.FunctionWrappersWrapper) - if alg isa OrdinaryDiffEqExponentialAlgorithm - return remake(alg, chunk_size = Val{1}()) - else - return remake(alg, chunk_size = Val{1}(), linsolve = linsolve) - end + return remake(alg, chunk_size = Val{1}()) end L = StaticArrayInterface.known_length(typeof(u0)) if L === nothing # dynamic sized - # If chunksize is zero, pick chunksize right at the start of solve and # then do function barrier to infer the full solve x = if prob.f.colorvec === nothing @@ -287,19 +255,10 @@ function DiffEqBase.prepare_alg(alg::Union{ end cs = ForwardDiff.pickchunksize(x) - - if alg isa OrdinaryDiffEqExponentialAlgorithm - return remake(alg, chunk_size = Val{cs}()) - else - return remake(alg, chunk_size = Val{cs}(), linsolve = linsolve) - end + return remake(alg, chunk_size = Val{cs}()) else # statically sized cs = pick_static_chunksize(Val{L}()) - if alg isa OrdinaryDiffEqExponentialAlgorithm - return remake(alg, chunk_size = cs) - else - return remake(alg, chunk_size = cs, linsolve = linsolve) - end + return remake(alg, chunk_size = cs) end end diff --git a/src/algorithms.jl b/src/algorithms.jl index a8959e0606..4d8d0a7009 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -63,13 +63,11 @@ function DiffEqBase.remake(thing::Union{ ST, CJ}, OrdinaryDiffEqImplicitAlgorithm{CS, AD, FDT, ST, CJ }, - DAEAlgorithm{CS, AD, FDT, ST, CJ}}; - linsolve, kwargs...) where {CS, AD, FDT, ST, CJ} + DAEAlgorithm{CS, AD, FDT, ST, CJ}}; kwargs...) where {CS, AD, FDT, ST, CJ} T = SciMLBase.remaker_of(thing) T(; SciMLBase.struct_as_namedtuple(thing)..., chunk_size = Val{CS}(), autodiff = Val{AD}(), standardtag = Val{ST}(), concrete_jac = CJ === nothing ? CJ : Val{CJ}(), - linsolve = linsolve, kwargs...) end diff --git a/test/interface/linear_solver_split_ode_test.jl b/test/interface/linear_solver_split_ode_test.jl index a920ba3341..f4e7403446 100644 --- a/test/interface/linear_solver_split_ode_test.jl +++ b/test/interface/linear_solver_split_ode_test.jl @@ -5,15 +5,15 @@ using LinearAlgebra, LinearSolve import OrdinaryDiffEq.dolinsolve n = 8 -dt = 1 / 16 +dt = 1 / 1000 u0 = ones(n) tspan = (0.0, 1.0) M1 = 2ones(n) |> Diagonal #|> Array M2 = 2ones(n) |> Diagonal #|> Array -f1 = M1 |> MatrixOperator -f2 = M2 |> MatrixOperator +f1 = (du,u,p,t) -> du .= M1 * u +f2 = (du,u,p,t) -> du .= M2 * u prob = SplitODEProblem(f1, f2, u0, tspan) for algname in (:SBDF2, @@ -21,72 +21,32 @@ for algname in (:SBDF2, :KenCarp47) @testset "$algname" begin alg0 = @eval $algname() - alg1 = @eval $algname(linsolve = GenericFactorization()) + alg1 = @eval $algname(linsolve = LUFactorization()) kwargs = (dt = dt,) - # expected error message - msg = "Split ODE problem do not work with factorization linear solvers. Bug detailed in https://github.com/SciML/OrdinaryDiffEq.jl/pull/1643. Defaulting to linsolve=KrylovJL()" - @test_logs (:warn, msg) solve(prob, alg0; kwargs...) + solve(prob, alg0; kwargs...) @test DiffEqBase.__solve(prob, alg0; kwargs...).retcode == ReturnCode.Success - @test_broken DiffEqBase.__solve(prob, alg1; kwargs...).retcode == ReturnCode.Success + @test DiffEqBase.__solve(prob, alg1; kwargs...).retcode == ReturnCode.Success end end -##### -# deep dive -##### - -alg0 = KenCarp47() # passing case -alg1 = KenCarp47(linsolve = GenericFactorization()) # failing case - -## objects -ig0 = SciMLBase.init(prob, alg0; dt = dt) -ig1 = SciMLBase.init(prob, alg1; dt = dt) - -nl0 = ig0.cache.nlsolver -nl1 = ig1.cache.nlsolver - -lc0 = nl0.cache.linsolve -lc1 = nl1.cache.linsolve - -W0 = lc0.A -W1 = lc1.A - -# perform first step -OrdinaryDiffEq.loopheader!(ig0) -OrdinaryDiffEq.loopheader!(ig1) - -OrdinaryDiffEq.perform_step!(ig0, ig0.cache) -OrdinaryDiffEq.perform_step!(ig1, ig1.cache) - -@test !OrdinaryDiffEq.nlsolvefail(nl0) -@test OrdinaryDiffEq.nlsolvefail(nl1) - -# check operators -@test W0._concrete_form != W1._concrete_form -@test_broken W0._func_cache == W1._func_cache - -# check operator application -b = ones(n) -@test W0 * b == W1 * b -@test mul!(rand(n), W0, b) == mul!(rand(n), W1, b) -#@test W0 \ b == W1 \ b - -# check linear solve -lc0.b .= 1.0 -lc1.b .= 1.0 - -solve(lc0) -solve(lc1) +f1 = M1 |> MatrixOperator +f2 = M2 |> MatrixOperator +prob = SplitODEProblem(f1, f2, u0, tspan) -@test_broken lc0.u == lc1.u +for algname in (:SBDF2, + :SBDF3, + :KenCarp47) + @testset "$algname" begin + alg0 = @eval $algname() -# solve contried problem using OrdinaryDiffEq machinery -linres0 = dolinsolve(ig0, lc0; A = W0, b = b, linu = ones(n), reltol = 1e-8) -linres1 = dolinsolve(ig1, lc1; A = W1, b = b, linu = ones(n), reltol = 1e-8) + kwargs = (dt = dt,) -@test_broken linres0 == linres1 + solve(prob, alg0; kwargs...) + @test DiffEqBase.__solve(prob, alg0; kwargs...).retcode == ReturnCode.Success + end +end ### # custom linsolve function @@ -101,6 +61,4 @@ end alg = KenCarp47(linsolve = LinearSolveFunction(linsolve)) -@test solve(prob, alg).retcode == ReturnCode.Success - -nothing +@test solve(prob, alg).retcode == ReturnCode.Success \ No newline at end of file From f3988e16f4aa677cd6892df83768a0185a5bb25a Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 7 Nov 2023 13:41:40 +0100 Subject: [PATCH 04/17] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 50a2546430..791d164f04 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEq" uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" authors = ["Chris Rackauckas ", "Yingbo Ma "] -version = "6.58.1" +version = "6.58.2" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 28a9d44bf22966493e2760bb0fa583f1ee9ce91e Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Tue, 7 Nov 2023 15:05:39 -0500 Subject: [PATCH 05/17] fix test zero mass matrix (i.e. all algebraic equations breaks rosenbrock solvers) --- test/interface/mass_matrix_tests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/interface/mass_matrix_tests.jl b/test/interface/mass_matrix_tests.jl index 5b1a008b97..dc2604c658 100644 --- a/test/interface/mass_matrix_tests.jl +++ b/test/interface/mass_matrix_tests.jl @@ -300,6 +300,7 @@ end x0 = zeros(n, n) M = zeros(n * n) |> Diagonal |> Matrix +M[1,1] = true # zero mass matrix breaks rosenbrock f = ODEFunction(dynamics!, mass_matrix = M) tspan = (0, 10.0) prob = ODEProblem(f, x0, tspan) From 05a53fb54c566a862ba03f3ba0f3f1a54860d828 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Tue, 7 Nov 2023 18:12:03 -0500 Subject: [PATCH 06/17] Update src/solve.jl Co-authored-by: Christopher Rackauckas --- src/solve.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solve.jl b/src/solve.jl index 13932fb9a3..2fcf4f6895 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -102,7 +102,7 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem, # technically this should also warn for zero operators but those are hard to check for alg isa Union{Rosenbrock23, Rosenbrock32} && error("Rosenbrock23 and Rosenbrock32 require at least one differential variable to produce valid solutions") if dense && verbose - @warn("Rosenbrock methods have questionable interpolations when applied to equations without differential states.") + @warn("Rosenbrock methods on equations without differential states do not bound the error on interpolations.") end end From 52565de04d723ac50ce77398fa3ec49683f13873 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Tue, 7 Nov 2023 18:26:29 -0500 Subject: [PATCH 07/17] address review --- src/solve.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solve.jl b/src/solve.jl index 2fcf4f6895..901f3cf730 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -101,7 +101,7 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem, all(isequal(0), prob.f.mass_matrix) # technically this should also warn for zero operators but those are hard to check for alg isa Union{Rosenbrock23, Rosenbrock32} && error("Rosenbrock23 and Rosenbrock32 require at least one differential variable to produce valid solutions") - if dense && verbose + if (dense || !isempty(saveat)) && verbose @warn("Rosenbrock methods on equations without differential states do not bound the error on interpolations.") end end From db2f79a0b3c436384f749e18cbdc0b6e5d23cbcf Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 8 Nov 2023 07:55:10 +0100 Subject: [PATCH 08/17] bump lower to v1.9 --- .buildkite/pipeline.yml | 1 - .github/workflows/CI.yml | 1 - .github/workflows/Downstream.yml | 2 +- Project.toml | 34 ++++++++++++++++---------------- 4 files changed, 18 insertions(+), 20 deletions(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index e26d007fa4..8feaee42f1 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -20,7 +20,6 @@ steps: matrix: setup: version: - - "1.6" - "1" group: - "Regression_I" diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index f558008b5d..ab02143c1e 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -26,7 +26,6 @@ jobs: - Multithreading version: - '1' - - '1.6' steps: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v1 diff --git a/.github/workflows/Downstream.yml b/.github/workflows/Downstream.yml index 495a9c8407..ae25ad90da 100644 --- a/.github/workflows/Downstream.yml +++ b/.github/workflows/Downstream.yml @@ -14,7 +14,7 @@ jobs: strategy: fail-fast: false matrix: - julia-version: [1,1.6] + julia-version: [1] os: [ubuntu-latest] package: - {user: SciML, repo: DelayDiffEq.jl, group: Interface} diff --git a/Project.toml b/Project.toml index 791d164f04..50d4ee7af9 100644 --- a/Project.toml +++ b/Project.toml @@ -45,46 +45,46 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77" [compat] -ADTypes = "0.1, 0.2" -Adapt = "1.1, 2.0, 3.0" -ArrayInterface = "6, 7" +ADTypes = "0.2" +Adapt = "3.0" +ArrayInterface = "7" DataStructures = "0.18" DiffEqBase = "6.128.2" -DocStringExtensions = "0.8, 0.9" +DocStringExtensions = "0.9" ExponentialUtilities = "1.22" -FastBroadcast = "0.1.9, 0.2" +FastBroadcast = "0.2" FastClosures = "0.3" FiniteDiff = "2" ForwardDiff = "0.10.3" FunctionWrappersWrappers = "0.1" IfElse = "0.1" -InteractiveUtils = "1.6" +InteractiveUtils = "1.9" LineSearches = "7" -LinearAlgebra = "1.6" +LinearAlgebra = "1.9" LinearSolve = "2.1.10" -Logging = "1.6" +Logging = "1.9" LoopVectorization = "0.12" MacroTools = "0.5" MuladdMacro = "0.2.1" NLsolve = "4.3" -NonlinearSolve = "1.1, 2" -Polyester = "0.3, 0.4, 0.5, 0.6, 0.7" -PreallocationTools = "0.2, 0.3, 0.4" +NonlinearSolve = "2" +Polyester = "0.7" +PreallocationTools = "0.4" PrecompileTools = "1" Preferences = "1.3" RecursiveArrayTools = "2.36" -Reexport = "0.2, 1.0" -SciMLBase = "1.94, 2" +Reexport = "1.0" +SciMLBase = "2" SciMLNLSolve = "0.1" -SciMLOperators = "0.2.12, 0.3" +SciMLOperators = "0.3" SimpleNonlinearSolve = "0.1.4" SimpleUnPack = "1" -SparseArrays = "1.6" +SparseArrays = "1.9" SparseDiffTools = "2.3" StaticArrayInterface = "1.2" -StaticArrays = "0.11, 0.12, 1.0" +StaticArrays = "1.0" TruncatedStacktraces = "1.2" -julia = "1.6" +julia = "1.9" [extras] AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" From 5251909cf79dec6c3ae577ad08d8f532444d05a2 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 8 Nov 2023 11:24:57 +0100 Subject: [PATCH 09/17] Update .github/workflows/Downstream.yml Co-authored-by: Hendrik Ranocha --- .github/workflows/Downstream.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/Downstream.yml b/.github/workflows/Downstream.yml index ae25ad90da..61beadf612 100644 --- a/.github/workflows/Downstream.yml +++ b/.github/workflows/Downstream.yml @@ -14,7 +14,7 @@ jobs: strategy: fail-fast: false matrix: - julia-version: [1] + julia-version: [1.9, 1] os: [ubuntu-latest] package: - {user: SciML, repo: DelayDiffEq.jl, group: Interface} From 7988f638bcae2f853094c082b756e12b997d53c4 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 8 Nov 2023 11:25:02 +0100 Subject: [PATCH 10/17] Update .github/workflows/CI.yml Co-authored-by: Hendrik Ranocha --- .github/workflows/CI.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index ab02143c1e..1524481770 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -25,6 +25,7 @@ jobs: - ODEInterfaceRegression - Multithreading version: + - '1.9' - '1' steps: - uses: actions/checkout@v4 From a96b68a53cbca1087b1f4b28c318015fe50e67b4 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 8 Nov 2023 11:25:07 +0100 Subject: [PATCH 11/17] Update .buildkite/pipeline.yml Co-authored-by: Hendrik Ranocha --- .buildkite/pipeline.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 8feaee42f1..cb1566aa93 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -20,6 +20,7 @@ steps: matrix: setup: version: + - "1.9" - "1" group: - "Regression_I" From 52dad4fe099c752c94a3725d6ceabd1d1c854156 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 8 Nov 2023 11:26:02 +0100 Subject: [PATCH 12/17] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 50d4ee7af9..8132810a05 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEq" uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" authors = ["Chris Rackauckas ", "Yingbo Ma "] -version = "6.58.2" +version = "6.59.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 11a1a20846a4a1e44ad4ade521f35c5b64766a29 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 11 Nov 2023 16:56:36 -0500 Subject: [PATCH 13/17] Fix up Rodas5(p) addsteps! --- src/dense/stiff_addsteps.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/dense/stiff_addsteps.jl b/src/dense/stiff_addsteps.jl index f61d0d2d75..87d1da4f35 100644 --- a/src/dense/stiff_addsteps.jl +++ b/src/dense/stiff_addsteps.jl @@ -522,7 +522,7 @@ end function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::Rosenbrock5ConstantCache, always_calc_begin = false, allow_calc_end = true, force_calc_end = false) - if length(k) < 2 || always_calc_begin + if length(k) < 3 || always_calc_begin @unpack tf, uf = cache @unpack a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, C21, C31, C32, C41, C42, C43, C51, C52, C53, C54, C61, C62, C63, C64, C65, C71, C72, C73, C74, C75, C76, C81, C82, C83, C84, C85, C86, C87, gamma, d1, d2, d3, d4, d5, c2, c3, c4, c5 = cache.tab @@ -644,7 +644,7 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::Rosenbrock5ConstantCach h48 * k8 copyat_or_push!(k, 1, k₁) copyat_or_push!(k, 2, k₂) - copyat_or_push!(k, 2, k₃) + copyat_or_push!(k, 3, k₃) end nothing end @@ -652,7 +652,7 @@ end function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::Rosenbrock5Cache, always_calc_begin = false, allow_calc_end = true, force_calc_end = false) - if length(k) < 2 || always_calc_begin + if length(k) < 3 || always_calc_begin @unpack du, du1, du2, tmp, k1, k2, k3, k4, k5, k6, k7, k8, dT, J, W, uf, tf, linsolve_tmp, jac_config, fsalfirst, weight = cache @unpack a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, C21, C31, C32, C41, C42, C43, C51, C52, C53, C54, C61, C62, C63, C64, C65, C71, C72, C73, C74, C75, C76, C81, C82, C83, C84, C85, C86, C87, gamma, d1, d2, d3, d4, d5, c2, c3, c4, c5 = cache.tab @@ -857,7 +857,7 @@ end function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::Rosenbrock5Cache{<:Array}, always_calc_begin = false, allow_calc_end = true, force_calc_end = false) - if length(k) < 2 || always_calc_begin + if length(k) < 3 || always_calc_begin @unpack du, du1, du2, k1, k2, k3, k4, k5, k6, k7, k8, dT, J, W, uf, tf, linsolve_tmp, jac_config, fsalfirst = cache @unpack a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, C21, C31, C32, C41, C42, C43, C51, C52, C53, C54, C61, C62, C63, C64, C65, C71, C72, C73, C74, C75, C76, C81, C82, C83, C84, C85, C86, C87, gamma, d1, d2, d3, d4, d5, c2, c3, c4, c5 = cache.tab From 16764b1d9c5f411be11bd37d611fa588d89af631 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Mon, 13 Nov 2023 17:35:56 +0100 Subject: [PATCH 14/17] 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 From 0765f08bfd3ab9de5a27da5094e2e7c89abd14f1 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 14 Nov 2023 01:33:50 -0500 Subject: [PATCH 15/17] Update test/integrators/event_detection_tests.jl --- test/integrators/event_detection_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/integrators/event_detection_tests.jl b/test/integrators/event_detection_tests.jl index 307baf7759..97885568e0 100644 --- a/test/integrators/event_detection_tests.jl +++ b/test/integrators/event_detection_tests.jl @@ -62,7 +62,7 @@ sol = solve(prob, Tsit5(), callback = cb2) @test minimum(Array(sol)) > -40 # https://github.com/SciML/OrdinaryDiffEq.jl/issues/2055 -for alg in (Rodas5(), Rodas5P()) +for alg in (Rodas4(),Rodas4P(),Rodas5(), Rodas5P()) sol2 = solve(prob, alg; callback = cb2) sol3 = appxtrue(sol, sol2) @test sol3.errors[:L2] < 1e-5 From 66317a8d1d2f23570a89e27ca63e923714524599 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 14 Nov 2023 04:34:53 -0500 Subject: [PATCH 16/17] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 8132810a05..13af420651 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEq" uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" authors = ["Chris Rackauckas ", "Yingbo Ma "] -version = "6.59.0" +version = "6.59.1" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 13671ff86b92ac6318bd6135e8ab5cfc45dbf577 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Tue, 14 Nov 2023 19:56:24 +0200 Subject: [PATCH 17/17] Fix page number --- src/algorithms.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index 4d8d0a7009..a30e544070 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -3006,7 +3006,7 @@ end for Alg in [:LawsonEuler, :NorsettEuler, :ETDRK2, :ETDRK3, :ETDRK4, :HochOst4] """ Hochbruck, Marlis, and Alexander Ostermann. “Exponential Integrators.” Acta - Numerica 19 (2010): 209–86. doi:10.1017/S0962492910000048. + Numerica 19 (2010): 209–286. doi:10.1017/S0962492910000048. """ @eval struct $Alg{CS, AD, FDT, ST, CJ} <: OrdinaryDiffEqExponentialAlgorithm{CS, AD, FDT, ST, CJ}