Merge branch 'master' into ADTypesSwitch

.github/workflows/CI.yml

@@ -90,7 +90,7 @@ jobs:
         env:
           GROUP: ${{ matrix.group }}
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v4
+      - uses: codecov/codecov-action@v5
         with:
           token: ${{ secrets.CODECOV_TOKEN }}
           file: lcov.info

.github/workflows/Documentation.yml

@@ -32,7 +32,7 @@ jobs:
           DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key
         run: julia --project=docs/ --code-coverage=user docs/make.jl
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v4
+      - uses: codecov/codecov-action@v5
         with:
           token: ${{ secrets.CODECOV_TOKEN }}
           file: lcov.info

.github/workflows/Downstream.yml

@@ -30,7 +30,10 @@ jobs:
           - {user: SciML, repo: SciMLSensitivity.jl, group: Core3}
           - {user: SciML, repo: SciMLSensitivity.jl, group: Core4}
           - {user: SciML, repo: SciMLSensitivity.jl, group: Core5}
-          - {user: SciML, repo: ModelingToolkit.jl, group: All}
+          - {user: SciML, repo: ModelingToolkit.jl, group: InterfaceI}
+          - {user: SciML, repo: ModelingToolkit.jl, group: InterfaceII}
+          - {user: SciML, repo: ModelingToolkit.jl, group: Initialization}
+          - {user: SciML, repo: ModelingToolkit.jl, group: SymbolicIndexingInterface}
           - {user: SciML, repo: DiffEqDevTools.jl, group: Core}
           - {user: nathanaelbosch, repo: ProbNumDiffEq.jl, group: Downstream}
           - {user: SKopecz, repo: PositiveIntegrators.jl, group: Downstream}
@@ -73,7 +76,7 @@ jobs:
             exit(0)  # Exit immediately, as a success
           end
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v4
+      - uses: codecov/codecov-action@v5
         with:
           token: ${{ secrets.CODECOV_TOKEN }}
           file: lcov.info

lib/OrdinaryDiffEqCore/Project.toml

@@ -1,7 +1,8 @@
 name = "OrdinaryDiffEqCore"
 uuid = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 authors = ["ParamThakkar123 <[email protected]>"]
-version = "1.11.0"
+version = "1.13.0"
+
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -49,7 +50,7 @@ Accessors = "0.1.36"
 Adapt = "3.0, 4"
 ArrayInterface = "7"
 DataStructures = "0.18"
-DiffEqBase = "6.159"
+DiffEqBase = "6.160"
 DiffEqDevTools = "2.44.4"
 DocStringExtensions = "0.9"
 EnumX = "1"
@@ -70,7 +71,7 @@ Random = "<0.0.1, 1"
 RecursiveArrayTools = "2.36, 3"
 Reexport = "1.0"
 SafeTestsets = "0.1.0"
-SciMLBase = "2.59.2"
+SciMLBase = "2.62"
 SciMLOperators = "0.3"
 SciMLStructures = "1"
 SimpleUnPack = "1"

lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl

@@ -60,9 +60,9 @@ using DiffEqBase: check_error!, @def, _vec, _reshape
 
 using FastBroadcast: @.., True, False
 
-using SciMLBase: NoInit, CheckInit, _unwrap_val
+using SciMLBase: NoInit, CheckInit, OverrideInit, AbstractDEProblem, _unwrap_val
 
-import SciMLBase: alg_order
+import SciMLBase: AbstractNonlinearProblem, alg_order
 
 import DiffEqBase: calculate_residuals,
                    calculate_residuals!, unwrap_cache,
@@ -76,7 +76,8 @@ import Accessors: @reset
 
 using SciMLStructures: canonicalize, Tunable, isscimlstructure
 
-using SymbolicIndexingInterface: parameter_values, is_variable, variable_index, symbolic_type, NotSymbolic
+using SymbolicIndexingInterface: state_values, parameter_values, is_variable, variable_index,
+                                 symbolic_type, NotSymbolic
 
 const CompiledFloats = Union{Float32, Float64}
 import Preferences

lib/OrdinaryDiffEqCore/src/alg_utils.jl

@@ -17,6 +17,9 @@ end
 
 SciMLBase.forwarddiffs_model_time(alg::RosenbrockAlgorithm) = true
 
+SciMLBase.allows_late_binding_tstops(::OrdinaryDiffEqAlgorithm) = true
+SciMLBase.allows_late_binding_tstops(::DAEAlgorithm) = true
+
 # isadaptive is defined below.
 
 ## OrdinaryDiffEq Internal Traits

lib/OrdinaryDiffEqCore/src/caches/basic_caches.jl

@@ -73,7 +73,8 @@ function alg_cache(alg::CompositeAlgorithm{CS, Tuple{A1, A2, A3, A4, A5, A6}}, u
     args = (u, rate_prototype, uEltypeNoUnits,
         uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2, f, t, dt,
         reltol, p, calck, Val(V))
-    argT = map(typeof, args)
+    # Core.Typeof to turn uEltypeNoUnits into Type{uEltypeNoUnits} rather than DataType
+    argT = map(Core.Typeof, args)
     T1 = Base.promote_op(alg_cache, A1, argT...)
     T2 = Base.promote_op(alg_cache, A2, argT...)
     T3 = Base.promote_op(alg_cache, A3, argT...)

lib/OrdinaryDiffEqCore/src/dense/generic_dense.jl

@@ -168,10 +168,12 @@ function default_ode_interpolant(
         return ode_interpolant(Θ, integrator.dt, integrator.uprev,
             integrator.u, integrator.k, cache.cache5, idxs,
             deriv, integrator.differential_vars)
-    else # alg_choice == 6
+    elseif alg_choice == 6
         return ode_interpolant(Θ, integrator.dt, integrator.uprev,
             integrator.u, integrator.k, cache.cache6, idxs,
             deriv, integrator.differential_vars)
+    else
+        error("DefaultCache invalid alg_choice. File an issue.")
     end
 end
 
@@ -227,6 +229,8 @@ end
             ode_interpolant!(val, Θ, integrator.dt, integrator.uprev, integrator.u,
                 integrator.k, integrator.cache.cache6,
                 idxs, deriv, integrator.differential_vars)
+        else
+            error("DefaultCache invalid alg_choice. File an issue.")
         end
     else
         ode_interpolant!(val, Θ, integrator.dt, integrator.uprev, integrator.u,
@@ -256,10 +260,12 @@ function default_ode_interpolant!(
         return ode_interpolant!(val, Θ, integrator.dt, integrator.uprev,
             integrator.u, integrator.k, cache.cache5, idxs,
             deriv, integrator.differential_vars)
-    else # alg_choice == 6
+    elseif alg_choice == 6
         return ode_interpolant!(val, Θ, integrator.dt, integrator.uprev,
             integrator.u, integrator.k, cache.cache6, idxs,
             deriv, integrator.differential_vars)
+    else
+        error("DefaultCache invalid alg_choice. File an issue.")
     end
 end
 
@@ -380,6 +386,8 @@ function default_ode_extrapolant!(
         ode_interpolant!(val, Θ, integrator.t - integrator.tprev,
             integrator.uprev2, integrator.uprev,
             integrator.k, cache.cache6, idxs, deriv, integrator.differential_vars)
+    else
+        error("DefaultCache invalid alg_choice. File an issue.")
     end
 end
 
@@ -444,6 +452,8 @@ function default_ode_extrapolant(
         ode_interpolant(Θ, integrator.t - integrator.tprev,
             integrator.uprev2, integrator.uprev,
             integrator.k, cache.cache6, idxs, deriv, integrator.differential_vars)
+    else
+        error("DefaultCache invalid alg_choice. File an issue.")
     end
 end
 
@@ -810,6 +820,8 @@ function ode_interpolation(tval::Number, id::I, idxs, deriv::D, p,
                     cache.cache6) # update the kcurrent
                 val = ode_interpolant(Θ, dt, timeseries[i₋], timeseries[i₊], ks[i₊],
                     cache.cache6, idxs, deriv, differential_vars)
+            else
+                error("DefaultCache invalid alg_choice. File an issue.")
             end
         else
             _ode_addsteps!(ks[i₊], ts[i₋], timeseries[i₋], timeseries[i₊], dt, f, p,
@@ -892,6 +904,8 @@ function ode_interpolation!(out, tval::Number, id::I, idxs, deriv::D, p,
                     cache.cache6) # update the kcurrent
                 ode_interpolant!(out, Θ, dt, timeseries[i₋], timeseries[i₊], ks[i₊],
                     cache.cache6, idxs, deriv, differential_vars)
+            else
+                error("DefaultCache invalid alg_choice. File an issue.")
             end
         else
             _ode_addsteps!(ks[i₊], ts[i₋], timeseries[i₋], timeseries[i₊], dt, f, p,

lib/OrdinaryDiffEqCore/src/initialize_dae.jl

@@ -20,16 +20,6 @@ function BrownFullBasicInit(; abstol = 1e-10, nlsolve = nothing)
 end
 BrownFullBasicInit(abstol) = BrownFullBasicInit(; abstol = abstol, nlsolve = nothing)
 
-struct OverrideInit{T, F} <: DiffEqBase.DAEInitializationAlgorithm
-    abstol::T
-    nlsolve::F
-end
-
-function OverrideInit(; abstol = 1e-10, nlsolve = nothing)
-    OverrideInit(abstol, nlsolve)
-end
-OverrideInit(abstol) = OverrideInit(; abstol = abstol, nlsolve = nothing)
-
 ## Notes
 
 #=
@@ -111,7 +101,7 @@ default_nlsolve(alg, isinplace, u, initprob, autodiff = false) = alg
 
 ## If the initialization is trivial just use nothing alg
 function default_nlsolve(
-        ::Nothing, isinplace::Val{true}, u::Nothing, ::NonlinearProblem, autodiff = false)
+        ::Nothing, isinplace::Val{true}, u::Nothing, ::AbstractNonlinearProblem, autodiff = false)
     nothing
 end
 
@@ -121,7 +111,7 @@ function default_nlsolve(
 end
 
 function default_nlsolve(
-        ::Nothing, isinplace::Val{false}, u::Nothing, ::NonlinearProblem, autodiff = false)
+        ::Nothing, isinplace::Val{false}, u::Nothing, ::AbstractNonlinearProblem, autodiff = false)
     nothing
 end
 
@@ -132,7 +122,7 @@ function default_nlsolve(
 end
 
 function OrdinaryDiffEqCore.default_nlsolve(
-        ::Nothing, isinplace, u, ::NonlinearProblem, autodiff = false)
+        ::Nothing, isinplace, u, ::AbstractNonlinearProblem, autodiff = false)
     error("This ODE requires a DAE initialization and thus a nonlinear solve but no nonlinear solve has been loaded. To solve this problem, do `using OrdinaryDiffEqNonlinearSolve` or pass a custom `nlsolve` choice into the `initializealg`.")
 end
 
@@ -143,130 +133,52 @@ end
 
 ## NoInit
 
-function _initialize_dae!(integrator, prob::Union{ODEProblem, DAEProblem},
+function _initialize_dae!(integrator, prob::AbstractDEProblem,
         alg::NoInit, x::Union{Val{true}, Val{false}})
 end
 
 ## OverrideInit
 
-function _initialize_dae!(integrator, prob::Union{ODEProblem, DAEProblem},
+function _initialize_dae!(integrator, prob::AbstractDEProblem,
         alg::OverrideInit, isinplace::Union{Val{true}, Val{false}})
-    initializeprob = prob.f.initializeprob
-
-    if SciMLBase.has_update_initializeprob!(prob.f)
-        prob.f.update_initializeprob!(initializeprob, prob)
-    end
+    initializeprob = prob.f.initialization_data.initializeprob
 
     # If it doesn't have autodiff, assume it comes from symbolic system like ModelingToolkit
     # Since then it's the case of not a DAE but has initializeprob
     # In which case, it should be differentiable
-    isAD = if initializeprob.u0 === nothing
+    iu0 = state_values(initializeprob)
+    isAD = if iu0 === nothing
         AutoForwardDiff
     elseif has_autodiff(integrator.alg)
         alg_autodiff(integrator.alg) isa AutoForwardDiff
     else
         true
     end
 
-    alg = default_nlsolve(alg.nlsolve, isinplace, initializeprob.u0, initializeprob, isAD)
-    nlsol = solve(initializeprob, alg, abstol = integrator.opts.abstol, reltol = integrator.opts.reltol)
+    nlsolve_alg = default_nlsolve(alg.nlsolve, isinplace, iu0, initializeprob, isAD)
+
+    u0, p, success = SciMLBase.get_initial_values(prob, integrator, prob.f, alg, isinplace; nlsolve_alg, abstol = integrator.opts.abstol, reltol = integrator.opts.reltol)
+
     if isinplace === Val{true}()
-        integrator.u .= prob.f.initializeprobmap(nlsol)
+        integrator.u .= u0
     elseif isinplace === Val{false}()
-        integrator.u = prob.f.initializeprobmap(nlsol)
+        integrator.u = u0
     else
         error("Unreachable reached. Report this error.")
     end
-    if SciMLBase.has_initializeprobpmap(prob.f)
-        integrator.p = prob.f.initializeprobpmap(prob, nlsol)
-        sol = integrator.sol
-        @reset sol.prob.p = integrator.p
-        integrator.sol = sol
-    end
+    integrator.p = p
+    sol = integrator.sol
+    @reset sol.prob.p = integrator.p
+    integrator.sol = sol
 
-    if nlsol.retcode != ReturnCode.Success
+    if !success
         integrator.sol = SciMLBase.solution_new_retcode(integrator.sol,
             ReturnCode.InitialFailure)
     end
 end
 
 ## CheckInit
-struct CheckInitFailureError <: Exception
-    normresid::Any
-    abstol::Any
-end
-
-function Base.showerror(io::IO, e::CheckInitFailureError)
-    print(io,
-        "CheckInit specified but initialization not satisifed. normresid = $(e.normresid) > abstol = $(e.abstol)")
-end
-
-function _initialize_dae!(integrator, prob::ODEProblem, alg::CheckInit,
-        isinplace::Val{true})
-    @unpack p, t, f = integrator
-    M = integrator.f.mass_matrix
-    tmp = first(get_tmp_cache(integrator))
-    u0 = integrator.u
-
-    algebraic_vars = [all(iszero, x) for x in eachcol(M)]
-    algebraic_eqs = [all(iszero, x) for x in eachrow(M)]
-    (iszero(algebraic_vars) || iszero(algebraic_eqs)) && return
-    update_coefficients!(M, u0, p, t)
-    f(tmp, u0, p, t)
-    tmp .= ArrayInterface.restructure(tmp, algebraic_eqs .* _vec(tmp))
-
-    normresid = integrator.opts.internalnorm(tmp, t)
-    if normresid > integrator.opts.abstol
-        throw(CheckInitFailureError(normresid, integrator.opts.abstol))
-    end
-end
-
-function _initialize_dae!(integrator, prob::ODEProblem, alg::CheckInit,
-        isinplace::Val{false})
-    @unpack p, t, f = integrator
-    u0 = integrator.u
-    M = integrator.f.mass_matrix
-
-    algebraic_vars = [all(iszero, x) for x in eachcol(M)]
-    algebraic_eqs = [all(iszero, x) for x in eachrow(M)]
-    (iszero(algebraic_vars) || iszero(algebraic_eqs)) && return
-    update_coefficients!(M, u0, p, t)
-    du = f(u0, p, t)
-    resid = _vec(du)[algebraic_eqs]
-
-    normresid = integrator.opts.internalnorm(resid, t)
-    if normresid > integrator.opts.abstol
-        throw(CheckInitFailureError(normresid, integrator.opts.abstol))
-    end
-end
-
-function _initialize_dae!(integrator, prob::DAEProblem,
-        alg::CheckInit, isinplace::Val{true})
-    @unpack p, t, f = integrator
-    u0 = integrator.u
-    resid = get_tmp_cache(integrator)[2]
-
-    f(resid, integrator.du, u0, p, t)
-    normresid = integrator.opts.internalnorm(resid, t)
-    if normresid > integrator.opts.abstol
-        throw(CheckInitFailureError(normresid, integrator.opts.abstol))
-    end
-end
-
-function _initialize_dae!(integrator, prob::DAEProblem,
-        alg::CheckInit, isinplace::Val{false})
-    @unpack p, t, f = integrator
-    u0 = integrator.u
-
-    nlequation_oop = u -> begin
-        f((u - u0) / dt, u, p, t)
-    end
-
-    nlequation = (u, _) -> nlequation_oop(u)
-
-    resid = f(integrator.du, u0, p, t)
-    normresid = integrator.opts.internalnorm(resid, t)
-    if normresid > integrator.opts.abstol
-        throw(CheckInitFailureError(normresid, integrator.opts.abstol))
-    end
+function _initialize_dae!(integrator, prob::AbstractDEProblem, alg::CheckInit,
+        isinplace::Union{Val{true}, Val{false}})
+    SciMLBase.get_initial_values(prob, integrator, prob.f, alg, isinplace; abstol = integrator.opts.abstol)
 end