Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Callback TerminateSteadyState() does not work with AutoVern7/9 and Rodas4/Rodas5P #764

Closed
jo-ap opened this issue Oct 14, 2024 · 1 comment
Labels

Comments

@jo-ap
Copy link

jo-ap commented Oct 14, 2024

Describe the bug 🐞

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):

julia> import Pkg; Pkg.status()
Status `~/.julia/environments/v1.11/Project.toml`
  [31a5f54b] Debugger v0.7.10
  [0c46a032] DifferentialEquations v7.14.0
julia> versioninfo()
Julia Version 1.11.0
Commit 501a4f25c2b (2024-10-07 11:40 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 4 × Intel(R) Core(TM) i5-6300U CPU @ 2.40GHz
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, skylake)
Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)
@jo-ap jo-ap added the bug label Oct 14, 2024
@jo-ap
Copy link
Author

jo-ap commented Oct 15, 2024

I came here via the docs and just now noticed that it is only the repo for the documentation. Closed it here and reported it again: SciML/DifferentialEquations.jl#1058.

@jo-ap jo-ap closed this as completed Oct 15, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

1 participant