Merge pull request #3070 from ChrisRackauckas-Claude/add-tsit5da-method

Add Tsit5DA: hybrid explicit/linear-implicit solver for DAEs

Author: ChrisRackauckas

lib/OrdinaryDiffEqRosenbrock/src/OrdinaryDiffEqRosenbrock.jl

@@ -23,7 +23,7 @@ import LinearSolve
 import LinearSolve: UniformScaling
 import ForwardDiff
 using FiniteDiff
-using LinearAlgebra: mul!, diag, diagm, I, Diagonal, norm, lu!
+using LinearAlgebra: mul!, diag, diagm, I, Diagonal, norm, lu, lu!
 using ADTypes
 import OrdinaryDiffEqCore, OrdinaryDiffEqDifferentiation
 
@@ -293,7 +293,7 @@ end
 
 export Rosenbrock23, Rosenbrock32, RosShamp4, Veldd4, Velds4, GRK4T, GRK4A,
     Ros4LStab, ROS3P, Rodas3, Rodas23W, Rodas3P, Rodas4, Rodas42, Rodas4P, Rodas4P2,
-    Rodas5, Rodas5P, Rodas5Pe, Rodas5Pr, Rodas6P,
+    Rodas5, Rodas5P, Rodas5Pe, Rodas5Pr, Rodas6P, HybridExplicitImplicitRK, Tsit5DA,
     RosenbrockW6S4OS, ROS34PW1a, ROS34PW1b, ROS34PW2, ROS34PW3, ROS34PRw, ROS3PRL,
     ROS3PRL2, ROK4a,
     ROS2, ROS2PR, ROS2S, ROS3, ROS3PR, Scholz4_7

lib/OrdinaryDiffEqRosenbrock/src/alg_utils.jl

@@ -34,6 +34,7 @@ alg_order(alg::Rodas5P) = 5
 alg_order(alg::Rodas5Pr) = 5
 alg_order(alg::Rodas5Pe) = 5
 alg_order(alg::Rodas6P) = 6
+alg_order(alg::HybridExplicitImplicitRK) = alg.order
 
 alg_adaptive_order(alg::Rosenbrock32) = 2
 alg_adaptive_order(alg::Rosenbrock23) = 3
@@ -61,12 +62,14 @@ isfsal(alg::Rodas42) = false
 isfsal(alg::Rodas4P) = false
 isfsal(alg::Rodas4P2) = false
 isfsal(alg::Rodas6P) = false
+isfsal(alg::HybridExplicitImplicitRK) = false
 
 function has_stiff_interpolation(
         ::Union{
             Rosenbrock23, Rosenbrock32, Rodas23W,
             Rodas3P, Rodas4, Rodas4P, Rodas4P2, Rodas5,
             Rodas5P, Rodas5Pe, Rodas5Pr, Rodas6P,
+            HybridExplicitImplicitRK,
         }
     )
     return true

lib/OrdinaryDiffEqRosenbrock/src/algorithms.jl

@@ -139,3 +139,68 @@ for (Alg, desc, refs, is_W) in [
         end
     end
 end
+
+################################################################################
+# HybridExplicitImplicitRK — generic tableau-based hybrid explicit/linear-implicit method
+################################################################################
+
+struct HybridExplicitImplicitRK{TabType, CS, AD, F, P, FDT, ST, CJ, StepLimiter, StageLimiter} <:
+    OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+    tab::TabType
+    order::Int
+    linsolve::F
+    precs::P
+    step_limiter!::StepLimiter
+    stage_limiter!::StageLimiter
+    autodiff::AD
+end
+
+function HybridExplicitImplicitRK(tab;
+        order,
+        chunk_size = Val{0}(), autodiff = AutoForwardDiff(),
+        standardtag = Val{true}(), concrete_jac = nothing,
+        diff_type = Val{:forward}(), linsolve = nothing,
+        precs = DEFAULT_PRECS, step_limiter! = trivial_limiter!,
+        stage_limiter! = trivial_limiter!
+    )
+    AD_choice, chunk_size, diff_type = _process_AD_choice(
+        autodiff, chunk_size, diff_type
+    )
+    return HybridExplicitImplicitRK{
+        typeof(tab), _unwrap_val(chunk_size), typeof(AD_choice), typeof(linsolve),
+        typeof(precs), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac), typeof(step_limiter!),
+        typeof(stage_limiter!),
+    }(
+        tab, order, linsolve, precs, step_limiter!,
+        stage_limiter!, AD_choice
+    )
+end
+
+# Keyword-only constructor for remake support
+function HybridExplicitImplicitRK(;
+        tab,
+        order,
+        chunk_size = Val{0}(), autodiff = AutoForwardDiff(),
+        standardtag = Val{true}(), concrete_jac = nothing,
+        diff_type = Val{:forward}(), linsolve = nothing,
+        precs = DEFAULT_PRECS, step_limiter! = trivial_limiter!,
+        stage_limiter! = trivial_limiter!
+    )
+    return HybridExplicitImplicitRK(tab;
+        order, chunk_size, autodiff, standardtag, concrete_jac,
+        diff_type, linsolve, precs, step_limiter!, stage_limiter!
+    )
+end
+
+"""
+A 12-stage order 5(4) hybrid explicit/linear-implicit method for semi-explicit index-1 DAEs.
+Differential variables are treated explicitly (like Tsit5), algebraic variables use Rosenbrock-type
+linear-implicit stages. Only the small algebraic Jacobian block needs factorization.
+For pure ODEs (no algebraic constraints), reduces to an explicit Runge-Kutta method.
+
+References:
+- Steinebach G., Rodas6P and Tsit5DA - two new Rosenbrock-type methods for DAEs.
+  arXiv:2511.21252, 2025.
+"""
+Tsit5DA(; kwargs...) = HybridExplicitImplicitRK(Tsit5DATableau; order = 5, kwargs...)

lib/OrdinaryDiffEqRosenbrock/src/interp_func.jl

@@ -21,6 +21,7 @@ function SciMLBase.interp_summary(
         Union{
             RosenbrockCombinedConstantCache, Rodas23WConstantCache, Rodas3PConstantCache,
             RosenbrockCache, Rodas23WCache, Rodas3PCache,
+            HybridExplicitImplicitConstantCache, HybridExplicitImplicitCache,
         },
     }
     return dense ? "specialized 3rd order \"free\" stiffness-aware interpolation" :

lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl

@@ -959,3 +959,159 @@ end
 ### RosenbrockW6S4O
 
 @RosenbrockW6S4OS(:cache)
+
+################################################################################
+
+### Tsit5DA - hybrid explicit/linear-implicit method for DAEs
+
+struct HybridExplicitImplicitConstantCache{TF, UF, Tab, JType, WType, F, AD} <: RosenbrockConstantCache
+    tf::TF
+    uf::UF
+    tab::Tab
+    J::JType
+    W::WType
+    linsolve::F
+    autodiff::AD
+    interp_order::Int
+end
+
+mutable struct HybridExplicitImplicitCache{
+        uType, rateType, uNoUnitsType, JType, WType, TabType,
+        TFType, UFType, F, JCType, GCType, RTolType, A,
+        StepLimiter, StageLimiter, DVType, AVType,
+        GZType, GYType, WZType, FZ,
+    } <: RosenbrockMutableCache
+    u::uType
+    uprev::uType
+    dense::Vector{rateType}
+    du::rateType
+    du1::rateType
+    du2::rateType
+    ks::Vector{rateType}
+    fsalfirst::rateType
+    fsallast::rateType
+    dT::rateType
+    J::JType
+    W::WType
+    tmp::rateType
+    atmp::uNoUnitsType
+    weight::uNoUnitsType
+    tab::TabType
+    tf::TFType
+    uf::UFType
+    linsolve_tmp::rateType
+    linsolve::F
+    jac_config::JCType
+    grad_config::GCType
+    reltol::RTolType
+    alg::A
+    step_limiter!::StepLimiter
+    stage_limiter!::StageLimiter
+    interp_order::Int
+    # DAE-specific fields
+    diff_vars::DVType
+    alg_vars::AVType
+    g_z::GZType         # n_g x n_g algebraic Jacobian block
+    g_y::GYType         # n_g x n_f coupling block
+    W_z::WZType         # -gamma * g_z (used for linear solve)
+    linsolve_tmp_z::FZ  # n_g-sized RHS for algebraic solve
+end
+function full_cache(c::HybridExplicitImplicitCache)
+    return [
+        c.u, c.uprev, c.dense..., c.du, c.du1, c.du2,
+        c.ks..., c.fsalfirst, c.fsallast, c.dT, c.tmp, c.atmp, c.weight, c.linsolve_tmp,
+    ]
+end
+
+function alg_cache(
+        alg::HybridExplicitImplicitRK, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+        dt, reltol, p, calck,
+        ::Val{false}, verbose
+    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    tf = TimeDerivativeWrapper(f, u, p)
+    uf = UDerivativeWrapper(f, t, p)
+    J, W = build_J_W(alg, u, uprev, p, t, dt, f, nothing, uEltypeNoUnits, Val(false))
+    tab = alg.tab(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
+    return HybridExplicitImplicitConstantCache(
+        tf, uf, tab, J, W, nothing, alg_autodiff(alg), size(tab.H, 1)
+    )
+end
+
+function alg_cache(
+        alg::HybridExplicitImplicitRK, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+        dt, reltol, p, calck,
+        ::Val{true}, verbose
+    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    tab = alg.tab(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
+    num_stages = size(tab.A, 1)
+    interp_order = size(tab.H, 1)
+
+    # Initialize vectors
+    dense = [zero(rate_prototype) for _ in 1:interp_order]
+    ks = [zero(rate_prototype) for _ in 1:num_stages]
+    du = zero(rate_prototype)
+    du1 = zero(rate_prototype)
+    du2 = zero(rate_prototype)
+    fsalfirst = zero(rate_prototype)
+    fsallast = zero(rate_prototype)
+    dT = zero(rate_prototype)
+    tmp = zero(rate_prototype)
+    atmp = similar(u, uEltypeNoUnits)
+    recursivefill!(atmp, false)
+    weight = similar(u, uEltypeNoUnits)
+    recursivefill!(weight, false)
+    linsolve_tmp = zero(rate_prototype)
+
+    tf = TimeGradientWrapper(f, uprev, p)
+    uf = UJacobianWrapper(f, t, p)
+
+    grad_config = build_grad_config(alg, f, tf, du1, t)
+    jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2)
+    J, W = build_J_W(alg, u, uprev, p, t, dt, f, jac_config, uEltypeNoUnits, Val(true))
+
+    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
+    Pl, Pr = wrapprecs(
+        alg.precs(
+            W, nothing, u, p, t, nothing, nothing, nothing,
+            nothing
+        )..., weight, tmp
+    )
+    linsolve = init(
+        linprob, alg.linsolve, alias = LinearAliasSpecifier(alias_A = true, alias_b = true),
+        Pl = Pl, Pr = Pr,
+        assumptions = LinearSolve.OperatorAssumptions(true),
+        verbose = verbose.linear_verbosity
+    )
+
+    # Detect algebraic variables from mass matrix
+    mass_matrix = f.mass_matrix
+    if mass_matrix === I
+        diff_vars = collect(1:length(u))
+        alg_vars = Int[]
+        g_z = zeros(eltype(u), 0, 0)
+        g_y = zeros(eltype(u), 0, 0)
+        W_z = zeros(eltype(u), 0, 0)
+        linsolve_tmp_z = zeros(eltype(u), 0)
+        linsolve_z = nothing
+    else
+        n = length(u)
+        diff_vars = findall(i -> mass_matrix[i, i] != 0, 1:n)
+        alg_vars = findall(i -> mass_matrix[i, i] == 0, 1:n)
+        n_g = length(alg_vars)
+        n_f = length(diff_vars)
+        g_z = zeros(eltype(u), n_g, n_g)
+        g_y = zeros(eltype(u), n_g, n_f)
+        W_z = zeros(eltype(u), n_g, n_g)
+        linsolve_tmp_z = zeros(eltype(u), n_g)
+    end
+
+    return HybridExplicitImplicitCache(
+        u, uprev, dense, du, du1, du2, ks,
+        fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf,
+        linsolve_tmp, linsolve, jac_config, grad_config, reltol, alg,
+        alg.step_limiter!, alg.stage_limiter!, interp_order,
+        diff_vars, alg_vars, g_z, g_y, W_z, linsolve_tmp_z
+    )
+end