CVODEIntegrator does not work with StepsizeLimiter from DiffEqCallbacks

[bgroenks96]

MWE:

using OrdinaryDiffEq, DiffEqCallbacks, Sundials
using ComponentArrays, LinearAlgebra

const nknots = 10
const h = 1.0/(nknots+1)
x = range(0, step=h, length=nknots)
u0 = sin.(π*x)

@inline function f(du,u,p,t)
  du .= zero(eltype(u))
  u₃ = @view u[3:end]
  u₂ = @view u[2:end-1]
  u₁ = @view u[1:end-2]
  k = p[1]
  @. du[2:end-1] = k*((u₃ - 2*u₂ + u₁)/(h^2.0))
  return nothing
end

p = ComponentArray(k=0.42)
jac_proto = Tridiagonal(similar(u0,nknots-1), similar(u0), similar(u0, nknots-1))
prob = ODEProblem(ODEFunction(f,jac_prototype=jac_proto), u0, (0.0,1.0), p)
@time sol = solve(prob, CVODE_BDF(), saveat=0.1, callback=StepsizeLimiter((u,p,t) -> 0.1))

Error:

ERROR: type CVODEIntegrator has no field dtcache
Stacktrace:
  [1] getproperty
    @ ~/.julia/packages/Sundials/k9hc3/src/common_interface/integrator_utils.jl:184 [inlined]
  [2] StepsizeLimiter_initialize(cb::DiscreteCallback{DiffEqCallbacks.var"#24#25", DiffEqCallbacks.StepsizeLimiterAffect{var"#7#8", Float64, Rational{Int64}, Bool}, typeof(DiffEqCallbacks.StepsizeLimiter_initialize), typeof(SciMLBase.FINALIZE_DEFAULT)}, u::Vector{Float64}, t::Float64, integrator::Sundials.CVODEIntegrator{Vector{Float64}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(k = 1,)}}}, Sundials.Handle{Sundials.CVODEMem}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Nothing, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(k = 1,)}}}, ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Tridiagonal{Float64, Vector{Float64}}, Tridiagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, CVODE_BDF{:Newton, :Dense, Nothing, Nothing}, SciMLBase.LinearInterpolation{Vector{Float64}, Vector{Vector{Float64}}}, DiffEqBase.DEStats}, CVODE_BDF{:Newton, :Dense, Nothing, Nothing}, ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Tridiagonal{Float64, Vector{Float64}}, Tridiagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}}, Sundials.FunJac{ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Tridiagonal{Float64, Vector{Float64}}, Tridiagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}}, Nothing, Nothing, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(k = 1,)}}}, Nothing, Nothing, Vector{Float64}, Nothing, Nothing, Nothing}, Nothing, Sundials.DEOptions{DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Float64, Vector{Float64}, Nothing, CallbackSet{Tuple{}, Tuple{DiscreteCallback{DiffEqCallbacks.var"#24#25", DiffEqCallbacks.StepsizeLimiterAffect{var"#7#8", Float64, Rational{Int64}, Bool}, typeof(DiffEqCallbacks.StepsizeLimiter_initialize), typeof(SciMLBase.FINALIZE_DEFAULT)}}}, Float64, Float64, typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE)}, Vector{Float64}, Tuple{Int64}, Vector{Float64}, Sundials.LinSolHandle{Sundials.Dense}, Sundials.MatrixHandle{Sundials.DenseMatrix}, Nothing})
    @ DiffEqCallbacks ~/.julia/packages/DiffEqCallbacks/ci89q/src/stepsizelimiters.jl:25
  [3] initialize!(u::Vector{Float64}, t::Float64, integrator::Sundials.CVODEIntegrator{Vector{Float64}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(k = 1,)}}}, Sundials.Handle{Sundials.CVODEMem}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Nothing, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(k = 1,)}}}, ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Tridiagonal{Float64, Vector{Float64}}, Tridiagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, CVODE_BDF{:Newton, :Dense, Nothing, Nothing}, SciMLBase.LinearInterpolation{Vector{Float64}, Vector{Vector{Float64}}}, DiffEqBase.DEStats}, CVODE_BDF{:Newton, :Dense, Nothing, Nothing}, ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Tridiagonal{Float64, Vector{Float64}}, Tridiagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}}, Sundials.FunJac{ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Tridiagonal{Float64, Vector{Float64}}, Tridiagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}}, Nothing, Nothing, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(k = 1,)}}}, Nothing, Nothing, Vector{Float64}, Nothing, Nothing, Nothing}, Nothing, Sundials.DEOptions{DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Float64, Vector{Float64}, Nothing, CallbackSet{Tuple{}, Tuple{DiscreteCallback{DiffEqCallbacks.var"#24#25", DiffEqCallbacks.StepsizeLimiterAffect{var"#7#8", Float64, Rational{Int64}, Bool}, typeof(DiffEqCallbacks.StepsizeLimiter_initialize), typeof(SciMLBase.FINALIZE_DEFAULT)}}}, Float64, Float64, typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE)}, Vector{Float64}, Tuple{Int64}, Vector{Float64}, Sundials.LinSolHandle{Sundials.Dense}, Sundials.MatrixHandle{Sundials.DenseMatrix}, Nothing}, any_modified::Bool, c::DiscreteCallback{DiffEqCallbacks.var"#24#25", DiffEqCallbacks.StepsizeLimiterAffect{var"#7#8", Float64, Rational{Int64}, Bool}, typeof(DiffEqCallbacks.StepsizeLimiter_initialize), typeof(SciMLBase.FINALIZE_DEFAULT)})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/S7V8q/src/callbacks.jl:17
  [4] initialize!(cb::CallbackSet{Tuple{}, Tuple{DiscreteCallback{DiffEqCallbacks.var"#24#25", DiffEqCallbacks.StepsizeLimiterAffect{var"#7#8", Float64, Rational{Int64}, Bool}, typeof(DiffEqCallbacks.StepsizeLimiter_initialize), typeof(SciMLBase.FINALIZE_DEFAULT)}}}, u::Vector{Float64}, t::Float64, integrator::Sundials.CVODEIntegrator{Vector{Float64}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(k = 1,)}}}, Sundials.Handle{Sundials.CVODEMem}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Nothing, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(k = 1,)}}}, ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Tridiagonal{Float64, Vector{Float64}}, Tridiagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, CVODE_BDF{:Newton, :Dense, Nothing, Nothing}, SciMLBase.LinearInterpolation{Vector{Float64}, Vector{Vector{Float64}}}, DiffEqBase.DEStats}, CVODE_BDF{:Newton, :Dense, Nothing, Nothing}, ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Tridiagonal{Float64, Vector{Float64}}, Tridiagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}}, Sundials.FunJac{ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Tridiagonal{Float64, Vector{Float64}}, Tridiagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}}, Nothing, Nothing, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(k = 1,)}}}, Nothing, Nothing, Vector{Float64}, Nothing, Nothing, Nothing}, Nothing, Sundials.DEOptions{DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Float64, Vector{Float64}, Nothing, CallbackSet{Tuple{}, Tuple{DiscreteCallback{DiffEqCallbacks.var"#24#25", DiffEqCallbacks.StepsizeLimiterAffect{var"#7#8", Float64, Rational{Int64}, Bool}, typeof(DiffEqCallbacks.StepsizeLimiter_initialize), typeof(SciMLBase.FINALIZE_DEFAULT)}}}, Float64, Float64, typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE)}, Vector{Float64}, Tuple{Int64}, Vector{Float64}, Sundials.LinSolHandle{Sundials.Dense}, Sundials.MatrixHandle{Sundials.DenseMatrix}, Nothing})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/S7V8q/src/callbacks.jl:7
  [5] initialize_callbacks!(integrator::Sundials.CVODEIntegrator{Vector{Float64}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(k = 1,)}}}, Sundials.Handle{Sundials.CVODEMem}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Nothing, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(k = 1,)}}}, ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Tridiagonal{Float64, Vector{Float64}}, Tridiagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, CVODE_BDF{:Newton, :Dense, Nothing, Nothing}, SciMLBase.LinearInterpolation{Vector{Float64}, Vector{Vector{Float64}}}, DiffEqBase.DEStats}, CVODE_BDF{:Newton, :Dense, Nothing, Nothing}, ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Tridiagonal{Float64, Vector{Float64}}, Tridiagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}}, Sundials.FunJac{ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Tridiagonal{Float64, Vector{Float64}}, Tridiagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}}, Nothing, Nothing, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(k = 1,)}}}, Nothing, Nothing, Vector{Float64}, Nothing, Nothing, Nothing}, Nothing, Sundials.DEOptions{DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Float64, Vector{Float64}, Nothing, CallbackSet{Tuple{}, Tuple{DiscreteCallback{DiffEqCallbacks.var"#24#25", DiffEqCallbacks.StepsizeLimiterAffect{var"#7#8", Float64, Rational{Int64}, Bool}, typeof(DiffEqCallbacks.StepsizeLimiter_initialize), typeof(SciMLBase.FINALIZE_DEFAULT)}}}, Float64, Float64, typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE)}, Vector{Float64}, Tuple{Int64}, Vector{Float64}, Sundials.LinSolHandle{Sundials.Dense}, Sundials.MatrixHandle{Sundials.DenseMatrix}, Nothing}, initialize_save::Bool)
    @ Sundials ~/.julia/packages/Sundials/k9hc3/src/common_interface/solve.jl:1729
  [6] initialize_callbacks!(integrator::Sundials.CVODEIntegrator{Vector{Float64}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(k = 1,)}}}, Sundials.Handle{Sundials.CVODEMem}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Nothing, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(k = 1,)}}}, ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Tridiagonal{Float64, Vector{Float64}}, Tridiagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, CVODE_BDF{:Newton, :Dense, Nothing, Nothing}, SciMLBase.LinearInterpolation{Vector{Float64}, Vector{Vector{Float64}}}, DiffEqBase.DEStats}, CVODE_BDF{:Newton, :Dense, Nothing, Nothing}, ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Tridiagonal{Float64, Vector{Float64}}, Tridiagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}}, Sundials.FunJac{ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Tridiagonal{Float64, Vector{Float64}}, Tridiagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}}, Nothing, Nothing, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(k = 1,)}}}, Nothing, Nothing, Vector{Float64}, Nothing, Nothing, Nothing}, Nothing, Sundials.DEOptions{DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Float64, Vector{Float64}, Nothing, CallbackSet{Tuple{}, Tuple{DiscreteCallback{DiffEqCallbacks.var"#24#25", DiffEqCallbacks.StepsizeLimiterAffect{var"#7#8", Float64, Rational{Int64}, Bool}, typeof(DiffEqCallbacks.StepsizeLimiter_initialize), typeof(SciMLBase.FINALIZE_DEFAULT)}}}, Float64, Float64, typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE)}, Vector{Float64}, Tuple{Int64}, Vector{Float64}, Sundials.LinSolHandle{Sundials.Dense}, Sundials.MatrixHandle{Sundials.DenseMatrix}, Nothing})
    @ Sundials ~/.julia/packages/Sundials/k9hc3/src/common_interface/solve.jl:1724
  [7] __init(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(k = 1,)}}}, ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Tridiagonal{Float64, Vector{Float64}}, Tridiagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::CVODE_BDF{:Newton, :Dense, Nothing, Nothing}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any}; verbose::Bool, callback::DiscreteCallback{DiffEqCallbacks.var"#24#25", DiffEqCallbacks.StepsizeLimiterAffect{var"#7#8", Float64, Rational{Int64}, Bool}, typeof(DiffEqCallbacks.StepsizeLimiter_initialize), typeof(SciMLBase.FINALIZE_DEFAULT)}, abstol::Float64, reltol::Float64, saveat::Float64, tstops::Vector{Float64}, maxiters::Int64, dt::Nothing, dtmin::Float64, dtmax::Float64, timeseries_errors::Bool, dense_errors::Bool, save_everystep::Bool, save_idxs::Nothing, save_on::Bool, save_start::Bool, save_end::Bool, dense::Bool, progress::Bool, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), save_timeseries::Nothing, advance_to_tstop::Bool, stop_at_next_tstop::Bool, userdata::Nothing, alias_u0::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Sundials ~/.julia/packages/Sundials/k9hc3/src/common_interface/solve.jl:476
  [8] __solve(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(k = 1,)}}}, ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Tridiagonal{Float64, Vector{Float64}}, Tridiagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::CVODE_BDF{:Newton, :Dense, Nothing, Nothing}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any}, recompile::Type{Val{true}}; calculate_error::Bool, kwargs::Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:saveat, :callback), Tuple{Float64, DiscreteCallback{DiffEqCallbacks.var"#24#25", DiffEqCallbacks.StepsizeLimiterAffect{var"#7#8", Float64, Rational{Int64}, Bool}, typeof(DiffEqCallbacks.StepsizeLimiter_initialize), typeof(SciMLBase.FINALIZE_DEFAULT)}}}})
    @ Sundials ~/.julia/packages/Sundials/k9hc3/src/common_interface/solve.jl:14
  [9] #solve_call#39
    @ ~/.julia/packages/DiffEqBase/S7V8q/src/solve.jl:221 [inlined]
 [10] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(k = 1,)}}}, ODEFunction{true, typeof(f), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Tridiagonal{Float64, Vector{Float64}}, Tridiagonal{Float64, Vector{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::ComponentVector{Float64}, args::CVODE_BDF{:Newton, :Dense, Nothing, Nothing}; kwargs::Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:saveat, :callback), Tuple{Float64, DiscreteCallback{DiffEqCallbacks.var"#24#25", DiffEqCallbacks.StepsizeLimiterAffect{var"#7#8", Float64, Rational{Int64}, Bool}, typeof(DiffEqCallbacks.StepsizeLimiter_initialize), typeof(SciMLBase.FINALIZE_DEFAULT)}}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/S7V8q/src/solve.jl:248
 [11] #solve#40
    @ ~/.julia/packages/DiffEqBase/S7V8q/src/solve.jl:234 [inlined]
 [12] top-level scope
    @ timing.jl:210

[ChrisRackauckas]

Thanks for sharing. Sundials indeed doesn't have a dtmax cache or dtcache so its integrator would have to be expanded. It's not hard, but it seems like your example would be better done with FBDF anyways, since Sundials wouldn't specialize on Tridiagonal and so FBDF would certainly be faster here. If you find a case that really needs this then I can handle it though, but with limited time this is easy to fix but probably not a major thing (I assume you already tried FBDF here?)

[bgroenks96]

CVODE_BDF has jac_upper and jac_lower fields which can be used to take advantage of the tridiagonal structure. I have not tried FBDF on this MWE, but on my actual nonlinear heat equation use cases, CVODE_BDF is pretty much always the fastest option (though not always stable).

[ChrisRackauckas]

Can you share an MWE? FBDF has been performing really well since the start of 2022, so I'd be happy to hone in on any issue there first. That would be a bit harder to solve, but would have more benefits so I'd prioritize it.

[bgroenks96]

So testing now on my "real" model (no MWE immediately available), CVODE_BDF is about 3-4x faster in wallclock execution time than FBDF using default tolerances; I haven't checked the error (compared to a low tolerance baseline), but the plots look pretty much identical.

However, I have had trouble with CVODE_BDF failing to converge (MaxIters) for some simulations with certain parameter settings, which is why I typically fall back to SSPRK22 or even Euler, both of which still work fine. I tried out one of these, and FBDF converged fine and was faster than Euler, so that's nice.

[ChrisRackauckas]

If you could share that model somehow with Yingbo and I that would be helpful.

[bgroenks96]

It is on GitHub and is installable via the package manager. I also have example scripts, and I can prepare one to compare the two solvers, if that would be helpful.