Merge branch 'master' into master

.github/workflows/CI.yml

@@ -28,7 +28,7 @@ jobs:
           - '1'
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
       - uses: actions/cache@v4
@@ -46,6 +46,8 @@ jobs:
         env:
           GROUP: ${{ matrix.group }}
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v3
+      - uses: codecov/codecov-action@v4
         with:
+          token: ${{ secrets.CODECOV_TOKEN }}
           file: lcov.info
+          fail_ci_if_error: true

.github/workflows/Documentation.yml

@@ -23,6 +23,8 @@ 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@v3
+      - uses: codecov/codecov-action@v4
         with:
+          token: ${{ secrets.CODECOV_TOKEN }}
           file: lcov.info
+          fail_ci_if_error: true

.github/workflows/Downstream.yml

@@ -36,7 +36,7 @@ jobs:
 
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.julia-version }}
           arch: x64
@@ -64,6 +64,8 @@ jobs:
             exit(0)  # Exit immediately, as a success
           end
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v3
+      - uses: codecov/codecov-action@v4
         with:
+          token: ${{ secrets.CODECOV_TOKEN }}
           file: lcov.info
+          fail_ci_if_error: true

.github/workflows/Invalidations.yml

@@ -16,7 +16,7 @@ jobs:
     if: github.base_ref == github.event.repository.default_branch
     runs-on: ubuntu-latest
     steps:
-    - uses: julia-actions/setup-julia@v1
+    - uses: julia-actions/setup-julia@v2
       with:
         version: '1'
     - uses: actions/checkout@v4

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.72.0"
+version = "6.74.1"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"

src/OrdinaryDiffEq.jl

@@ -163,6 +163,7 @@ include("caches/extrapolation_caches.jl")
 include("caches/prk_caches.jl")
 include("caches/pdirk_caches.jl")
 include("caches/dae_caches.jl")
+include("caches/qprk_caches.jl")
 
 include("tableaus/low_order_rk_tableaus.jl")
 include("tableaus/high_order_rk_tableaus.jl")
@@ -174,6 +175,7 @@ include("tableaus/sdirk_tableaus.jl")
 include("tableaus/firk_tableaus.jl")
 include("tableaus/rkn_tableaus.jl")
 include("tableaus/rkc_tableaus.jl")
+include("tableaus/qprk_tableaus.jl")
 
 include("integrators/type.jl")
 include("integrators/controllers.jl")
@@ -209,6 +211,7 @@ include("perform_step/extrapolation_perform_step.jl")
 include("perform_step/prk_perform_step.jl")
 include("perform_step/pdirk_perform_step.jl")
 include("perform_step/dae_perform_step.jl")
+include("perform_step/qprk_perform_step.jl")
 
 include("dense/generic_dense.jl")
 include("dense/interpolants.jl")
@@ -351,7 +354,7 @@ export FunctionMap, Euler, Heun, Ralston, Midpoint, RK4, ExplicitRK, OwrenZen3,
        OwrenZen5,
        BS3, BS5, RK46NL, DP5, Tsit5, DP8, Vern6, Vern7, Vern8, TanYam7, TsitPap8,
        Vern9, Feagin10, Feagin12, Feagin14, CompositeAlgorithm, Anas5, RKO65, FRK65, PFRK87,
-       RKM, MSRK5, MSRK6, Stepanov5, SIR54
+       RKM, MSRK5, MSRK6, Stepanov5, SIR54, QPRK98
 
 export SSPRK22, SSPRK33, KYKSSPRK42, SSPRK53, SSPRK53_2N1, SSPRK53_2N2, SSPRK53_H, SSPRK63,
        SSPRK73, SSPRK83, SSPRK43, SSPRK432,
@@ -389,7 +392,7 @@ export Rosenbrock23, Rosenbrock32, RosShamp4, Veldd4, Velds4, GRK4T, GRK4A,
        Rodas5, Rodas5P,
        RosenbrockW6S4OS, ROS34PW1a, ROS34PW1b, ROS34PW2, ROS34PW3, ROS34PRw, ROS3PRL,
        ROS3PRL2,
-       ROS2PR, ROS2S, ROS3, ROS3PR, Scholz4_7
+       ROS2, ROS2PR, ROS2S, ROS3, ROS3PR, Scholz4_7
 
 export LawsonEuler, NorsettEuler, ETD1, ETDRK2, ETDRK3, ETDRK4, HochOst4, Exp4, EPIRK4s3A,
        EPIRK4s3B,

src/alg_utils.jl

@@ -284,6 +284,13 @@ function DiffEqBase.prepare_alg(alg::CompositeAlgorithm, u0, p, prob)
     CompositeAlgorithm(algs, alg.choice_function)
 end
 
+has_autodiff(alg::OrdinaryDiffEqAlgorithm) = false
+function has_autodiff(alg::Union{
+        OrdinaryDiffEqAdaptiveImplicitAlgorithm, OrdinaryDiffEqImplicitAlgorithm,
+        CompositeAlgorithm, OrdinaryDiffEqExponentialAlgorithm, DAEAlgorithm})
+    true
+end
+
 # Extract AD type parameter from algorithm, returning as Val to ensure type stability for boolean options.
 function _alg_autodiff(alg::OrdinaryDiffEqAlgorithm)
     error("This algorithm does not have an autodifferentiation option defined.")
@@ -626,6 +633,7 @@ alg_order(alg::Feagin12) = 12
 alg_order(alg::Feagin14) = 14
 alg_order(alg::PFRK87) = 8
 
+alg_order(alg::ROS2) = 2
 alg_order(alg::ROS2PR) = 2
 alg_order(alg::ROS2S) = 2
 alg_order(alg::ROS3) = 3
@@ -710,6 +718,8 @@ alg_order(alg::Alshina2) = 2
 alg_order(alg::Alshina3) = 3
 alg_order(alg::Alshina6) = 6
 
+alg_order(alg::QPRK98) = 9
+
 alg_maximum_order(alg) = alg_order(alg)
 alg_maximum_order(alg::CompositeAlgorithm) = maximum(alg_order(x) for x in alg.algs)
 alg_maximum_order(alg::ExtrapolationMidpointDeuflhard) = 2(alg.max_order + 1)

src/algorithms.jl

@@ -2900,6 +2900,11 @@ end
 - Shampine L.F. and Reichelt M., (1997) The MATLAB ODE Suite, SIAM Journal of
 Scientific Computing, 18 (1), pp. 1-22.
 
+#### ROS2
+
+- J. G. Verwer et al. (1999): A second-order Rosenbrock method applied to photochemical dispersion problems
+  https://doi.org/10.1137/S1064827597326651
+
 #### ROS3P
 
 - Lang, J. & Verwer, ROS3P—An Accurate Third-Order Rosenbrock Solver Designed for
@@ -2969,6 +2974,7 @@ University of Geneva, Switzerland.
 for Alg in [
     :Rosenbrock23,
     :Rosenbrock32,
+    :ROS2,
     :ROS2PR,
     :ROS2S,
     :ROS3,

src/algorithms/explicit_rk.jl

@@ -652,7 +652,23 @@ Base.@kwdef struct RKN4{StageLimiter, StepLimiter, Thread} <: OrdinaryDiffEqAlgo
     step_limiter!::StepLimiter = trivial_limiter!
     thread::Thread = False()
 end
-# for backwards compatibility
+                      
 function RKN4(stage_limiter!, step_limiter! = trivial_limiter!)
     RKN4(stage_limiter!, step_limiter!, False())
 end
+
+@doc explicit_rk_docstring(
+    "Runge–Kutta pairs of orders 9(8) for use in quadruple precision computations", "QPRK98",
+    references = "Kovalnogov VN, Fedorov RV, Karpukhina TV, Simos TE, Tsitouras C. Runge–Kutta pairs 
+    of orders 9 (8) for use in quadruple precision computations. Numerical Algorithms, 2023. 
+    doi: https://doi.org/10.1007/s11075-023-01632-8")
+Base.@kwdef struct QPRK98{StageLimiter, StepLimiter, Thread} <:
+                   OrdinaryDiffEqAdaptiveAlgorithm
+    stage_limiter!::StageLimiter = trivial_limiter!
+    step_limiter!::StepLimiter = trivial_limiter!
+    thread::Thread = False()
+end
+
+function QPRK98(stage_limiter!, step_limiter! = trivial_limiter!)
+    QPRK98(stage_limiter!, step_limiter!, False())
+end

src/caches/bdf_caches.jl

@@ -47,8 +47,11 @@ function alg_cache(alg::ABDF2, u, rate_prototype, ::Type{uEltypeNoUnits},
     fsalfirstprev = zero(rate_prototype)
     atmp = similar(u, uEltypeNoUnits)
     recursivefill!(atmp, false)
+    algebraic_vars = f.mass_matrix === I ? nothing :
+                     [all(iszero, x) for x in eachcol(f.mass_matrix)]
 
-    eulercache = ImplicitEulerCache(u, uprev, uprev2, fsalfirst, atmp, nlsolver)
+    eulercache = ImplicitEulerCache(
+        u, uprev, uprev2, fsalfirst, atmp, nlsolver, algebraic_vars)
 
     dtₙ₋₁ = one(dt)
     zₙ₋₁ = zero(u)

src/caches/qprk_caches.jl

@@ -0,0 +1,68 @@
+struct QPRK98ConstantCache <: OrdinaryDiffEqConstantCache end
+
+@cache struct QPRK98Cache{
+    uType, rateType, uNoUnitsType, StageLimiter, StepLimiter, Thread} <:
+              OrdinaryDiffEqMutableCache
+    u::uType
+    uprev::uType
+    fsalfirst::rateType
+    k2::rateType
+    k3::rateType
+    k4::rateType
+    k5::rateType
+    k6::rateType
+    k7::rateType
+    k8::rateType
+    k9::rateType
+    k10::rateType
+    k11::rateType
+    k12::rateType
+    k13::rateType
+    k14::rateType
+    k15::rateType
+    k16::rateType
+    utilde::uType
+    tmp::uType
+    atmp::uNoUnitsType
+    k::rateType
+    stage_limiter!::StageLimiter
+    step_limiter!::StepLimiter
+    thread::Thread
+end
+
+function alg_cache(alg::QPRK98, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
+        uprev, uprev2, f, t, dt, reltol, p, calck,
+        ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    k1 = zero(rate_prototype)
+    k2 = zero(rate_prototype)
+    k3 = zero(rate_prototype)
+    k4 = zero(rate_prototype)
+    k5 = zero(rate_prototype)
+    k6 = zero(rate_prototype)
+    k7 = zero(rate_prototype)
+    k8 = zero(rate_prototype)
+    k9 = zero(rate_prototype)
+    k10 = zero(rate_prototype)
+    k11 = zero(rate_prototype)
+    k12 = zero(rate_prototype)
+    k13 = zero(rate_prototype)
+    k14 = zero(rate_prototype)
+    k15 = zero(rate_prototype)
+    k16 = zero(rate_prototype)
+    utilde = zero(u)
+    tmp = zero(u)
+    atmp = similar(u, uEltypeNoUnits)
+    k = zero(rate_prototype)
+    recursivefill!(atmp, false)
+    QPRK98Cache(u, uprev, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, k12, k13, k14, k15,
+        k16, utilde, tmp, atmp, k, alg.stage_limiter!, alg.step_limiter!,
+        alg.thread)
+end
+
+function alg_cache(::QPRK98, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+        dt, reltol, p, calck,
+        ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    QPRK98ConstantCache()
+end

src/caches/rosenbrock_caches.jl

@@ -407,6 +407,12 @@ end
 
 ################################################################################
 
+### ROS2 methods
+
+@ROS2(:cache)
+
+################################################################################
+
 ### ROS23 methods
 
 @ROS23(:cache)

src/caches/sdirk_caches.jl

@@ -1,13 +1,14 @@
 abstract type SDIRKMutableCache <: OrdinaryDiffEqMutableCache end
 
-@cache mutable struct ImplicitEulerCache{uType, rateType, uNoUnitsType, N} <:
+@cache mutable struct ImplicitEulerCache{uType, rateType, uNoUnitsType, N, AV} <:
                       SDIRKMutableCache
     u::uType
     uprev::uType
     uprev2::uType
     fsalfirst::rateType
     atmp::uNoUnitsType
     nlsolver::N
+    algebraic_vars::AV
 end
 
 function alg_cache(alg::ImplicitEuler, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -22,7 +23,10 @@ function alg_cache(alg::ImplicitEuler, u, rate_prototype, ::Type{uEltypeNoUnits}
     atmp = similar(u, uEltypeNoUnits)
     recursivefill!(atmp, false)
 
-    ImplicitEulerCache(u, uprev, uprev2, fsalfirst, atmp, nlsolver)
+    algebraic_vars = f.mass_matrix === I ? nothing :
+                     [all(iszero, x) for x in eachcol(f.mass_matrix)]
+
+    ImplicitEulerCache(u, uprev, uprev2, fsalfirst, atmp, nlsolver, algebraic_vars)
 end
 
 mutable struct ImplicitEulerConstantCache{N} <: OrdinaryDiffEqConstantCache