Merge pull request #19 from Shreyas-Ekanathan/upstream

Upstream

.github/workflows/CI.yml

@@ -34,7 +34,6 @@ jobs:
           - OrdinaryDiffEqExponentialRK
           - OrdinaryDiffEqExtrapolation
           - OrdinaryDiffEqFIRK
-          - OrdinaryDiffEqFIRKGenerator
           - OrdinaryDiffEqFeagin
           - OrdinaryDiffEqFunctionMap
           - OrdinaryDiffEqHighOrderRK
@@ -91,8 +90,8 @@ 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
-          fail_ci_if_error: true
+          fail_ci_if_error: false

.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

@@ -73,7 +73,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,7 @@
 name = "OrdinaryDiffEqCore"
 uuid = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 authors = ["ParamThakkar123 <[email protected]>"]
-version = "1.11.0"
+version = "1.12.1"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -70,7 +70,7 @@ Random = "<0.0.1, 1"
 RecursiveArrayTools = "2.36, 3"
 Reexport = "1.0"
 SafeTestsets = "0.1.0"
-SciMLBase = "2.60"
+SciMLBase = "2.62"
 SciMLOperators = "0.3"
 SciMLStructures = "1"
 SimpleUnPack = "1"

lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl

@@ -60,7 +60,7 @@ 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
 

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
 
 #=
@@ -143,19 +133,15 @@ 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
@@ -168,105 +154,30 @@ function _initialize_dae!(integrator, prob::Union{ODEProblem, DAEProblem},
         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, initializeprob.u0, 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

lib/OrdinaryDiffEqDefault/test/default_solver_tests.jl

@@ -40,6 +40,8 @@ end
 prob_rober = ODEProblem(rober, [1.0, 0.0, 0.0], (0.0, 1e3), (0.04, 3e7, 1e4))
 sol = solve(prob_rober)
 rosensol = solve(prob_rober, AutoTsit5(Rosenbrock23(autodiff = false)))
+#test that cache is type stable
+@test typeof(sol.interp.cache.cache3) == typeof(rosensol.interp.cache.caches[2])
 # test that default has the same performance as AutoTsit5(Rosenbrock23()) (which we expect it to use for this).
 @test sol.stats.naccept == rosensol.stats.naccept
 @test sol.stats.nf == rosensol.stats.nf
@@ -50,6 +52,8 @@ rosensol = solve(prob_rober, AutoTsit5(Rosenbrock23(autodiff = false)))
 sol = solve(prob_rober, reltol = 1e-7, abstol = 1e-7)
 rosensol = solve(
     prob_rober, AutoVern7(Rodas5P(autodiff = false)), reltol = 1e-7, abstol = 1e-7)
+#test that cache is type stable
+@test typeof(sol.interp.cache.cache4) == typeof(rosensol.interp.cache.caches[2])
 # test that default has the same performance as AutoTsit5(Rosenbrock23()) (which we expect it to use for this).
 @test sol.stats.naccept == rosensol.stats.naccept
 @test sol.stats.nf == rosensol.stats.nf

lib/OrdinaryDiffEqFIRK/Project.toml

@@ -1,11 +1,12 @@
 name = "OrdinaryDiffEqFIRK"
 uuid = "5960d6e9-dd7a-4743-88e7-cf307b64f125"
 authors = ["ParamThakkar123 <[email protected]>"]
-version = "1.4.0"
+version = "1.5.0"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FastPower = "a4df4552-cc26-4903-aec0-212e50a0e84b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
@@ -15,12 +16,14 @@ OrdinaryDiffEqDifferentiation = "4302a76b-040a-498a-8c04-15b101fed76b"
 OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 
 [compat]
 DiffEqBase = "6.152.2"
 DiffEqDevTools = "2.44.4"
 FastBroadcast = "0.3.5"
+FastGaussQuadrature = "1.0.2"
 FastPower = "1"
 LinearAlgebra = "<0.0.1, 1"
 LinearSolve = "2.32.0"
@@ -33,16 +36,18 @@ Random = "<0.0.1, 1"
 RecursiveArrayTools = "3.27.0"
 Reexport = "1.2.2"
 SafeTestsets = "0.1.0"
+SciMLBase = "2.60.0"
 SciMLOperators = "0.3.9"
 Test = "<0.0.1, 1"
 julia = "1.10"
 
 [extras]
 DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e"
 ODEProblemLibrary = "fdc4e326-1af4-4b90-96e7-779fcce2daa5"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test", "ODEProblemLibrary"]
+test = ["DiffEqDevTools", "GenericSchur", "Random", "SafeTestsets", "Test", "ODEProblemLibrary"]