Vsvo embedded bdf

lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl

@@ -154,6 +154,6 @@ end
 
 export ABDF2, QNDF1, QBDF1, QNDF2, QBDF2, QNDF, QBDF, FBDF,
     SBDF, SBDF2, SBDF3, SBDF4, MEBDF2, IMEXEuler, IMEXEulerARK,
-    DABDF2, DImplicitEuler, DFBDF
+    DABDF2, DImplicitEuler, DFBDF, MOOSE234
 
 end

lib/OrdinaryDiffEqBDF/src/alg_utils.jl

@@ -40,4 +40,12 @@ alg_order(alg::DFBDF) = 1 #dummy value
 
 isfsal(alg::DImplicitEuler) = false
 
-has_stiff_interpolation(::Union{QNDF, FBDF, DFBDF}) = true
+has_stiff_interpolation(::Union{QNDF, FBDF, DFBDF, MOOSE234}) = true
+
+alg_extrapolates(alg::MOOSE234) = true
+alg_order(alg::MOOSE234) = 2
+isfsal(alg::MOOSE234) = false
+get_current_alg_order(alg::MOOSE234, cache) = cache.order
+get_current_adaptive_order(alg::MOOSE234, cache) = cache.order
+qsteady_min_default(alg::MOOSE234) = 9 // 10
+qsteady_max_default(alg::MOOSE234) = 2 // 1

lib/OrdinaryDiffEqBDF/src/algorithms.jl

@@ -15,10 +15,10 @@ function BDF_docstring(
         """ * "\n" * extra_keyword_default
 
     keyword_default_description = """
-        - `autodiff`: Uses [ADTypes.jl](https://sciml.github.io/ADTypes.jl/stable/) 
+        - `autodiff`: Uses [ADTypes.jl](https://sciml.github.io/ADTypes.jl/stable/)
             to specify whether to use automatic differentiation via
             [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) or finite
-            differencing via [FiniteDiff.jl](https://github.com/JuliaDiff/FiniteDiff.jl). 
+            differencing via [FiniteDiff.jl](https://github.com/JuliaDiff/FiniteDiff.jl).
             Defaults to `AutoForwardDiff()` for automatic differentiation, which by default uses
             `chunksize = 0`, and thus uses the internal ForwardDiff.jl algorithm for the choice.
             To use `FiniteDiff.jl`, the `AutoFiniteDiff()` ADType can be used, which has a keyword argument
@@ -630,6 +630,63 @@ end
 
 @truncate_stacktrace FBDF
 
+@doc BDF_docstring(
+    "A variable-stepsize, variable-order (2/3/4) implicit method based on the MOOSE234
+(Multiple Order One Solve Embedded 234) scheme of DeCaria, Guzel, Layton, and Li.
+Performs a single BDF3 nonlinear solve per step, then applies two cheap linear
+time-filter post-processing steps to produce embedded solutions at orders 2 and 4.
+This yields error estimates at three orders, enabling VSVO behavior at the cost of
+a single implicit solve per step.",
+    "MOOSE234",
+    references = """@article{decaria2018variable,
+    title={A variable stepsize, variable order family of low complexity},
+    author={DeCaria, Victor and Guzel, Abigail and Layton, William and Li, Yanghong},
+    journal={arXiv preprint arXiv:1810.06670},
+    year={2018}}""",
+    extra_keyword_description = """
+    - `nlsolve`: TBD
+    - `extrapolant`: TBD
+    - `controller`: TBD
+    - `step_limiter!`: function of the form `limiter!(u, integrator, p, t)`
+    """,
+    extra_keyword_default = """
+    nlsolve = NLNewton(),
+    extrapolant = :linear,
+    controller = :Standard,
+    step_limiter! = trivial_limiter!,
+    """
+)
+struct MOOSE234{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
+    OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+    linsolve::F
+    nlsolve::F2
+    precs::P
+    extrapolant::Symbol
+    controller::Symbol
+    step_limiter!::StepLimiter
+    autodiff::AD
+end
+
+function MOOSE234(;
+        chunk_size = Val{0}(), autodiff = AutoForwardDiff(), standardtag = Val{true}(),
+        concrete_jac = nothing, diff_type = Val{:forward}(),
+        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        extrapolant = :linear, controller = :Standard, step_limiter! = trivial_limiter!
+    )
+    AD_choice, chunk_size, diff_type = _process_AD_choice(autodiff, chunk_size, diff_type)
+
+    return MOOSE234{
+        _unwrap_val(chunk_size), typeof(AD_choice), typeof(linsolve),
+        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac), typeof(step_limiter!),
+    }(
+        linsolve, nlsolve, precs, extrapolant,
+        controller, step_limiter!, AD_choice
+    )
+end
+
+@truncate_stacktrace MOOSE234
+
 """
 QBDF1: Multistep Method
 

lib/OrdinaryDiffEqBDF/src/bdf_caches.jl

@@ -792,3 +792,164 @@ function alg_cache(
         iters_from_event, dense, alg.step_limiter!, fd_weights, stald
     )
 end
+
+# MOOSE234 — Variable-stepsize, variable-order (2/3/4) embedded method
+# DeCaria et al., arXiv:1810.06670v1
+
+@cache mutable struct MOOSE234ConstantCache{
+        N, tsType, tType, uType, uuType,
+        EEstType, rType, wType,
+    } <:
+    OrdinaryDiffEqConstantCache
+    nlsolver::N
+    ts::tsType
+    ts_tmp::tsType
+    t_old::tType
+    u_history::uuType
+    order::Int
+    prev_order::Int
+    u_corrector::uType
+    nconsteps::Int
+    consfailcnt::Int
+    qwait::Int
+    terkm1::EEstType
+    terk::EEstType
+    terkp1::EEstType
+    r::rType
+    weights::wType
+    iters_from_event::Int
+end
+
+function alg_cache(
+        alg::MOOSE234, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+        dt, reltol, p, calck,
+        ::Val{false}, verbose
+    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    γ, c = Int64(1) // 1, 1
+    nlsolver = build_nlsolver(
+        alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits,
+        uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(false), verbose
+    )
+
+    hist_size = 6 # max_order(4) + 2
+    ts = zero(Vector{typeof(t)}(undef, hist_size))
+    ts_tmp = similar(ts)
+
+    if u isa Number
+        u_history = zeros(eltype(u), hist_size)
+        u_corrector = zeros(eltype(u), hist_size)
+    else
+        u_history = [zero(u) for _ in 1:hist_size]
+        u_corrector = [zero(u) for _ in 1:hist_size]
+    end
+
+    order = 2
+    prev_order = 2
+    terkm1 = tTypeNoUnits(1)
+    terk = tTypeNoUnits(1)
+    terkp1 = tTypeNoUnits(1)
+    r = zero(Vector{typeof(t)}(undef, hist_size))
+    weights = zero(Vector{typeof(t)}(undef, hist_size))
+    weights[1] = 1
+    nconsteps = 0
+    consfailcnt = 0
+    qwait = 4 # order + 2
+    t_old = zero(t)
+    iters_from_event = 0
+
+    return MOOSE234ConstantCache(
+        nlsolver, ts, ts_tmp, t_old, u_history, order, prev_order,
+        u_corrector, nconsteps, consfailcnt, qwait, terkm1,
+        terk, terkp1, r, weights, iters_from_event
+    )
+end
+
+@cache mutable struct MOOSE234Cache{
+        N, rateType, uNoUnitsType, tsType, tType, uType, uuType,
+        EEstType, rType, wType, StepLimiter, fdWeightsType,
+    } <:
+    BDFMutableCache
+    fsalfirst::rateType
+    nlsolver::N
+    ts::tsType
+    ts_tmp::tsType
+    t_old::tType
+    u_history::uuType
+    order::Int
+    prev_order::Int
+    u_corrector::uuType
+    u₀::uType
+    nconsteps::Int
+    consfailcnt::Int
+    qwait::Int
+    tmp::uType
+    atmp::uNoUnitsType
+    terkm1::EEstType
+    terk::EEstType
+    terkp1::EEstType
+    terk_tmp::uType
+    terkp1_tmp::uType
+    r::rType
+    weights::wType
+    equi_ts::tsType
+    iters_from_event::Int
+    dense::Vector{uType}
+    step_limiter!::StepLimiter
+    fd_weights::fdWeightsType
+end
+
+@truncate_stacktrace MOOSE234Cache 1
+
+function alg_cache(
+        alg::MOOSE234, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+        dt, reltol, p, calck,
+        ::Val{true}, verbose
+    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    γ, c = Int64(1) // 1, 1
+    fsalfirst = zero(rate_prototype)
+    nlsolver = build_nlsolver(
+        alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits,
+        uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(true), verbose
+    )
+
+    hist_size = 6
+    ts = Vector{typeof(t)}(undef, hist_size)
+    u_history = [zero(u) for _ in 1:hist_size]
+    order = 2
+    prev_order = 2
+    u_corrector = [zero(u) for _ in 1:hist_size]
+    recursivefill!(ts, zero(t))
+    terkm1 = tTypeNoUnits(1)
+    terk = tTypeNoUnits(1)
+    terkp1 = tTypeNoUnits(1)
+    terk_tmp = similar(u)
+    terkp1_tmp = similar(u)
+    r = Vector{typeof(t)}(undef, hist_size)
+    weights = Vector{typeof(t)}(undef, hist_size)
+    recursivefill!(r, zero(t))
+    recursivefill!(weights, zero(t))
+    weights[1] = 1
+    nconsteps = 0
+    consfailcnt = 0
+    qwait = 4 # order + 2
+    t_old = zero(t)
+    atmp = similar(u, uEltypeNoUnits)
+    recursivefill!(atmp, zero(uEltypeNoUnits))
+    u₀ = similar(u)
+    equi_ts = similar(ts)
+    tmp = similar(u)
+    ts_tmp = similar(ts)
+    iters_from_event = 0
+
+    dense = [zero(u) for _ in 1:10] # 2*(max_order+1) = 2*5 = 10
+    fd_weights = zeros(typeof(t), 5, 5) # (max_order+1) × (max_order+1)
+
+    return MOOSE234Cache(
+        fsalfirst, nlsolver, ts, ts_tmp, t_old, u_history, order, prev_order,
+        u_corrector, u₀, nconsteps, consfailcnt, qwait, tmp, atmp,
+        terkm1, terk, terkp1, terk_tmp, terkp1_tmp, r, weights, equi_ts,
+        iters_from_event, dense, alg.step_limiter!, fd_weights
+    )
+end