Add NordsieckBDF and DNordsieckBDF: Newton-based Nordsieck BDF solvers

lib/OrdinaryDiffEqBDF/Project.toml

@@ -33,6 +33,7 @@ LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
 DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 ODEProblemLibrary = "fdc4e326-1af4-4b90-96e7-779fcce2daa5"
+OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -43,6 +44,7 @@ DiffEqDevTools = {path = "../DiffEqDevTools"}
 OrdinaryDiffEqCore = {path = "../OrdinaryDiffEqCore"}
 OrdinaryDiffEqDifferentiation = {path = "../OrdinaryDiffEqDifferentiation"}
 OrdinaryDiffEqNonlinearSolve = {path = "../OrdinaryDiffEqNonlinearSolve"}
+OrdinaryDiffEqRosenbrock = {path = "../OrdinaryDiffEqRosenbrock"}
 OrdinaryDiffEqSDIRK = {path = "../OrdinaryDiffEqSDIRK"}
 
 [compat]
@@ -72,11 +74,12 @@ julia = "1.10"
 ADTypes = "1.22.0"
 RecursiveArrayTools = "4.2.0"
 ODEProblemLibrary = "1"
+OrdinaryDiffEqRosenbrock = "2"
 OrdinaryDiffEqNonlinearSolve = "2"
 DiffEqBase = "7"
 Reexport = "1.2.2"
 SafeTestsets = "0.1.0"
 StaticArrays = "1.9.18"
 
 [targets]
-test = ["DiffEqDevTools", "DifferentiationInterface", "ForwardDiff", "Random", "SafeTestsets", "StaticArrays", "Test", "ODEProblemLibrary", "NonlinearSolve", "Enzyme", "LinearSolve", "Pkg", "SciMLTesting"]
+test = ["DiffEqDevTools", "DifferentiationInterface", "ForwardDiff", "Random", "SafeTestsets", "StaticArrays", "Test", "ODEProblemLibrary", "NonlinearSolve", "Enzyme", "LinearSolve", "Pkg", "SciMLTesting", "OrdinaryDiffEqRosenbrock"]

lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl

@@ -66,6 +66,7 @@ include("algorithms.jl")
 include("alg_utils.jl")
 include("bdf_utils.jl")
 include("stald.jl")
+include("nordsieck_utils.jl")
 include("bdf_caches.jl")
 include("dae_caches.jl")
 include("controllers.jl")
@@ -90,7 +91,7 @@ end
 PrecompileTools.@compile_workload begin
     lorenz = OrdinaryDiffEqCore.lorenz
     lorenz_oop = OrdinaryDiffEqCore.lorenz_oop
-    solver_list = [FBDF()]
+    solver_list = [FBDF(), NordsieckBDF()]
     prob_list = []
 
     if Preferences.@load_preference("PrecompileDefaultSpecialize", true)
@@ -182,6 +183,7 @@ end
 
 export ABDF2, QNDF1, QBDF1, QNDF2, QBDF2, QNDF, QBDF, FBDF,
     SBDF, SBDF2, SBDF3, SBDF4, MEBDF2, IMEXEuler, IMEXEulerARK,
-    DABDF2, DImplicitEuler, DFBDF
+    DABDF2, DImplicitEuler, DFBDF,
+    NordsieckBDF, DNordsieckBDF
 
 end

lib/OrdinaryDiffEqBDF/src/alg_utils.jl

@@ -47,3 +47,19 @@ alg_order(alg::DFBDF) = 1 #dummy value
 isfsal(alg::DImplicitEuler) = false
 
 has_stiff_interpolation(::Union{QNDF, FBDF, DFBDF}) = true
+
+############################################ NordsieckBDF / DNordsieckBDF
+alg_order(alg::NordsieckBDF) = 1  # dummy: the running order lives in the cache
+alg_order(alg::DNordsieckBDF) = 1
+isadaptive(alg::DNordsieckBDF) = true
+get_current_alg_order(alg::NordsieckBDFAlgs, cache) = cache.order
+get_current_adaptive_order(alg::NordsieckBDFAlgs, cache) = cache.order
+has_stiff_interpolation(::NordsieckBDFAlgs) = true
+
+# The Newton increment norm is scaled by tq[2], which converts it into the units of
+# the local error test. That makes `NLNewton(κ = ...)` mean CVODE's NLSCOEF: the
+# fraction of the error-test budget the corrector may consume.
+has_special_newton_error(alg::NordsieckBDFAlgs) = true
+
+
+# The step-size logic is CVODE's (`cvSetEta` keeps h unless eta >= 1.5), so the

lib/OrdinaryDiffEqBDF/src/algorithms.jl

@@ -764,3 +764,141 @@ function DFBDF(;
 end
 
 @truncate_stacktrace DFBDF
+
+############################################ NordsieckBDF / DNordsieckBDF
+# ================================================================= algorithms
+@doc generic_solver_docstring(
+    """
+    An adaptive-order, adaptive-time BDF method on a propagated **Nordsieck**
+    history array `zn[j] = h^j/j! * y^(j)(t_n)`, following SUNDIALS CVODE.
+
+    Where `FBDF` stores raw `(t_i, u_i)` history and rebuilds the predictor and the
+    error/order estimates from it every step, this method propagates one array:
+    predicting is a Pascal-triangle shift, changing the step size is
+    `zn[j] *= eta^j`, and accepting a step is a rank-1 update. Because nothing is
+    reconstructed from stored points, a loose nonlinear solve cannot be amplified
+    into a step-size collapse, so the corrector can be solved to a fraction of the
+    local error budget (`nlsolve = NLNewton(κ = …)`, CVODE's NLSCOEF) instead of to
+    full accuracy. On stiff benchmarks that is worth roughly a factor of two in
+    f-evaluations relative to `FBDF`.
+
+    Dense output is the Nordsieck polynomial itself and is free.
+    """,
+    "NordsieckBDF",
+    "Multistep Method",
+    """
+    @article{byrne1975polyalgorithm,
+    title={A polyalgorithm for the numerical solution of ordinary differential equations},
+    author={Byrne, George D and Hindmarsh, Alan C},
+    journal={ACM Transactions on Mathematical Software},
+    volume={1}, number={1}, pages={71--96}, year={1975}}
+    @article{hindmarsh2005sundials,
+    title={{SUNDIALS}: Suite of nonlinear and differential/algebraic equation solvers},
+    author={Hindmarsh, Alan C and Brown, Peter N and Grant, Keith E and Lee, Steven L
+            and Serban, Radu and Shumaker, Dan E and Woodward, Carol S},
+    journal={ACM Transactions on Mathematical Software},
+    volume={31}, number={3}, pages={363--396}, year={2005}}""",
+    """
+    - `nlsolve`: nonlinear solver for the implicit stage. Its `κ` acts as CVODE's
+      NLSCOEF, i.e. the fraction of the local error budget the corrector is allowed
+      to consume, because the increment norm is scaled by the test quantity `tq[2]`.
+    - `max_order`: maximum BDF order (1–5).
+    - `step_limiter!`: function of the form `limiter!(u, integrator, p, t)`.
+    """,
+    """
+    nlsolve = NLNewton(),
+    extrapolant = :linear,
+    max_order::Val{MO} = Val{5}(),
+    step_limiter! = trivial_limiter!,
+    """
+)
+struct NordsieckBDF{MO, AD, F, F2, T, StepLimiter, CJ, QT} <:
+    OrdinaryDiffEqNewtonAdaptiveAlgorithm
+    max_order::Val{MO}
+    linsolve::F
+    nlsolve::F2
+    tol::T
+    extrapolant::Symbol
+    step_limiter!::StepLimiter
+    autodiff::AD
+    concrete_jac::CJ
+    stald::Bool
+    qmax::QT
+    qsteady_min::QT
+    qsteady_max::QT
+end
+
+function NordsieckBDF(;
+        max_order::Val{MO} = Val{5}(),
+        autodiff = AutoForwardDiff(), concrete_jac = nothing,
+        linsolve = nothing, nlsolve = NLNewton(), tol = nothing,
+        extrapolant = :linear, step_limiter! = trivial_limiter!, stald = false,
+        qsteady_min = 1 // 1, qsteady_max = 1 // 1, qmax = 10 // 1
+    ) where {MO}
+    autodiff = _fixup_ad(autodiff)
+    return NordsieckBDF(
+        max_order, linsolve, nlsolve, tol, extrapolant, step_limiter!,
+        autodiff, _unwrap_val(concrete_jac), stald, qmax, qsteady_min, qsteady_max
+    )
+end
+
+@truncate_stacktrace NordsieckBDF
+
+@doc generic_solver_docstring(
+    """
+    Fully implicit DAE solver: the [`NordsieckBDF`](@ref) method applied to
+    `f(du, u, p, t) = 0`. The Nordsieck array supplies both the state predictor and
+    the derivative `du = zn[1]/h`, so the corrector solves
+    `f((zn[1] + l[1]*acor)/h, ypred + acor, p, t) = 0` with leading coefficient
+    `cj = l[1]/h` — the same role IDA's `cj` plays.
+    """,
+    "DNordsieckBDF",
+    "Fully Implicit Multistep Method",
+    """
+    @article{hindmarsh2005sundials,
+    title={{SUNDIALS}: Suite of nonlinear and differential/algebraic equation solvers},
+    author={Hindmarsh, Alan C and Brown, Peter N and Grant, Keith E and Lee, Steven L
+            and Serban, Radu and Shumaker, Dan E and Woodward, Carol S},
+    journal={ACM Transactions on Mathematical Software},
+    volume={31}, number={3}, pages={363--396}, year={2005}}""",
+    """
+    - `nlsolve`: nonlinear solver for the implicit stage; its `κ` acts as NLSCOEF.
+    - `max_order`: maximum BDF order (1–5).
+    """,
+    """
+    nlsolve = NLNewton(),
+    extrapolant = :linear,
+    max_order::Val{MO} = Val{5}(),
+    """
+)
+struct DNordsieckBDF{MO, AD, F, F2, T, CJ, QT} <: DAEAlgorithm
+    max_order::Val{MO}
+    linsolve::F
+    nlsolve::F2
+    tol::T
+    extrapolant::Symbol
+    autodiff::AD
+    concrete_jac::CJ
+    stald::Bool
+    qmax::QT
+    qsteady_min::QT
+    qsteady_max::QT
+end
+
+function DNordsieckBDF(;
+        max_order::Val{MO} = Val{5}(),
+        autodiff = AutoForwardDiff(), concrete_jac = nothing,
+        linsolve = nothing, nlsolve = NLNewton(), tol = nothing,
+        extrapolant = :linear, stald = false,
+        qsteady_min = 1 // 1, qsteady_max = 1 // 1, qmax = 10 // 1
+    ) where {MO}
+    autodiff = _fixup_ad(autodiff)
+    return DNordsieckBDF(
+        max_order, linsolve, nlsolve, tol, extrapolant, autodiff,
+        _unwrap_val(concrete_jac), stald, qmax, qsteady_min, qsteady_max
+    )
+end
+
+@truncate_stacktrace DNordsieckBDF
+
+const NordsieckBDFAlgs = Union{NordsieckBDF, DNordsieckBDF}

lib/OrdinaryDiffEqBDF/src/bdf_caches.jl

@@ -800,3 +800,124 @@ function alg_cache(
         iters_from_event, dense, alg.step_limiter!, fd_weights, stald
     )
 end
+
+############################################ NordsieckBDF
+# ================================================================= caches
+@cache mutable struct NordsieckBDFCache{
+        MO, N, rateType, uNoUnitsType, uType, tType, coeffType, staldType, StepLimiter,
+    } <: BDFMutableCache
+    fsalfirst::rateType
+    nlsolver::N
+    zn::Vector{uType}
+    ypred::uType
+    acor::uType
+    tempv::uType
+    tmp::uType
+    atmp::uNoUnitsType
+    l::coeffType
+    tau::coeffType
+    tq::coeffType
+    order::Int
+    qprime::Int
+    qwait::Int
+    nst::Int
+    nef::Int
+    ncf::Int
+    index_acor::Int
+    max_order::Val{MO}
+    max_order_int::Int
+    hscale::tType
+    eta::tType
+    etamax::tType
+    etaq::tType
+    etaqm1::tType
+    etaqp1::tType
+    saved_tq5::tType
+    predicted::Bool
+    stald::staldType
+    step_limiter!::StepLimiter
+end
+
+@truncate_stacktrace NordsieckBDFCache 1
+
+mutable struct NordsieckBDFConstantCache{MO, N, uType, tType, coeffType, staldType} <:
+    OrdinaryDiffEqConstantCache
+    nlsolver::N
+    zn::Vector{uType}
+    ypred::uType
+    acor::uType
+    l::coeffType
+    tau::coeffType
+    tq::coeffType
+    order::Int
+    qprime::Int
+    qwait::Int
+    nst::Int
+    nef::Int
+    ncf::Int
+    index_acor::Int
+    max_order::Val{MO}
+    max_order_int::Int
+    hscale::tType
+    eta::tType
+    etamax::tType
+    etaq::tType
+    etaqm1::tType
+    etaqp1::tType
+    saved_tq5::tType
+    predicted::Bool
+    stald::staldType
+end
+
+# DAE caches carry `u₀` because `get_dae_uprev` uses it as the predictor the
+# correction `z` is measured against.
+
+function alg_cache(
+        alg::NordsieckBDF{MO}, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+        dt, reltol, p, calck, ::Val{true}, verbose
+    ) where {MO, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    γ, c = one(tTypeNoUnits), one(tTypeNoUnits)
+    nlsolver = build_nlsolver(
+        alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits,
+        uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(true), verbose
+    )
+    zn = [zero(u) for _ in 1:(MO + 1)]
+    coeffs() = zeros(typeof(t), MO + 3)
+    tq = zeros(typeof(t), 6)
+    stald = StabilityLimitDetectionState(real(uBottomEltypeNoUnits); enabled = alg.stald)
+    return NordsieckBDFCache{
+        MO, typeof(nlsolver), typeof(rate_prototype),
+        typeof(similar(u, uEltypeNoUnits)), typeof(u), typeof(t),
+        typeof(coeffs()), typeof(stald), typeof(alg.step_limiter!),
+    }(
+        zero(rate_prototype), nlsolver, zn, zero(u), zero(u), zero(u), zero(u),
+        similar(u, uEltypeNoUnits), coeffs(), coeffs(), tq,
+        1, 1, 2, 0, 0, 0, MO, Val(MO), MO,
+        zero(t), one(t), typeof(t)(NORD_ETA_MAX_FS), one(t), zero(t), zero(t),
+        zero(t), false, stald, alg.step_limiter!
+    )
+end
+
+function alg_cache(
+        alg::NordsieckBDF{MO}, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+        dt, reltol, p, calck, ::Val{false}, verbose
+    ) where {MO, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    γ, c = one(tTypeNoUnits), one(tTypeNoUnits)
+    nlsolver = build_nlsolver(
+        alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits,
+        uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(false), verbose
+    )
+    zn = [zero(u) for _ in 1:(MO + 1)]
+    coeffs() = zeros(typeof(t), MO + 3)
+    stald = StabilityLimitDetectionState(real(uBottomEltypeNoUnits); enabled = alg.stald)
+    return NordsieckBDFConstantCache{
+        MO, typeof(nlsolver), typeof(u), typeof(t), typeof(coeffs()), typeof(stald),
+    }(
+        nlsolver, zn, zero(u), zero(u), coeffs(), coeffs(), zeros(typeof(t), 6),
+        1, 1, 2, 0, 0, 0, MO, Val(MO), MO,
+        zero(t), one(t), typeof(t)(NORD_ETA_MAX_FS), one(t), zero(t), zero(t),
+        zero(t), false, stald
+    )
+end