Using TerminateSteadyState() as a callback to solve for an ODE problem with AutoVern7 or AutoVern9 and Rodas4 or Rodas5P (which is suggested by https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/) fails.
I traced the problem to allDerivPass, which calls DiffEqBase.get_du!. There, out is updated to integrator.fsallast, which is not set. Therefore, get_du! returns nothing for the RHS and this results in errors further down.
Expected behavior
The callback TerminateSteadyState() should work with AutoVern7 or AutoVern9 together with Rodas4 or Rodas5P. I guess this boils down to DiffEqBase.get_du! returning the RHS for the integrator correctly.
If the error is to be expected, this should be documented and a warning thrown when calling solve.
Minimal Reproducible Example 👇
using DifferentialEquations
# Rosenzweig-MacArthur Predator-Prey Model
function f!(du, u, p, t)
x, y = u
α, β, δ, γ, η, ρ = p
du[1] = ρ*x*(1-x) - α*(η + β)*x*y/(δ + γ*x)
du[2] = -η*y + γ*(η + β)*x*y/(δ + γ*x)
end
u0 = [0.1,0.1]
p = [0.1, 0.2, 1.5, 0.4, 0.2, 0.3]
prob = ODEProblem(f!, u0, (0.0, 1000.0), p)
# error comes from integrator.fsallast not being set when calling DiffEqBase.get_du! in allDerivPass
# occurs for all combinations of AutoVern7/AutoVern9 and Rodas4/Rodas5P
solve(prob, AutoVern9(Rodas4()); callback = TerminateSteadyState())
Error & Stacktrace ⚠️
julia> solve(prob, AutoVern9(Rodas4()); callback = TerminateSteadyState())
ERROR: MethodError: Cannot `convert` an object of type Nothing to an object of type Float64
The function `convert` exists, but no method is defined for this combination of argument types.
Closest candidates are:
convert(::Type{T}, ::VectorizationBase.AbstractSIMD) where T<:Union{Bool, Float16, Float32, Float64, Int16, Int32, Int64, Int8, UInt16, UInt32, UInt64, UInt8, SIMDTypes.Bit}
@ VectorizationBase ~/.julia/packages/VectorizationBase/LqJbS/src/base_defs.jl:201
convert(::Type{<:Real}, ::T) where T<:SparseConnectivityTracer.AbstractTracer
@ SparseConnectivityTracer ~/.julia/packages/SparseConnectivityTracer/hOFV1/src/overloads/conversion.jl:10
convert(::Type{T}, ::T) where T
@ Base Base.jl:126
...
Stacktrace:
[1] fill!(dest::Vector{Float64}, x::Nothing)
@ Base ./array.jl:327
[2] copyto!
@ ./broadcast.jl:928 [inlined]
[3] materialize!
@ ./broadcast.jl:878 [inlined]
[4] materialize!
@ ./broadcast.jl:875 [inlined]
[5] get_du!
@ ~/.julia/packages/OrdinaryDiffEqCore/YWZVL/src/integrators/integrator_interface.jl:82 [inlined]
[6] allDerivPass(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, abstol::Float64, reltol::Float64, min_t::Nothing)
@ DiffEqCallbacks ~/.julia/packages/DiffEqCallbacks/n5zrr/src/terminatesteadystate.jl:7
[7] (::DiffEqCallbacks.var"#104#106"{…})(u::Vector{…}, t::Float64, integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
@ DiffEqCallbacks ~/.julia/packages/DiffEqCallbacks/n5zrr/src/terminatesteadystate.jl:68
[8] apply_discrete_callback!
@ ~/.julia/packages/DiffEqBase/ODi5x/src/callbacks.jl:611 [inlined]
[9] handle_callbacks!
@ ~/.julia/packages/OrdinaryDiffEqCore/YWZVL/src/integrators/integrator_utils.jl:355 [inlined]
[10] _loopfooter!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
@ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/YWZVL/src/integrators/integrator_utils.jl:243
[11] loopfooter!
@ ~/.julia/packages/OrdinaryDiffEqCore/YWZVL/src/integrators/integrator_utils.jl:207 [inlined]
[12] solve!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
@ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/YWZVL/src/solve.jl:552
[13] #__solve#75
@ ~/.julia/packages/OrdinaryDiffEqCore/YWZVL/src/solve.jl:7 [inlined]
[14] __solve
@ ~/.julia/packages/OrdinaryDiffEqCore/YWZVL/src/solve.jl:1 [inlined]
[15] #solve_call#44
@ ~/.julia/packages/DiffEqBase/ODi5x/src/solve.jl:612 [inlined]
[16] solve_call
@ ~/.julia/packages/DiffEqBase/ODi5x/src/solve.jl:569 [inlined]
[17] #solve_up#53
@ ~/.julia/packages/DiffEqBase/ODi5x/src/solve.jl:1092 [inlined]
[18] solve_up
@ ~/.julia/packages/DiffEqBase/ODi5x/src/solve.jl:1078 [inlined]
[19] #solve#51
@ ~/.julia/packages/DiffEqBase/ODi5x/src/solve.jl:1015 [inlined]
[20] top-level scope
@ REPL[2]:1
Some type information was truncated. Use `show(err)` to see complete types.
Environment (please complete the following information):
Describe the bug 🐞
Using
TerminateSteadyState()
as a callback tosolve
for an ODE problem withAutoVern7
orAutoVern9
andRodas4
orRodas5P
(which is suggested by https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/) fails. I traced the problem toallDerivPass
, which callsDiffEqBase.get_du!
. There,out
is updated tointegrator.fsallast
, which is not set. Therefore,get_du!
returnsnothing
for the RHS and this results in errors further down.Expected behavior
The callback
TerminateSteadyState()
should work withAutoVern7
orAutoVern9
together withRodas4
orRodas5P
. I guess this boils down toDiffEqBase.get_du!
returning the RHS for the integrator correctly. If the error is to be expected, this should be documented and a warning thrown when callingsolve
.Minimal Reproducible Example 👇
Error & Stacktrace ⚠️
Environment (please complete the following information):