Merge branch 'SciML:master' into master

.buildkite/pipeline.yml

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

.github/workflows/CI.yml

@@ -25,7 +25,6 @@ jobs:
           - ODEInterfaceRegression
           - Multithreading
         version:
-          - '1.9'
           - '1'
     steps:
       - uses: actions/checkout@v4

.github/workflows/Downstream.yml

@@ -14,7 +14,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        julia-version: [1.9, 1]
+        julia-version: [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.63.0"
+version = "6.69.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -23,7 +23,6 @@ LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
-LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
@@ -63,18 +62,17 @@ LineSearches = "7"
 LinearAlgebra = "1.9"
 LinearSolve = "2.1.10"
 Logging = "1.9"
-LoopVectorization = "0.12"
 MacroTools = "0.5"
 MuladdMacro = "0.2.1"
 NLsolve = "4"
-NonlinearSolve = "3"
+NonlinearSolve = "3.3"
 Polyester = "0.7"
-PreallocationTools = "0.4"
+PreallocationTools = "0.4.15"
 PrecompileTools = "1"
 Preferences = "1.3"
 RecursiveArrayTools = "2.36, 3"
 Reexport = "1.0"
-SciMLBase = "2"
+SciMLBase = "2.17"
 SciMLOperators = "0.3"
 SimpleNonlinearSolve = "1"
 SimpleUnPack = "1"
@@ -83,7 +81,7 @@ SparseDiffTools = "2.3"
 StaticArrayInterface = "1.2"
 StaticArrays = "1.0"
 TruncatedStacktraces = "1.2"
-julia = "1.9"
+julia = "1.10"
 
 [extras]
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"

README.md

@@ -15,6 +15,7 @@ ordinary differential equation solvers and utilities. While completely independe
 and usable on its own, users interested in using this
 functionality should check out [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl).
 
+
 ## Installation
 
 Assuming that you already have Julia correctly installed, it suffices to import

src/OrdinaryDiffEq.jl

@@ -18,8 +18,6 @@ using MuladdMacro, SparseArrays, FastClosures
 
 using LinearAlgebra
 
-using LoopVectorization
-
 import StaticArrayInterface
 
 import InteractiveUtils

src/dense/generic_dense.jl

@@ -360,6 +360,7 @@ function ode_interpolation(tvals, id::I, idxs, deriv::D, p,
             i₋ = max(1, _searchsortedlast(ts, t, i₋, tdir > 0))
             i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋
         end
+        id.sensitivitymode && error(SENSITIVITY_INTERP_MESSAGE)
         i₋₊ref[] = (i₋, i₊)
         dt = ts[i₊] - ts[i₋]
         Θ = iszero(dt) ? oneunit(t) / oneunit(dt) : (t - ts[i₋]) / dt
@@ -405,6 +406,7 @@ function ode_interpolation!(vals, tvals, id::I, idxs, deriv::D, p,
             i₋ = max(1, _searchsortedlast(ts, t, i₋, tdir > 0))
             i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋
         end
+        id.sensitivitymode && error(SENSITIVITY_INTERP_MESSAGE)
 
         dt = ts[i₊] - ts[i₋]
         Θ = iszero(dt) ? oneunit(t) / oneunit(dt) : (t - ts[i₋]) / dt
@@ -477,6 +479,7 @@ function ode_interpolation(tval::Number, id::I, idxs, deriv::D, p,
         i₋ = max(1, _searchsortedlast(ts, tval, 1, tdir > 0))
         i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋
     end
+    id.sensitivitymode && error(SENSITIVITY_INTERP_MESSAGE)
 
     @inbounds begin
         dt = ts[i₊] - ts[i₋]
@@ -525,6 +528,7 @@ function ode_interpolation!(out, tval::Number, id::I, idxs, deriv::D, p,
         i₋ = max(1, _searchsortedlast(ts, tval, 1, tdir > 0))
         i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋
     end
+    id.sensitivitymode && error(SENSITIVITY_INTERP_MESSAGE)
 
     @inbounds begin
         dt = ts[i₊] - ts[i₋]
@@ -687,8 +691,13 @@ Herimte Interpolation, chosen if no other dispatch for ode_interpolant
 @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{false}}, idxs::Nothing,
     T::Type{Val{0}}, differential_vars) # Default interpolant is Hermite
     #@.. broadcast=false (1-Θ)*y₀+Θ*y₁+Θ*(Θ-1)*((1-2Θ)*(y₁-y₀)+(Θ-1)*dt*k[1] + Θ*dt*k[2])
-    @inbounds (1 - Θ) * y₀ + Θ * y₁ +
-              differential_vars .* (Θ * (Θ - 1) * ((1 - 2Θ) * (y₁ - y₀) + (Θ - 1) * dt * k[1] + Θ * dt * k[2]))
+    if all(differential_vars)
+        @inbounds (1 - Θ) * y₀ + Θ * y₁ +
+                (Θ * (Θ - 1) * ((1 - 2Θ) * (y₁ - y₀) + (Θ - 1) * dt * k[1] + Θ * dt * k[2]))
+    else
+        @inbounds (1 - Θ) * y₀ + Θ * y₁ +
+                differential_vars .* (Θ * (Θ - 1) * ((1 - 2Θ) * (y₁ - y₀) + (Θ - 1) * dt * k[1] + Θ * dt * k[2]))
+    end
 end
 
 @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{true}}, idxs::Nothing,
@@ -755,10 +764,17 @@ Herimte Interpolation, chosen if no other dispatch for ode_interpolant
 @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{false}}, idxs::Nothing,
     T::Type{Val{1}}, differential_vars) # Default interpolant is Hermite
     #@.. broadcast=false k[1] + Θ*(-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(3*dt*k[1] + 3*dt*k[2] + 6*y₀ - 6*y₁) + 6*y₁)/dt
-    @inbounds (.!differential_vars).*(y₁ - y₀)/dt + differential_vars .*(
-              k[1] +
-              Θ * (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ +
-               Θ * (3 * dt * k[1] + 3 * dt * k[2] + 6 * y₀ - 6 * y₁) + 6 * y₁) / dt)
+    if all(differential_vars)
+        @inbounds  (
+            k[1] +
+            Θ * (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ +
+            Θ * (3 * dt * k[1] + 3 * dt * k[2] + 6 * y₀ - 6 * y₁) + 6 * y₁) / dt)
+    else
+        @inbounds (.!differential_vars).*(y₁ - y₀)/dt + differential_vars .*(
+                k[1] +
+                Θ * (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ +
+                Θ * (3 * dt * k[1] + 3 * dt * k[2] + 6 * y₀ - 6 * y₁) + 6 * y₁) / dt)
+    end
 end
 
 @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{true}}, idxs::Nothing,
@@ -826,8 +842,13 @@ Herimte Interpolation, chosen if no other dispatch for ode_interpolant
 @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{false}}, idxs::Nothing,
     T::Type{Val{2}}, differential_vars) # Default interpolant is Hermite
     #@.. broadcast=false (-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁) + 6*y₁)/(dt*dt)
-    @inbounds differential_vars .* (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ +
-               Θ * (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) + 6 * y₁) / (dt * dt)
+    if all(differential_vars)
+        @inbounds (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ +
+        Θ * (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) + 6 * y₁) / (dt * dt)
+    else
+        @inbounds differential_vars .* (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ +
+                Θ * (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) + 6 * y₁) / (dt * dt)
+    end
 end
 
 @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{true}}, idxs::Nothing,
@@ -887,7 +908,11 @@ Herimte Interpolation, chosen if no other dispatch for ode_interpolant
 @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{false}}, idxs::Nothing,
     T::Type{Val{3}}, differential_vars) # Default interpolant is Hermite
     #@.. broadcast=false (6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁)/(dt*dt*dt)
-    @inbounds differential_vars .* (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) / (dt * dt * dt)
+    if all(differential_vars)
+        @inbounds (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) / (dt * dt * dt)
+    else
+        @inbounds differential_vars .* (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) / (dt * dt * dt)
+    end
 end
 
 @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{true}}, idxs::Nothing,

src/dense/interpolants.jl

@@ -30,16 +30,28 @@ end
 
 @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k,
     cache::Union{FunctionMapConstantCache, FunctionMapCache},
-    idxs, T::Type{Val{0}}, differential_vars::Nothing)
+    idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing)
     y₀
 end
 
-@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k,
+@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k,
     cache::Union{FunctionMapConstantCache, FunctionMapCache},
     idxs, T::Type{Val{0}}, differential_vars::Nothing)
+    y₀[idxs]
+end
+
+@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k,
+    cache::Union{FunctionMapConstantCache, FunctionMapCache},
+    idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing)
     recursivecopy!(out, y₀)
 end
 
+@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k,
+    cache::Union{FunctionMapConstantCache, FunctionMapCache},
+    idxs, T::Type{Val{0}}, differential_vars::Nothing)
+    @views out[idxs] .= y₀[idxs]
+end
+
 """
 Hairer Norsett Wanner Solving Ordinary Differential Euations I - Nonstiff Problems Page 192
 """

src/dense/rosenbrock_interpolants.jl

@@ -314,4 +314,70 @@ end
                                     y₀[idxs] + y₁[idxs]) / dt
     out
 end
-#-
+
+# Second Derivative
+@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rosenbrock5ConstantCache,
+    idxs::Nothing, T::Type{Val{2}}, differential_vars)
+    @inbounds (-2 * k[1] + 2 * k[2] + Θ * (-6 * k[2] + 6 * k[3] - 12 * Θ * k[3])) / dt^2
+end
+
+@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rosenbrock5Cache, idxs::Nothing,
+    T::Type{Val{2}}, differential_vars)
+    @inbounds @.. broadcast=false (-2 * k[1] + 2 * k[2] +
+                                    Θ * (-6 * k[2] + 6 * k[3] - 12 * Θ * k[3]))/dt^2
+end
+
+@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k,
+    cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache},
+    idxs, T::Type{Val{2}}, differential_vars)
+    @.. broadcast=false (-2 * k[1][idxs] + 2 * k[2][idxs] +
+                          Θ * (-6 * k[2][idxs] + 6 * k[3][idxs] - 12 * Θ * k[3][idxs]))/dt^2
+end
+
+@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k,
+    cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache},
+    idxs::Nothing, T::Type{Val{2}}, differential_vars)
+    @.. broadcast=false out= (-2 * k[1] + 2 * k[2] +
+                              Θ * (-6 * k[2] + 6 * k[3] - 12 * Θ * k[3])) / dt^2
+    out
+end
+
+@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k,
+    cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache},
+    idxs, T::Type{Val{2}}, differential_vars)
+    @views @.. broadcast=false out= (-2 * k[1][idxs] + 2 * k[2][idxs] +
+                                     Θ *
+                                     (-6 * k[2][idxs] + 6 * k[3][idxs] - 12 * Θ * k[3][idxs])) / dt^2
+    out
+end
+
+# Third Derivative
+@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rosenbrock5ConstantCache,
+    idxs::Nothing, T::Type{Val{3}}, differential_vars)
+    @inbounds (-6 * k[2] + 6 * k[3] - 24 * Θ * k[3]) / dt^3
+end
+
+@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rosenbrock5Cache, idxs::Nothing,
+    T::Type{Val{3}}, differential_vars)
+    @inbounds @.. broadcast=false (-6 * k[2] + 6 * k[3] - 24 * Θ * k[3])/dt^3
+end
+
+@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k,
+    cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache},
+    idxs, T::Type{Val{3}}, differential_vars)
+    @.. broadcast=false (-6 * k[2][idxs] + 6 * k[3][idxs] - 24 * Θ * k[3][idxs])/dt^3
+end
+
+@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k,
+    cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache},
+    idxs::Nothing, T::Type{Val{3}}, differential_vars)
+    @.. broadcast=false out= (-6 * k[2] + 6 * k[3] - 24 * Θ * k[3]) / dt^3
+    out
+end
+
+@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k,
+    cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache},
+    idxs, T::Type{Val{3}}, differential_vars)
+    @views @.. broadcast=false out= (-6 * k[2][idxs] + 6 * k[3][idxs] - 24 * Θ * k[3][idxs]) / dt^3
+    out
+end

src/derivative_wrappers.jl

@@ -6,7 +6,7 @@ const FIRST_AUTODIFF_TGRAD_MESSAGE = """
                                1. Turn off automatic differentiation (e.g. Rosenbrock23() becomes
                                   Rosenbrock23(autodiff=false)). More details can be found at
                                   https://docs.sciml.ai/DiffEqDocs/stable/features/performance_overloads/
-                               2. Improving the compatibility of `f` with ForwardDiff.jl automatic 
+                               2. Improving the compatibility of `f` with ForwardDiff.jl automatic
                                   differentiation (using tools like PreallocationTools.jl). More details
                                   can be found at https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#Autodifferentiation-and-Dual-Numbers
                                3. Defining analytical Jacobians and time gradients. More details can be
@@ -48,7 +48,7 @@ const FIRST_AUTODIFF_JAC_MESSAGE = """
                                1. Turn off automatic differentiation (e.g. Rosenbrock23() becomes
                                   Rosenbrock23(autodiff=false)). More details can befound at
                                   https://docs.sciml.ai/DiffEqDocs/stable/features/performance_overloads/
-                               2. Improving the compatibility of `f` with ForwardDiff.jl automatic 
+                               2. Improving the compatibility of `f` with ForwardDiff.jl automatic
                                   differentiation (using tools like PreallocationTools.jl). More details
                                   can be found at https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#Autodifferentiation-and-Dual-Numbers
                                3. Defining analytical Jacobians. More details can be