Merge branch 'SciML:master' into SharpFineRKN6

.buildkite/pipeline.yml

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

.github/workflows/CI.yml

@@ -25,10 +25,10 @@ jobs:
           - ODEInterfaceRegression
           - Multithreading
         version:
+          - '1.9'
           - '1'
-          - '1.6'
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}

.github/workflows/Documentation.yml

@@ -11,7 +11,7 @@ jobs:
   build:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@latest
         with:
           version: '1'

.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}
@@ -35,14 +35,14 @@ jobs:
           - {user: nathanaelbosch, repo: ProbNumDiffEq.jl, group: Downstream}
 
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.julia-version }}
           arch: x64
       - uses: julia-actions/julia-buildpkg@latest
       - name: Clone Downstream
-        uses: actions/checkout@v3
+        uses: actions/checkout@v4
         with:
           repository: ${{ matrix.package.user }}/${{ matrix.package.repo }}
           path: downstream

.github/workflows/FormatCheck.yml

@@ -21,7 +21,7 @@ jobs:
         with:
           version: ${{ matrix.julia-version }}
 
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - name: Install JuliaFormatter and format
         # This will use the latest version by default but you can set the version like so:
         #

.github/workflows/Invalidations.yml

@@ -19,12 +19,12 @@ jobs:
     - uses: julia-actions/setup-julia@v1
       with:
         version: '1'
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - uses: julia-actions/julia-buildpkg@v1
     - uses: julia-actions/julia-invalidations@v1
       id: invs_pr
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
       with:
         ref: ${{ github.event.repository.default_branch }}
     - uses: julia-actions/julia-buildpkg@v1

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.55.0"
+version = "6.59.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -45,42 +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.9"
 LineSearches = "7"
+LinearAlgebra = "1.9"
 LinearSolve = "2.1.10"
+Logging = "1.9"
 LoopVectorization = "0.12"
 MacroTools = "0.5"
 MuladdMacro = "0.2.1"
 NLsolve = "4.3"
-NonlinearSolve = "1.1"
-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"
+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.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/adams_utils.jl

@@ -7,7 +7,7 @@ function ϕ_and_ϕstar!(cache, du, k)
         ξ = dt = dts[1]
         ξ0 = zero(dt)
         β[1] = one(dt)
-        if typeof(cache) <: OrdinaryDiffEqMutableCache
+        if cache isa OrdinaryDiffEqMutableCache
             ϕ_n[1] .= du
             ϕstar_n[1] .= du
         else
@@ -18,7 +18,7 @@ function ϕ_and_ϕstar!(cache, du, k)
             ξ0 += dts[i]
             β[i] = β[i - 1] * ξ / ξ0
             ξ += dts[i]
-            if typeof(cache) <: OrdinaryDiffEqMutableCache
+            if cache isa OrdinaryDiffEqMutableCache
                 @.. broadcast=false ϕ_n[i]=ϕ_n[i - 1] - ϕstar_nm1[i - 1]
                 @.. broadcast=false ϕstar_n[i]=β[i] * ϕ_n[i]
             else
@@ -35,7 +35,7 @@ function ϕ_and_ϕstar!(cache::Union{VCABMConstantCache, VCABMCache}, du, k)
         ξ = dt = dts[1]
         ξ0 = zero(dt)
         β[1] = one(dt)
-        if typeof(cache) <: OrdinaryDiffEqMutableCache
+        if cache isa OrdinaryDiffEqMutableCache
             ϕ_n[1] .= du
             ϕstar_n[1] .= du
         else
@@ -46,7 +46,7 @@ function ϕ_and_ϕstar!(cache::Union{VCABMConstantCache, VCABMCache}, du, k)
             ξ0 += dts[i]
             β[i] = β[i - 1] * ξ / ξ0
             ξ += dts[i]
-            if typeof(cache) <: OrdinaryDiffEqMutableCache
+            if cache isa OrdinaryDiffEqMutableCache
                 @.. broadcast=false ϕ_n[i]=ϕ_n[i - 1] - ϕstar_nm1[i - 1]
                 @.. broadcast=false ϕstar_n[i]=β[i] * ϕ_n[i]
             else
@@ -63,7 +63,7 @@ function expand_ϕ_and_ϕstar!(cache, i)
     @unpack ξ, ξ0, β, dts, ϕstar_nm1, ϕ_n, ϕstar_n = cache
     ξ0 += dts[i]
     β[i] = β[i - 1] * ξ / ξ0
-    if typeof(cache) <: OrdinaryDiffEqMutableCache
+    if cache isa OrdinaryDiffEqMutableCache
         @.. broadcast=false ϕ_n[i]=ϕ_n[i - 1] - ϕstar_nm1[i - 1]
         @.. broadcast=false ϕstar_n[i]=β[i] * ϕ_n[i]
     else
@@ -77,13 +77,13 @@ function ϕ_np1!(cache, du_np1, k)
         @unpack ϕ_np1, ϕstar_n = cache
         for i in 1:k
             if i != 1
-                if typeof(cache) <: OrdinaryDiffEqMutableCache
+                if cache isa OrdinaryDiffEqMutableCache
                     @.. broadcast=false ϕ_np1[i]=ϕ_np1[i - 1] - ϕstar_n[i - 1]
                 else
                     ϕ_np1[i] = ϕ_np1[i - 1] - ϕstar_n[i - 1]
                 end
             else
-                if typeof(cache) <: OrdinaryDiffEqMutableCache
+                if cache isa OrdinaryDiffEqMutableCache
                     ϕ_np1[i] .= du_np1
                 else
                     ϕ_np1[i] = du_np1

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 &&
-                 !(typeof(prob.f.jac_prototype) <: 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
 
@@ -794,7 +753,7 @@ function default_controller(alg::Union{ExtrapolationMidpointDeuflhard,
 end
 
 function _digest_beta1_beta2(alg, cache, ::Val{QT}, _beta1, _beta2) where {QT}
-    if typeof(alg) <: OrdinaryDiffEqCompositeAlgorithm
+    if alg isa OrdinaryDiffEqCompositeAlgorithm
         beta2 = _beta2 === nothing ?
                 _composite_beta2_default(alg.algs, cache.current, Val(QT)) : _beta2
         beta1 = _beta1 === nothing ?
@@ -967,10 +926,10 @@ alg_can_repeat_jac(alg::OrdinaryDiffEqNewtonAdaptiveAlgorithm) = true
 alg_can_repeat_jac(alg::IRKC) = false
 
 function unwrap_alg(alg::SciMLBase.DEAlgorithm, is_stiff)
-    iscomp = typeof(alg) <: CompositeAlgorithm
+    iscomp = alg isa CompositeAlgorithm
     if !iscomp
         return alg
-    elseif typeof(alg.choice_function) <: AutoSwitchCache
+    elseif alg.choice_function isa AutoSwitchCache
         if is_stiff === nothing
             throwautoswitch(alg)
         end
@@ -987,10 +946,10 @@ end
 
 function unwrap_alg(integrator, is_stiff)
     alg = integrator.alg
-    iscomp = typeof(alg) <: CompositeAlgorithm
+    iscomp = alg isa CompositeAlgorithm
     if !iscomp
         return alg
-    elseif typeof(alg.choice_function) <: AutoSwitchCache
+    elseif alg.choice_function isa AutoSwitchCache
         if is_stiff === nothing
             throwautoswitch(alg)
         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
 

src/caches/extrapolation_caches.jl

@@ -1074,7 +1074,7 @@ function alg_cache(alg::ImplicitDeuflhardExtrapolation, u, rate_prototype,
     end
 
     #Update stage_number by the jacobian size
-    jac_dim = typeof(rate_prototype) <: Union{CompiledFloats, BigFloat} ? 1 :
+    jac_dim = rate_prototype isa Union{CompiledFloats, BigFloat} ? 1 :
               sum(size(rate_prototype))
     stage_number = stage_number .+ jac_dim
 

src/caches/kencarp_kvaerno_caches.jl

@@ -43,7 +43,7 @@ function alg_cache(alg::KenCarp3, u, rate_prototype, ::Type{uEltypeNoUnits},
         uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(true))
     fsalfirst = zero(rate_prototype)
 
-    if typeof(f) <: SplitFunction
+    if f isa SplitFunction
         k1 = zero(u)
         k2 = zero(u)
         k3 = zero(u)
@@ -228,7 +228,7 @@ function alg_cache(alg::KenCarp4, u, rate_prototype, ::Type{uEltypeNoUnits},
         uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(true))
     fsalfirst = zero(rate_prototype)
 
-    if typeof(f) <: SplitFunction
+    if f isa SplitFunction
         k1 = zero(u)
         k2 = zero(u)
         k3 = zero(u)
@@ -368,7 +368,7 @@ function alg_cache(alg::KenCarp5, u, rate_prototype, ::Type{uEltypeNoUnits},
         uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(true))
     fsalfirst = zero(rate_prototype)
 
-    if typeof(f) <: SplitFunction
+    if f isa SplitFunction
         k1 = zero(u)
         k2 = zero(u)
         k3 = zero(u)
@@ -455,7 +455,7 @@ function alg_cache(alg::KenCarp47, u, rate_prototype, ::Type{uEltypeNoUnits},
         uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(true))
     fsalfirst = zero(rate_prototype)
 
-    if typeof(f) <: SplitFunction
+    if f isa SplitFunction
         k1 = zero(u)
         k2 = zero(u)
         k3 = zero(u)
@@ -542,7 +542,7 @@ function alg_cache(alg::KenCarp58, u, rate_prototype, ::Type{uEltypeNoUnits},
         uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(true))
     fsalfirst = zero(rate_prototype)
 
-    if typeof(f) <: SplitFunction
+    if f isa SplitFunction
         k1 = zero(u)
         k2 = zero(u)
         k3 = zero(u)