ModelingToolkit.jl icon indicating copy to clipboard operation
ModelingToolkit.jl copied to clipboard

Inifinite loop when solving an ODE with an unused NaN parameter

Open hersle opened this issue 1 year ago • 0 comments

Suppose I solve an ODE with a parameter p that I set to NaN, but where p does not "destroy" any unknowns in the system:

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using DifferentialEquations

@parameters p
@variables x(t) y(t)
@named sys = ODESystem([
    D(x) ~ 0.0
    y ~ x * p
], t)
sys = structural_simplify(sys)

prob = ODEProblem(sys, [sys.x => 0.0], (0.0, 1.0), [sys.p => NaN])
sol = solve(prob, Tsit5())

I know this probably looks stupid from this artificial example. But I would still argue that the ODE should solve successfully, and that sol[x] is "good output", while sol[y] is "bad output" full of NaNs. Instead, promote_u0() gets stuck in an infinite loop:

ERROR: LoadError: StackOverflowError:
Stacktrace:
  [1] poptask(W::Base.IntrusiveLinkedListSynchronized{Task})
    @ Base ./task.jl:999
  [2] wait()
    @ Base ./task.jl:1008
  [3] uv_write(s::Base.TTY, p::Ptr{UInt8}, n::UInt64)
    @ Base ./stream.jl:1048
  [4] unsafe_write(s::Base.TTY, p::Ptr{UInt8}, n::UInt64)
    @ Base ./stream.jl:1120
  [5] write
    @ ./strings/io.jl:248 [inlined]
  [6] print
    @ ./strings/io.jl:250 [inlined]
  [7] print(::Base.TTY, ::String, ::String, ::Vararg{String})
    @ Base ./strings/io.jl:46
  [8] println(::Base.TTY, ::String, ::Vararg{String})
    @ Base ./strings/io.jl:75
  [9] println(xs::String)
    @ Base ./coreio.jl:4
 [10] promote_u0(u0::Vector{Float64}, p::Vector{Float64}, t0::Float64)
    @ DiffEqBase ~/.julia/dev/DiffEqBase/src/forwarddiff.jl:355
 [11] promote_u0(u0::Vector{Float64}, p::Vector{Float64}, t0::Float64) (repeats 36884 times)
    @ DiffEqBase ~/.julia/dev/DiffEqBase/src/forwarddiff.jl:361
 [12] promote_u0 (repeats 2 times)
    @ ~/.julia/dev/DiffEqBase/src/forwarddiff.jl:361 [inlined]
 [13] get_concrete_problem(prob::ODEProblem{…}, isadapt::Bool; kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/dev/DiffEqBase/src/solve.jl:1171
 [14] get_concrete_problem
    @ ~/.julia/dev/DiffEqBase/src/solve.jl:1167 [inlined]
 [15] #solve_up#53
    @ ~/.julia/dev/DiffEqBase/src/solve.jl:1074 [inlined]
 [16] solve_up
    @ ~/.julia/dev/DiffEqBase/src/solve.jl:1066 [inlined]
 [17] #solve#51
    @ ~/.julia/dev/DiffEqBase/src/solve.jl:1003 [inlined]
 [18] solve(prob::ODEProblem{…}, args::Tsit5{…})
    @ DiffEqBase ~/.julia/dev/DiffEqBase/src/solve.jl:993
 [19] top-level scope
    @ ~/.julia/dev/SymBoltz/bug.jl:30
in expression starting at /home/hermasl/.julia/dev/SymBoltz/bug.jl:30
Some type information was truncated. Use `show(err)` to see complete types.

I am quite sure this worked before, but broke with a recent release.

hersle avatar Sep 04 '24 13:09 hersle