Added ExponentialRK, IMEXMultistep and Linear Solvers

.github/workflows/CI.yml

@@ -33,6 +33,7 @@ jobs:
           - SSPRK
           - LowStorageRK
           - QPRK
+          - Linear
         version:
           - '1'
     steps:

lib/OrdinaryDiffEqExponentialRK/Project.toml

@@ -0,0 +1,25 @@
+name = "OrdinaryDiffEqExponentialRK"
+uuid = "e0540318-69ee-4070-8777-9e2de6de23de"
+authors = ["ParamThakkar123 <[email protected]>"]
+version = "1.0.0"
+
+[deps]
+FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
+MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+
+[compat]
+julia = "1.10"
+
+[extras]
+DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"]

lib/OrdinaryDiffEqExponentialRK/src/OrdinaryDiffEqExponentialRK.jl

@@ -0,0 +1,23 @@
+module OrdinaryDiffEqExponentialRK
+
+import OrdinaryDiffEq: alg_order, alg_adaptive_order, ismultistep, OrdinaryDiffEqExponentialAlgorithm,
+                       _unwrap_val, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache,
+                       build_jac_config, UJacobianWrapper, @cache, alg_cache, UDerivativeWrapper,
+                       initialize!, perform_step!, @unpack, unwrap_alg, calc_J, calc_J!,
+                       OrdinaryDiffEqAdaptiveExponentialAlgorithm, CompositeAlgorithm,
+                       ExponentialAlgorithm, fsal_typeof, isdtchangeable
+using RecursiveArrayTools
+using MuladdMacro, FastBroadcast
+using LinearAlgebra: axpy!, mul!
+using DiffEqBase, SciMLBase
+using ExponentialUtilities
+
+include("algorithms.jl")
+include("alg_utils.jl")
+include("exponential_rk_caches.jl")
+include("exponential_rk_perform_step.jl")
+
+export LawsonEuler, NorsettEuler, ETD1, ETDRK2, ETDRK3, ETDRK4, HochOst4, Exp4, EPIRK4s3A,
+       EPIRK4s3B,
+       EPIRK5s3, EXPRB53s3, EPIRK5P1, EPIRK5P2, ETD2, Exprb32, Exprb43
+end

lib/OrdinaryDiffEqExponentialRK/src/alg_utils.jl

@@ -0,0 +1,39 @@
+function isdtchangeable(alg::Union{LawsonEuler, NorsettEuler, ETDRK2, ETDRK3, ETDRK4, HochOst4, ETD2})
+    false
+end # due to caching
+
+alg_order(alg::LawsonEuler) = 1
+alg_order(alg::NorsettEuler) = 1
+alg_order(alg::ETDRK2) = 2
+alg_order(alg::ETDRK3) = 3
+alg_order(alg::ETDRK4) = 4
+alg_order(alg::HochOst4) = 4
+alg_order(alg::Exp4) = 4
+alg_order(alg::EPIRK4s3A) = 4
+alg_order(alg::EPIRK4s3B) = 4
+alg_order(alg::EPIRK5s3) = 5
+alg_order(alg::EPIRK5P1) = 5
+alg_order(alg::EPIRK5P2) = 5
+alg_order(alg::EXPRB53s3) = 5
+alg_order(alg::ETD2) = 2
+alg_order(alg::Exprb32) = 3
+alg_order(alg::Exprb43) = 4
+
+alg_adaptive_order(alg::Exprb32) = 2
+alg_adaptive_order(alg::Exprb43) = 4
+
+function DiffEqBase.prepare_alg(
+        alg::ETD2,
+        u0::AbstractArray,
+        p, prob)
+    alg
+end
+
+fsal_typeof(alg::ETD2, rate_prototype) = ETD2Fsal{typeof(rate_prototype)}
+function fsal_typeof(alg::CompositeAlgorithm, rate_prototype)
+    fsal = map(x -> fsal_typeof(x, rate_prototype), alg.algs)
+    @assert length(unique(fsal))==1 "`fsal_typeof` must be consistent"
+    return fsal[1]
+end
+
+ismultistep(alg::ETD2) = true

lib/OrdinaryDiffEqExponentialRK/src/algorithms.jl

@@ -0,0 +1,62 @@
+for Alg in [:LawsonEuler, :NorsettEuler, :ETDRK2, :ETDRK3, :ETDRK4, :HochOst4]
+    """
+    Hochbruck, Marlis, and Alexander Ostermann. “Exponential Integrators.” Acta
+      Numerica 19 (2010): 209–286. doi:10.1017/S0962492910000048.
+    """
+    @eval struct $Alg{CS, AD, FDT, ST, CJ} <:
+                 OrdinaryDiffEqExponentialAlgorithm{CS, AD, FDT, ST, CJ}
+        krylov::Bool
+        m::Int
+        iop::Int
+    end
+    @eval function $Alg(; krylov = false, m = 30, iop = 0, autodiff = true,
+            standardtag = Val{true}(), concrete_jac = nothing,
+            chunk_size = Val{0}(),
+            diff_type = Val{:forward})
+        $Alg{_unwrap_val(chunk_size), _unwrap_val(autodiff),
+            diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac)}(krylov,
+            m,
+            iop)
+    end
+end
+
+const ETD1 = NorsettEuler # alias
+
+for Alg in [:Exprb32, :Exprb43]
+    @eval struct $Alg{CS, AD, FDT, ST, CJ} <:
+                 OrdinaryDiffEqAdaptiveExponentialAlgorithm{CS, AD, FDT, ST, CJ}
+        m::Int
+        iop::Int
+    end
+    @eval function $Alg(; m = 30, iop = 0, autodiff = true, standardtag = Val{true}(),
+            concrete_jac = nothing, chunk_size = Val{0}(),
+            diff_type = Val{:forward})
+        $Alg{_unwrap_val(chunk_size), _unwrap_val(autodiff),
+            diff_type, _unwrap_val(standardtag),
+            _unwrap_val(concrete_jac)}(m,
+            iop)
+    end
+end
+for Alg in [:Exp4, :EPIRK4s3A, :EPIRK4s3B, :EPIRK5s3, :EXPRB53s3, :EPIRK5P1, :EPIRK5P2]
+    @eval struct $Alg{CS, AD, FDT, ST, CJ} <:
+                 OrdinaryDiffEqExponentialAlgorithm{CS, AD, FDT, ST, CJ}
+        adaptive_krylov::Bool
+        m::Int
+        iop::Int
+    end
+    @eval function $Alg(; adaptive_krylov = true, m = 30, iop = 0, autodiff = true,
+            standardtag = Val{true}(), concrete_jac = nothing,
+            chunk_size = Val{0}(), diff_type = Val{:forward})
+        $Alg{_unwrap_val(chunk_size), _unwrap_val(autodiff), diff_type,
+            _unwrap_val(standardtag), _unwrap_val(concrete_jac)}(adaptive_krylov,
+            m,
+            iop)
+    end
+end
+
+"""
+ETD2: Exponential Runge-Kutta Method
+Second order Exponential Time Differencing method (in development).
+"""
+struct ETD2 <:
+       OrdinaryDiffEqExponentialAlgorithm{0, false, Val{:forward}, Val{true}, nothing} end

src/caches/linear_nonlinear_caches.jl → lib/OrdinaryDiffEqExponentialRK/src/exponential_rk_caches.jl

@@ -901,4 +901,4 @@ function alg_cache(alg::ETD2, u, rate_prototype, ::Type{uEltypeNoUnits},
     ETD2Cache(
         u, uprev, zero(u), zero(rate_prototype), zero(rate_prototype), Phi[1], Phi[2],
         Phi[2] + Phi[3], -Phi[3])
-end
+end

src/perform_step/exponential_rk_perform_step.jl → lib/OrdinaryDiffEqExponentialRK/src/exponential_rk_perform_step.jl

@@ -1,5 +1,3 @@
-using LinearAlgebra: axpy!
-
 # Helper function to compute the G_nj factors for the classical ExpRK methods
 @inline _compute_nl(f::SplitFunction, u, p, t, A) = f.f2(u, p, t)
 @inline _compute_nl(f, u, p, t, A) = f(u, p, t) - A * u
@@ -1635,4 +1633,4 @@ function perform_step!(integrator, cache::ETD2Cache, repeat_step = false)
     integrator.stats.nf += 1
     integrator.stats.nf2 += 1
     @.. broadcast=false integrator.k[2]=fsallast.lin + fsallast.nl
-end
+end

test/algconvergence/linear_nonlinear_krylov_tests.jl → lib/OrdinaryDiffEqExponentialRK/test/linear_nonlinear_krylov_tests.jl

@@ -136,4 +136,4 @@ end
     prob = ODEProblem(exp_fun, u0, (0.0, 1.0))
     sol = solve(prob, LawsonEuler(krylov = true, m = N); dt = 0.1)
     @test sol(1.0) ≈ exp_fun.analytic(u0, nothing, 1.0)
-end
+end

lib/OrdinaryDiffEqExponentialRK/test/runtests.jl

@@ -0,0 +1,3 @@
+using SafeTestsets
+
+@time @safetestset "Linear-Nonlinear Krylov Methods Tests" include("linear_nonlinear_krylov_tests.jl")

lib/OrdinaryDiffEqIMEXMultistep/Project.toml

@@ -0,0 +1,20 @@
+name = "OrdinaryDiffEqIMEXMultistep"
+uuid = "9f002381-b378-40b7-97a6-27a27c83f129"
+authors = ["ParamThakkar123 <[email protected]>"]
+version = "1.0.0"
+
+[deps]
+FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+
+[compat]
+julia = "1.10"
+
+[extras]
+DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"]

lib/OrdinaryDiffEqIMEXMultistep/src/OrdinaryDiffEqIMEXMultistep.jl

@@ -0,0 +1,16 @@
+module OrdinaryDiffEqIMEXMultistep
+
+import OrdinaryDiffEq: alg_order, issplit, OrdinaryDiffEqNewtonAlgorithm, _unwrap_val, dolinsolve,
+                       DEFAULT_PRECS, NLNewton, OrdinaryDiffEqConstantCache, OrdinaryDiffEqMutableCache,
+                       build_nlsolver, @cache, alg_cache, initialize!, perform_step!, @unpack,
+                       markfirststage!, nlsolve!, nlsolvefail, du_alias_or_new
+using FastBroadcast
+
+include("algorithms.jl")
+include("alg_utils.jl")
+include("imex_multistep_caches.jl")
+include("imex_multistep_perform_step.jl")
+
+export CNAB2, CNLF2
+
+end

lib/OrdinaryDiffEqIMEXMultistep/src/alg_utils.jl

@@ -0,0 +1,4 @@
+issplit(alg::Union{CNAB2, CNLF2}) = true
+
+alg_order(alg::CNAB2) = 2
+alg_order(alg::CNLF2) = 2

lib/OrdinaryDiffEqIMEXMultistep/src/algorithms.jl

@@ -0,0 +1,42 @@
+# IMEX Multistep methods
+
+struct CNAB2{CS, AD, F, F2, P, FDT, ST, CJ} <:
+    OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ}
+ linsolve::F
+ nlsolve::F2
+ precs::P
+ extrapolant::Symbol
+end
+
+function CNAB2(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(),
+     concrete_jac = nothing, diff_type = Val{:forward},
+     linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+     extrapolant = :linear)
+ CNAB2{
+     _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve),
+     typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac)}(
+     linsolve,
+     nlsolve,
+     precs,
+     extrapolant)
+end
+
+struct CNLF2{CS, AD, F, F2, P, FDT, ST, CJ} <:
+    OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ}
+ linsolve::F
+ nlsolve::F2
+ precs::P
+ extrapolant::Symbol
+end
+function CNLF2(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(),
+     concrete_jac = nothing, diff_type = Val{:forward},
+     linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+     extrapolant = :linear)
+ CNLF2{
+     _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve),
+     typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac)}(
+     linsolve,
+     nlsolve,
+     precs,
+     extrapolant)
+end

lib/OrdinaryDiffEqLinear/Project.toml

@@ -0,0 +1,24 @@
+name = "OrdinaryDiffEqLinear"
+uuid = "521117fe-8c41-49f8-b3b6-30780b3f0fb5"
+authors = ["ParamThakkar123 <[email protected]>"]
+version = "1.0.0"
+
+[deps]
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18"
+
+[compat]
+julia = "1.10"
+
+[extras]
+DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"]

lib/OrdinaryDiffEqLinear/src/OrdinaryDiffEqLinear.jl

@@ -0,0 +1,25 @@
+module OrdinaryDiffEqLinear
+
+import OrdinaryDiffEq: alg_order, alg_extrapolates, dt_required, OrdinaryDiffEqLinearExponentialAlgorithm,
+                       OrdinaryDiffEqAdaptiveAlgorithm, OrdinaryDiffEqAlgorithm,
+                       OrdinaryDiffEqExponentialAlgorithm,
+                       OrdinaryDiffEqMutableCache, @cache, alg_cache, OrdinaryDiffEqConstantCache,
+                       initialize!, perform_step!, @unpack, unwrap_alg, calculate_residuals!,
+                       _vec, isdtchangeable
+using LinearAlgebra: mul!
+using DiffEqBase
+using SciMLOperators: AbstractSciMLOperator
+using ExponentialUtilities
+using RecursiveArrayTools
+
+include("algorithms.jl")
+include("alg_utils.jl")
+include("linear_caches.jl")
+include("integrator_interface.jl")
+include("linear_perform_step.jl")
+
+export MagnusMidpoint, LinearExponential, MagnusLeapfrog, LieEuler, CayleyEuler,
+       MagnusGauss4, MagnusNC6, MagnusGL6, MagnusGL8, MagnusNC8, MagnusGL4,
+       MagnusAdapt4, RKMK2, RKMK4, LieRK4, CG2, CG3, CG4a
+
+end

lib/OrdinaryDiffEqLinear/src/alg_utils.jl

@@ -0,0 +1,33 @@
+alg_order(alg::MagnusMidpoint) = 2
+alg_order(alg::RKMK2) = 2
+alg_order(alg::RKMK4) = 4
+alg_order(alg::LieRK4) = 4
+alg_order(alg::CG3) = 3
+alg_order(alg::CG2) = 2
+alg_order(alg::CG4a) = 4
+alg_order(alg::MagnusAdapt4) = 4
+alg_order(alg::MagnusGauss4) = 4
+alg_order(alg::MagnusNC6) = 6
+alg_order(alg::MagnusGL6) = 6
+alg_order(alg::MagnusGL8) = 8
+alg_order(alg::MagnusNC8) = 8
+alg_order(alg::MagnusGL4) = 4
+alg_order(alg::LinearExponential) = 1
+alg_order(alg::MagnusLeapfrog) = 2
+alg_order(alg::LieEuler) = 1
+alg_order(alg::CayleyEuler) = 2
+
+alg_extrapolates(alg::MagnusLeapfrog) = true
+
+dt_required(alg::LinearExponential) = false
+
+function DiffEqBase.prepare_alg(
+        alg::LinearExponential,
+        u0::AbstractArray,
+        p, prob)
+    alg
+end
+
+function isdtchangeable(alg::Union{LieEuler, MagnusGauss4, CayleyEuler})
+    false
+end # due to caching