Added High Order and Low Order RK methods

Project.toml

@@ -18,7 +18,6 @@ FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf"
-IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -60,7 +59,6 @@ FillArrays = "1.9"
 FiniteDiff = "2"
 ForwardDiff = "0.10.3"
 FunctionWrappersWrappers = "0.1"
-IfElse = "0.1"
 InteractiveUtils = "1.9"
 LineSearches = "7"
 LinearAlgebra = "1.9"

lib/OrdinaryDiffEqAdamsBashforthMoulton/Project.toml

@@ -0,0 +1,25 @@
+name = "OrdinaryDiffEqAdamsBashforthMoulton"
+uuid = "89bda076-bce5-4f1c-845f-551c83cdda9a"
+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"
+Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
+ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
+
+[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/OrdinaryDiffEqAdamsBashforthMoulton/src/OrdinaryDiffEqAdamsBashforthMoulton.jl

@@ -0,0 +1,22 @@
+module OrdinaryDiffEqAdamsBashforthMoulton
+
+import OrdinaryDiffEq: OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, @cache, alg_cache,
+                       initialize!, @unpack, perform_step!, alg_order, isstandard, OrdinaryDiffEqAlgorithm,
+                       OrdinaryDiffEqAdaptiveAlgorithm, OrdinaryDiffEqAdamsVarOrderVarStepAlgorithm,
+                       constvalue, calculate_residuals, calculate_residuals!, trivial_limiter!
+import OrdinaryDiffEq: BS3ConstantCache, BS3Cache, RK4ConstantCache, RK4Cache
+import RecursiveArrayTools: recursivefill!
+using MuladdMacro, FastBroadcast
+import Static: False
+import ADTypes: AutoForwardDiff
+
+include("algorithms.jl")
+include("alg_utils.jl")
+include("adams_bashforth_moulton_caches.jl")
+include("adams_utils.jl")
+include("adams_bashforth_moulton_perform_step.jl")
+
+export AB3, AB4, AB5, ABM32, ABM43, ABM54, VCAB3,
+       VCAB4, VCAB5, VCABM3, VCABM4, VCABM5, VCABM
+
+end

src/caches/adams_bashforth_moulton_caches.jl → lib/OrdinaryDiffEqAdamsBashforthMoulton/src/adams_bashforth_moulton_caches.jl

@@ -1019,120 +1019,3 @@ function alg_cache(alg::VCABM, u, rate_prototype, ::Type{uEltypeNoUnits},
         max_order, atmp, tmp, ξ, ξ0, utilde, utildem1, utildem2, utildep1, atmpm1,
         atmpm2, atmpp1, 1)
 end
-
-# IMEX Multistep methods
-
-# CNAB2
-
-@cache mutable struct CNAB2ConstantCache{rateType, N, uType, tType} <:
-                      OrdinaryDiffEqConstantCache
-    k2::rateType
-    nlsolver::N
-    uprev3::uType
-    tprev2::tType
-end
-
-@cache mutable struct CNAB2Cache{uType, rateType, N, tType} <: OrdinaryDiffEqMutableCache
-    u::uType
-    uprev::uType
-    uprev2::uType
-    fsalfirst::rateType
-    k1::rateType
-    k2::rateType
-    du₁::rateType
-    nlsolver::N
-    uprev3::uType
-    tprev2::tType
-end
-
-function alg_cache(alg::CNAB2, u, rate_prototype, ::Type{uEltypeNoUnits},
-        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
-        dt, reltol, p, calck,
-        ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
-    γ, c = 1 // 2, 1
-    nlsolver = build_nlsolver(alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits,
-        uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(false))
-
-    k2 = rate_prototype
-    uprev3 = u
-    tprev2 = t
-
-    CNAB2ConstantCache(k2, nlsolver, uprev3, tprev2)
-end
-
-function alg_cache(alg::CNAB2, u, rate_prototype, ::Type{uEltypeNoUnits},
-        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
-        dt, reltol, p, calck,
-        ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
-    γ, c = 1 // 2, 1
-    nlsolver = build_nlsolver(alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits,
-        uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(true))
-    fsalfirst = zero(rate_prototype)
-
-    k1 = zero(rate_prototype)
-    k2 = zero(rate_prototype)
-    du₁ = zero(rate_prototype)
-    uprev3 = zero(u)
-    tprev2 = t
-
-    CNAB2Cache(u, uprev, uprev2, fsalfirst, k1, k2, du₁, nlsolver, uprev3, tprev2)
-end
-
-# CNLF2
-
-@cache mutable struct CNLF2ConstantCache{rateType, N, uType, tType} <:
-                      OrdinaryDiffEqConstantCache
-    k2::rateType
-    nlsolver::N
-    uprev2::uType
-    uprev3::uType
-    tprev2::tType
-end
-
-@cache mutable struct CNLF2Cache{uType, rateType, N, tType} <: OrdinaryDiffEqMutableCache
-    u::uType
-    uprev::uType
-    uprev2::uType
-    fsalfirst::rateType
-    k1::rateType
-    k2::rateType
-    du₁::rateType
-    nlsolver::N
-    uprev3::uType
-    tprev2::tType
-end
-
-function alg_cache(alg::CNLF2, u, rate_prototype, ::Type{uEltypeNoUnits},
-        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
-        dt, reltol, p, calck,
-        ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
-    γ, c = 1 // 1, 1
-    nlsolver = build_nlsolver(alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits,
-        uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(false))
-
-    k2 = rate_prototype
-    uprev2 = u
-    uprev3 = u
-    tprev2 = t
-
-    CNLF2ConstantCache(k2, nlsolver, uprev2, uprev3, tprev2)
-end
-
-function alg_cache(alg::CNLF2, u, rate_prototype, ::Type{uEltypeNoUnits},
-        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
-        dt, reltol, p, calck,
-        ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
-    γ, c = 1 // 1, 1
-    nlsolver = build_nlsolver(alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits,
-        uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(true))
-    fsalfirst = zero(rate_prototype)
-
-    k1 = zero(rate_prototype)
-    k2 = zero(rate_prototype)
-    du₁ = zero(rate_prototype)
-    uprev2 = zero(u)
-    uprev3 = zero(u)
-    tprev2 = t
-
-    CNLF2Cache(u, uprev, uprev2, fsalfirst, k1, k2, du₁, nlsolver, uprev3, tprev2)
-end

src/perform_step/adams_bashforth_moulton_perform_step.jl → lib/OrdinaryDiffEqAdamsBashforthMoulton/src/adams_bashforth_moulton_perform_step.jl

@@ -1540,209 +1540,3 @@ end
         return nothing
     end # inbounds
 end
-
-# CNAB2
-
-function initialize!(integrator, cache::CNAB2ConstantCache)
-    integrator.kshortsize = 2
-    integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
-    integrator.fsalfirst = integrator.f.f1(integrator.uprev, integrator.p, integrator.t) +
-                           integrator.f.f2(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
-    integrator.stats.nf += 1
-    integrator.stats.nf2 += 1
-
-    # Avoid undefined entries if k is an array of arrays
-    integrator.fsallast = zero(integrator.fsalfirst)
-    integrator.k[1] = integrator.fsalfirst
-    integrator.k[2] = integrator.fsallast
-end
-
-function perform_step!(integrator, cache::CNAB2ConstantCache, repeat_step = false)
-    @unpack t, dt, uprev, u, f, p = integrator
-    @unpack k2, nlsolver = cache
-    cnt = integrator.iter
-    f1 = integrator.f.f1
-    f2 = integrator.f.f2
-    du₁ = f1(uprev, p, t)
-    integrator.stats.nf += 1
-    k1 = integrator.fsalfirst - du₁
-    # Explicit part
-    if cnt == 1
-        tmp = uprev + dt * k1
-    else
-        tmp = uprev + dt * (3 // 2 * k1 - 1 // 2 * k2)
-    end
-    nlsolver.tmp = tmp
-    # Implicit part
-    # precalculations
-    γ = 1 // 2
-    γdt = γ * dt
-
-    # initial guess
-    zprev = dt * du₁
-    nlsolver.z = z = zprev # Constant extrapolation
-
-    nlsolver.tmp += γ * zprev
-    markfirststage!(nlsolver)
-    z = nlsolve!(nlsolver, integrator, cache, repeat_step)
-    nlsolvefail(nlsolver) && return
-    u = nlsolver.tmp + 1 // 2 * z
-
-    cache.k2 = k1
-    integrator.fsallast = f1(u, p, t + dt) + f2(u, p, t + dt)
-    integrator.stats.nf += 1
-    integrator.stats.nf2 += 1
-    integrator.k[1] = integrator.fsalfirst
-    integrator.k[2] = integrator.fsallast
-    integrator.u = u
-end
-
-function initialize!(integrator, cache::CNAB2Cache)
-    integrator.kshortsize = 2
-    integrator.fsalfirst = cache.fsalfirst
-    integrator.fsallast = du_alias_or_new(cache.nlsolver, integrator.fsalfirst)
-    resize!(integrator.k, integrator.kshortsize)
-    integrator.k[1] = integrator.fsalfirst
-    integrator.k[2] = integrator.fsallast
-    integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t)
-    integrator.stats.nf += 1
-end
-
-function perform_step!(integrator, cache::CNAB2Cache, repeat_step = false)
-    @unpack t, dt, uprev, u, f, p = integrator
-    @unpack k1, k2, du₁, nlsolver = cache
-    @unpack z, tmp = nlsolver
-    @unpack f1 = f
-    cnt = integrator.iter
-
-    f1(du₁, uprev, p, t)
-    integrator.stats.nf += 1
-    @.. broadcast=false k1=integrator.fsalfirst - du₁
-    # Explicit part
-    if cnt == 1
-        @.. broadcast=false tmp=uprev + dt * k1
-    else
-        @.. broadcast=false tmp=uprev + dt * (3 // 2 * k1 - 1 // 2 * k2)
-    end
-    # Implicit part
-    # precalculations
-    γ = 1 // 2
-    γdt = γ * dt
-
-    # initial guess
-    @.. broadcast=false z=dt * du₁
-    @.. broadcast=false tmp+=γ * z
-    markfirststage!(nlsolver)
-    z = nlsolve!(nlsolver, integrator, cache, repeat_step)
-    nlsolvefail(nlsolver) && return
-    @.. broadcast=false u=tmp + 1 // 2 * z
-
-    cache.k2 .= k1
-    f(integrator.fsallast, u, p, t + dt)
-    integrator.stats.nf += 1
-end
-
-# CNLF2
-
-function initialize!(integrator, cache::CNLF2ConstantCache)
-    integrator.kshortsize = 2
-    integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
-    integrator.fsalfirst = integrator.f.f1(integrator.uprev, integrator.p, integrator.t) +
-                           integrator.f.f2(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
-    integrator.stats.nf += 1
-    integrator.stats.nf2 += 1
-
-    # Avoid undefined entries if k is an array of arrays
-    integrator.fsallast = zero(integrator.fsalfirst)
-    integrator.k[1] = integrator.fsalfirst
-    integrator.k[2] = integrator.fsallast
-end
-
-function perform_step!(integrator, cache::CNLF2ConstantCache, repeat_step = false)
-    @unpack t, dt, uprev, u, f, p = integrator
-    @unpack k2, uprev2, nlsolver = cache
-    cnt = integrator.iter
-    f1 = integrator.f.f1
-    f2 = integrator.f.f2
-    du₁ = f1(uprev, p, t)
-    integrator.stats.nf += 1
-    # Explicit part
-    if cnt == 1
-        tmp = uprev + dt * (integrator.fsalfirst - du₁)
-    else
-        tmp = uprev2 + 2 // 1 * dt * (integrator.fsalfirst - du₁)
-    end
-    # Implicit part
-    # precalculations
-    γ = 1 // 1
-    if cnt != 1
-        tmp += γ * dt * k2
-    end
-    γdt = γ * dt
-    nlsolver.tmp = tmp
-
-    # initial guess
-    zprev = dt * du₁
-    nlsolver.z = z = zprev # Constant extrapolation
-
-    markfirststage!(nlsolver)
-    z = nlsolve!(nlsolver, integrator, cache, repeat_step)
-    nlsolvefail(nlsolver) && return
-    u = nlsolver.tmp + γ * z
-
-    cache.uprev2 = uprev
-    cache.k2 = du₁
-    integrator.fsallast = f1(u, p, t + dt) + f2(u, p, t + dt)
-    integrator.stats.nf += 1
-    integrator.stats.nf2 += 1
-    integrator.k[1] = integrator.fsalfirst
-    integrator.k[2] = integrator.fsallast
-    integrator.u = u
-end
-
-function initialize!(integrator, cache::CNLF2Cache)
-    integrator.kshortsize = 2
-    integrator.fsalfirst = cache.fsalfirst
-    integrator.fsallast = du_alias_or_new(cache.nlsolver, integrator.fsalfirst)
-    resize!(integrator.k, integrator.kshortsize)
-    integrator.k[1] = integrator.fsalfirst
-    integrator.k[2] = integrator.fsallast
-    integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t)
-    integrator.stats.nf += 1
-end
-
-function perform_step!(integrator, cache::CNLF2Cache, repeat_step = false)
-    @unpack t, dt, uprev, u, f, p = integrator
-    @unpack uprev2, k2, du₁, nlsolver = cache
-    @unpack z, tmp = nlsolver
-    @unpack f1 = f
-    cnt = integrator.iter
-
-    f1(du₁, uprev, p, t)
-    integrator.stats.nf += 1
-    # Explicit part
-    if cnt == 1
-        @.. broadcast=false tmp=uprev + dt * (integrator.fsalfirst - du₁)
-    else
-        @.. broadcast=false tmp=uprev2 + 2 // 1 * dt * (integrator.fsalfirst - du₁)
-    end
-    # Implicit part
-    # precalculations
-    γ = 1 // 1
-    if cnt != 1
-        @.. broadcast=false tmp+=γ * dt * k2
-    end
-    γdt = γ * dt
-
-    # initial guess
-    @.. broadcast=false z=dt * du₁
-    markfirststage!(nlsolver)
-    z = nlsolve!(nlsolver, integrator, cache, repeat_step)
-    nlsolvefail(nlsolver) && return
-    @.. broadcast=false u=tmp + γ * z
-
-    cache.uprev2 .= uprev
-    cache.k2 .= du₁
-    f(integrator.fsallast, u, p, t + dt)
-    integrator.stats.nf += 1
-end