Merge branch 'master' into fix-build-grad-config

.buildkite/pipeline.yml

@@ -20,7 +20,7 @@ steps:
     matrix:
       setup:
         version:
-          - "1.6"
+          - "1.9"
           - "1"
         group:
           - "Regression_I"

.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

.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}

Project.toml

@@ -1,7 +1,7 @@
 name = "OrdinaryDiffEq"
 uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 authors = ["Chris Rackauckas <[email protected]>", "Yingbo Ma <[email protected]>"]
-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"

src/alg_utils.jl

@@ -234,50 +234,18 @@ 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
     if !(alg_autodiff(alg) isa AutoForwardDiff) ||
        (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
 

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}

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,15 +644,15 @@ 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
 
 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]

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

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

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