Merge pull request #2092 from gstein3m/master

New methods Rodas23W / Rodas3P with error test for interpolation, see issue 2054 #2054
SciML · Jan 19, 2024 · 90a4bf6 · 90a4bf6
2 parents 88e374e + 526cefb
commit 90a4bf6
Show file tree

Hide file tree

Showing 11 changed files with 1,179 additions and 23 deletions.
diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl
@@ -385,7 +385,7 @@ export MagnusMidpoint, LinearExponential, MagnusLeapfrog, LieEuler, CayleyEuler,
     MagnusAdapt4, RKMK2, RKMK4, LieRK4, CG2, CG3, CG4a
 
 export Rosenbrock23, Rosenbrock32, RosShamp4, Veldd4, Velds4, GRK4T, GRK4A,
-    Ros4LStab, ROS3P, Rodas3, Rodas4, Rodas42, Rodas4P, Rodas4P2, Rodas5, Rodas5P,
+    Ros4LStab, ROS3P, Rodas3, Rodas23W, Rodas3P, Rodas4, Rodas42, Rodas4P, Rodas4P2, Rodas5, Rodas5P,
     RosenbrockW6S4OS, ROS34PW1a, ROS34PW1b, ROS34PW2, ROS34PW3
 
 export LawsonEuler, NorsettEuler, ETD1, ETDRK2, ETDRK3, ETDRK4, HochOst4, Exp4, EPIRK4s3A,

diff --git a/src/alg_utils.jl b/src/alg_utils.jl
@@ -30,6 +30,8 @@ isfsal(tab::DiffEqBase.ExplicitRKTableau) = tab.fsal
 
 # isfsal(alg::CompositeAlgorithm) = isfsal(alg.algs[alg.current])
 isfsal(alg::FunctionMap) = false
+isfsal(alg::Rodas3P) = false
+isfsal(alg::Rodas23W) = false
 isfsal(alg::Rodas5) = false
 isfsal(alg::Rodas5P) = false
 isfsal(alg::Rodas4) = false
@@ -620,9 +622,11 @@ alg_order(alg::Feagin14) = 14
 alg_order(alg::PFRK87) = 8
 
 alg_order(alg::Rosenbrock23) = 2
+alg_order(alg::Rodas23W) = 3
 alg_order(alg::Rosenbrock32) = 3
 alg_order(alg::ROS3P) = 3
 alg_order(alg::Rodas3) = 3
+alg_order(alg::Rodas3P) = 3
 alg_order(alg::ROS34PW1a) = 3
 alg_order(alg::ROS34PW1b) = 3
 alg_order(alg::ROS34PW2) = 3

diff --git a/src/algorithms.jl b/src/algorithms.jl
@@ -2910,6 +2910,11 @@ Scientific Computing, 18 (1), pp. 1-22.
 - Kaps, P. & Rentrop, Generalized Runge-Kutta methods of order four with stepsize control
   for stiff ordinary differential equations. P. Numer. Math. (1979) 33: 55. doi:10.1007/BF01396495
 
+#### Rodas23W, Rodas3P
+
+- Steinebach G., Rodas23W / Rodas32P - a Rosenbrock-type method for DAEs with additional error estimate for dense output and Julia implementation,
+ in progress
+
 #### Rodas4P
 
 - Steinebach G. Order-reduction of ROW-methods for DAEs and method of lines
@@ -2921,10 +2926,14 @@ Scientific Computing, 18 (1), pp. 1-22.
   Differential-Algebraic Equations Forum. Springer, Cham. https://doi.org/10.1007/978-3-030-53905-4_6
 
 #### Rodas5
-
 - Di Marzo G. RODAS5(4) – Méthodes de Rosenbrock d’ordre 5(4) adaptées aux problemes
 différentiels-algébriques. MSc mathematics thesis, Faculty of Science,
 University of Geneva, Switzerland.
+
+#### Rodas5P
+- Steinebach G.   Construction of Rosenbrock–Wanner method Rodas5P and numerical benchmarks within the Julia Differential Equations package. 
+ In: BIT Numerical Mathematics, 63(2), 2023
+
 =#
 
 for Alg in [
@@ -2942,6 +2951,8 @@ for Alg in [
     :GRK4T,
     :GRK4A,
     :Ros4LStab,
+    :Rodas23W,
+    :Rodas3P,
     :Rodas4,
     :Rodas42,
     :Rodas4P,

diff --git a/src/caches/rosenbrock_caches.jl b/src/caches/rosenbrock_caches.jl
@@ -418,6 +418,215 @@ jac_cache(c::Rosenbrock4Cache) = (c.J, c.W)
 
 ### Rodas methods
 
+struct Rodas23WConstantCache{TF, UF, Tab, JType, WType, F, AD} <: OrdinaryDiffEqConstantCache
+    tf::TF
+    uf::UF
+    tab::Tab
+    J::JType
+    W::WType
+    linsolve::F
+    autodiff::AD
+end
+
+struct Rodas3PConstantCache{TF, UF, Tab, JType, WType, F, AD} <: OrdinaryDiffEqConstantCache
+    tf::TF
+    uf::UF
+    tab::Tab
+    J::JType
+    W::WType
+    linsolve::F
+    autodiff::AD
+end
+
+@cache mutable struct Rodas23WCache{uType, rateType, uNoUnitsType, JType, WType, TabType,
+    TFType, UFType, F, JCType, GCType, RTolType, A} <:
+                      RosenbrockMutableCache
+    u::uType
+    uprev::uType
+    dense1::rateType
+    dense2::rateType
+    dense3::rateType
+    du::rateType
+    du1::rateType
+    du2::rateType
+    k1::rateType
+    k2::rateType
+    k3::rateType
+    k4::rateType
+    k5::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
+end
+
+@cache mutable struct Rodas3PCache{uType, rateType, uNoUnitsType, JType, WType, TabType,
+    TFType, UFType, F, JCType, GCType, RTolType, A} <:
+                      RosenbrockMutableCache
+    u::uType
+    uprev::uType
+    dense1::rateType
+    dense2::rateType
+    dense3::rateType
+    du::rateType
+    du1::rateType
+    du2::rateType
+    k1::rateType
+    k2::rateType
+    k3::rateType
+    k4::rateType
+    k5::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
+end
+
+function alg_cache(alg::Rodas23W, u, rate_prototype, ::Type{uEltypeNoUnits},
+    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+    dt, reltol, p, calck,
+    ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    dense1 = zero(rate_prototype)
+    dense2 = zero(rate_prototype)
+    dense3 = zero(rate_prototype)
+    du = zero(rate_prototype)
+    du1 = zero(rate_prototype)
+    du2 = zero(rate_prototype)
+    k1 = zero(rate_prototype)
+    k2 = zero(rate_prototype)
+    k3 = zero(rate_prototype)
+    k4 = zero(rate_prototype)
+    k5 = zero(rate_prototype)
+    fsalfirst = zero(rate_prototype)
+    fsallast = zero(rate_prototype)
+    dT = zero(rate_prototype)
+    J, W = build_J_W(alg, u, uprev, p, t, dt, f, uEltypeNoUnits, Val(true))
+    tmp = zero(rate_prototype)
+    atmp = similar(u, uEltypeNoUnits)
+    recursivefill!(atmp, false)
+    weight = similar(u, uEltypeNoUnits)
+    recursivefill!(weight, false)
+    tab = Rodas3PTableau(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
+
+    tf = TimeGradientWrapper(f, uprev, p)
+    uf = UJacobianWrapper(f, t, p)
+    linsolve_tmp = zero(rate_prototype)
+    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_A = true, alias_b = true,
+        Pl = Pl, Pr = Pr,
+        assumptions = LinearSolve.OperatorAssumptions(true))
+    grad_config = build_grad_config(alg, f, tf, du1, t)
+    jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2)
+    Rodas23WCache(u, uprev, dense1, dense2, dense3, du, du1, du2, k1, k2, k3, k4, k5,
+        fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf, linsolve_tmp,
+        linsolve, jac_config, grad_config, reltol, alg)
+end
+
+TruncatedStacktraces.@truncate_stacktrace Rodas23WCache 1
+function alg_cache(alg::Rodas3P, u, rate_prototype, ::Type{uEltypeNoUnits},
+    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+    dt, reltol, p, calck,
+    ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    dense1 = zero(rate_prototype)
+    dense2 = zero(rate_prototype)
+    dense3 = zero(rate_prototype)
+    du = zero(rate_prototype)
+    du1 = zero(rate_prototype)
+    du2 = zero(rate_prototype)
+    k1 = zero(rate_prototype)
+    k2 = zero(rate_prototype)
+    k3 = zero(rate_prototype)
+    k4 = zero(rate_prototype)
+    k5 = zero(rate_prototype)
+    fsalfirst = zero(rate_prototype)
+    fsallast = zero(rate_prototype)
+    dT = zero(rate_prototype)
+    J, W = build_J_W(alg, u, uprev, p, t, dt, f, uEltypeNoUnits, Val(true))
+    tmp = zero(rate_prototype)
+    atmp = similar(u, uEltypeNoUnits)
+    recursivefill!(atmp, false)
+    weight = similar(u, uEltypeNoUnits)
+    recursivefill!(weight, false)
+    tab = Rodas3PTableau(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
+
+    tf = TimeGradientWrapper(f, uprev, p)
+    uf = UJacobianWrapper(f, t, p)
+    linsolve_tmp = zero(rate_prototype)
+    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_A = true, alias_b = true,
+        Pl = Pl, Pr = Pr,
+        assumptions = LinearSolve.OperatorAssumptions(true))
+    grad_config = build_grad_config(alg, f, tf, du1, t)
+    jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2)
+    Rodas3PCache(u, uprev, dense1, dense2, dense3, du, du1, du2, k1, k2, k3, k4, k5,
+        fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf, linsolve_tmp,
+        linsolve, jac_config, grad_config, reltol, alg)
+end
+
+TruncatedStacktraces.@truncate_stacktrace Rodas3PCache 1
+
+function alg_cache(alg::Rodas23W, u, rate_prototype, ::Type{uEltypeNoUnits},
+    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+    dt, reltol, p, calck,
+    ::Val{false}) 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, uEltypeNoUnits, Val(false))
+    linprob = nothing #LinearProblem(W,copy(u); u0=copy(u))
+    linsolve = nothing #init(linprob,alg.linsolve,alias_A=true,alias_b=true)
+    Rodas23WConstantCache(tf, uf,
+        Rodas3PTableau(constvalue(uBottomEltypeNoUnits),
+            constvalue(tTypeNoUnits)), J, W, linsolve,
+        alg_autodiff(alg))
+end
+
+function alg_cache(alg::Rodas3P, u, rate_prototype, ::Type{uEltypeNoUnits},
+    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+    dt, reltol, p, calck,
+    ::Val{false}) 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, uEltypeNoUnits, Val(false))
+    linprob = nothing #LinearProblem(W,copy(u); u0=copy(u))
+    linsolve = nothing #init(linprob,alg.linsolve,alias_A=true,alias_b=true)
+    Rodas3PConstantCache(tf, uf,
+        Rodas3PTableau(constvalue(uBottomEltypeNoUnits),
+            constvalue(tTypeNoUnits)), J, W, linsolve,
+        alg_autodiff(alg))
+end
+
+### Rodas4 methods
+
 struct Rodas4ConstantCache{TF, UF, Tab, JType, WType, F, AD} <: OrdinaryDiffEqConstantCache
     tf::TF
     uf::UF

diff --git a/src/dense/rosenbrock_interpolants.jl b/src/dense/rosenbrock_interpolants.jl
@@ -2,6 +2,8 @@
 
 ROSENBROCKS_WITH_INTERPOLATIONS = Union{Rosenbrock23ConstantCache, Rosenbrock23Cache,
                                         Rosenbrock32ConstantCache, Rosenbrock32Cache,
+                                        Rodas23WConstantCache, Rodas3PConstantCache, 
+                                        Rodas23WCache, Rodas3PCache, 
                                         Rodas4ConstantCache, Rosenbrock5ConstantCache,
                                         Rodas4Cache, Rosenbrock5Cache}
 
@@ -134,34 +136,35 @@ end
 """
 From MATLAB ODE Suite by Shampine
 """
-@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rodas4ConstantCache,
+@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Rodas4ConstantCache, Rodas23WConstantCache, Rodas3PConstantCache},
     idxs::Nothing, T::Type{Val{0}}, differential_vars)
     Θ1 = 1 - Θ
     @inbounds Θ1 * y₀ + Θ * (y₁ + Θ1 * (k[1] + Θ * k[2]))
 end
 
-@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rodas4Cache, idxs::Nothing,
-    T::Type{Val{0}}, differential_vars)
+@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Rodas4Cache, Rodas23WCache, Rodas3PCache}, 
+    idxs::Nothing, T::Type{Val{0}}, differential_vars)
     Θ1 = 1 - Θ
     @inbounds @.. broadcast=false Θ1 * y₀+Θ * (y₁ + Θ1 * (k[1] + Θ * k[2]))
 end
 
 @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k,
-    cache::Union{Rodas4ConstantCache, Rodas4Cache}, idxs,
-    T::Type{Val{0}}, differential_vars)
+    cache::Union{Rodas4ConstantCache, Rodas4Cache, Rodas23WConstantCache, Rodas23WCache, Rodas3PConstantCache, Rodas3PCache}, 
+    idxs,  T::Type{Val{0}}, differential_vars)
     Θ1 = 1 - Θ
     @.. broadcast=false Θ1 * y₀[idxs]+Θ * (y₁[idxs] + Θ1 * (k[1][idxs] + Θ * k[2][idxs]))
 end
 
 @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k,
-    cache::Union{Rodas4ConstantCache, Rodas4Cache},
+    cache::Union{Rodas4ConstantCache, Rodas4Cache, Rodas23WConstantCache, Rodas23WCache, Rodas3PConstantCache, Rodas3PCache},
     idxs::Nothing, T::Type{Val{0}}, differential_vars)
     Θ1 = 1 - Θ
     @.. broadcast=false out=Θ1 * y₀ + Θ * (y₁ + Θ1 * (k[1] + Θ * k[2]))
     out
 end
 
-@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Rodas4Cache{<:Array},
+@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, 
+    cache::Union{Rodas4Cache{<:Array}, Rodas23WCache{<:Array}, Rodas3PCache{<:Array}},
     idxs::Nothing, T::Type{Val{0}}, differential_vars)
     Θ1 = 1 - Θ
     @inbounds @simd ivdep for i in eachindex(out)
@@ -171,45 +174,45 @@ end
 end
 
 @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k,
-    cache::Union{Rodas4ConstantCache, Rodas4Cache}, idxs,
-    T::Type{Val{0}}, differential_vars)
+    cache::Union{Rodas4ConstantCache, Rodas4Cache, Rodas23WConstantCache, Rodas23WCache, Rodas3PConstantCache, Rodas3PCache}, 
+    idxs, T::Type{Val{0}}, differential_vars)
     Θ1 = 1 - Θ
     @views @.. broadcast=false out=Θ1 * y₀[idxs] +
                                    Θ * (y₁[idxs] + Θ1 * (k[1][idxs] + Θ * k[2][idxs]))
     out
 end
 
 # First Derivative
-@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rodas4ConstantCache,
+@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Rodas4ConstantCache, Rodas23WConstantCache, Rodas3PConstantCache},
     idxs::Nothing, T::Type{Val{1}}, differential_vars)
     @inbounds (k[1] + Θ * (-2 * k[1] + 2 * k[2] - 3 * k[2] * Θ) - y₀ + y₁) / dt
 end
 
-@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rodas4Cache, idxs::Nothing,
-    T::Type{Val{1}}, differential_vars)
+@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Rodas4Cache, Rodas23WCache, Rodas3PCache}, 
+    idxs::Nothing,  T::Type{Val{1}}, differential_vars)
     @inbounds @.. broadcast=false (k[1] + Θ * (-2 * k[1] + 2 * k[2] - 3 * k[2] * Θ) - y₀ +
                                    y₁)/dt
 end
 
 @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k,
-    cache::Union{Rodas4ConstantCache, Rodas4Cache}, idxs,
-    T::Type{Val{1}}, differential_vars)
+    cache::Union{Rodas4ConstantCache, Rodas4Cache, Rodas23WConstantCache, Rodas23WCache, Rodas3PConstantCache, Rodas3PCache}, 
+    idxs, T::Type{Val{1}}, differential_vars)
     @.. broadcast=false (k[1][idxs] +
                          Θ * (-2 * k[1][idxs] + 2 * k[2][idxs] - 3 * k[2][idxs] * Θ) -
                          y₀[idxs] + y₁[idxs])/dt
 end
 
 @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k,
-    cache::Union{Rodas4ConstantCache, Rodas4Cache},
+    cache::Union{Rodas4ConstantCache, Rodas4Cache, Rodas23WConstantCache, Rodas23WCache, Rodas3PConstantCache, Rodas3PCache},
     idxs::Nothing, T::Type{Val{1}}, differential_vars)
     @.. broadcast=false out=(k[1] + Θ * (-2 * k[1] + 2 * k[2] - 3 * k[2] * Θ) - y₀ + y₁) /
                             dt
     out
 end
 
 @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k,
-    cache::Union{Rodas4ConstantCache, Rodas4Cache}, idxs,
-    T::Type{Val{1}}, differential_vars)
+    cache::Union{Rodas4ConstantCache, Rodas4Cache, Rodas23WConstantCache, Rodas23WCache, Rodas3PConstantCache, Rodas3PCache}, 
+    idxs, T::Type{Val{1}}, differential_vars)
     @views @.. broadcast=false out=(k[1][idxs] +
                                     Θ *
                                     (-2 * k[1][idxs] + 2 * k[2][idxs] -