diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index e26d007fa4..cb1566aa93 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -20,7 +20,7 @@ steps: matrix: setup: version: - - "1.6" + - "1.9" - "1" group: - "Regression_I" diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index f558008b5d..1524481770 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -25,8 +25,8 @@ jobs: - ODEInterfaceRegression - Multithreading version: + - '1.9' - '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..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,1.6] + julia-version: [1.9, 1] os: [ubuntu-latest] package: - {user: SciML, repo: DelayDiffEq.jl, group: Interface} diff --git a/Project.toml b/Project.toml index 17f9afe04f..4a2c0f03a6 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.59.1" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" @@ -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.139" -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" 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..a30e544070 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 @@ -3008,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} diff --git a/src/dense/stiff_addsteps.jl b/src/dense/stiff_addsteps.jl index f61d0d2d75..2dfbe9ba2f 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 @@ -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 @@ -857,7 +860,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 @@ -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/src/solve.jl b/src/solve.jl index 9f6cd6228f..9d0036aa02 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -96,6 +96,16 @@ 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 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 || !isempty(saveat)) && verbose + @warn("Rosenbrock methods on equations without differential states do not bound the error on interpolations.") + 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 diff --git a/test/integrators/event_detection_tests.jl b/test/integrators/event_detection_tests.jl index 9caddeacda..97885568e0 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 (Rodas4(),Rodas4P(),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 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 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)