search

SciML / ModelingToolkit.jl

An acausal modeling framework for automatically parallelized scientific machine learning (SciML) in Julia. A computer algebra system for integrated symbolics for physics-informed machine learning and automated transformations of differential equations
https://mtk.sciml.ai/dev/
Other
1.41k stars 203 forks source link

Events not working with automatic solver selection #2741

Open TorkelE opened 3 months ago

TorkelE commented 3 months ago

MWE:

using ModelingToolkit, OrdinaryDiffEq
@variables t X(t)
@parameters p d
D = Differential(t)

eqs = [D(X) ~ p - d*X]
discrete_events = [[1.0] => [p ~ 2p]]
@mtkbuild osys = ODESystem(eqs, t; discrete_events)

oprob = ODEProblem(osys, [X => 1.0], (0.0, 5.0), [p => 2.0, d => 1.0])
sol = solve(oprob)

yields a 


ERROR: No matching function wrapper was found!
Stacktrace:
  [1] _call(::Tuple{}, arg::Tuple{…}, fww::FunctionWrappersWrappers.FunctionWrappersWrapper{…})
    @ FunctionWrappersWrappers ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:23
  [2] _call(fw::Tuple{…}, arg::Tuple{…}, fww::FunctionWrappersWrappers.FunctionWrappersWrapper{…})
    @ FunctionWrappersWrappers ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:13
  [3] _call(fw::Tuple{…}, arg::Tuple{…}, fww::FunctionWrappersWrappers.FunctionWrappersWrapper{…})
    @ FunctionWrappersWrappers ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:13
  [4] _call(fw::Tuple{…}, arg::Tuple{…}, fww::FunctionWrappersWrappers.FunctionWrappersWrapper{…})
    @ FunctionWrappersWrappers ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:13
  [5] _call(fw::Tuple{…}, arg::Tuple{…}, fww::FunctionWrappersWrappers.FunctionWrappersWrapper{…})
    @ FunctionWrappersWrappers ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:13
  [6] (::FunctionWrappersWrappers.FunctionWrappersWrapper{…})(::Vector{…}, ::ModelingToolkit.MTKParameters{…}, ::Float64)
    @ FunctionWrappersWrappers ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:10
  [7] (::ODEFunction{…})(::Vector{…}, ::Vararg{…})
    @ SciMLBase ~/.julia/packages/SciMLBase/JUp1I/src/scimlfunctions.jl:2296
  [8] reset_fsal!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/k912u/src/integrators/integrator_utils.jl:493
  [9] apply_step!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/k912u/src/integrators/integrator_utils.jl:404
 [10] loopheader!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/k912u/src/integrators/integrator_utils.jl:14
 [11] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/k912u/src/solve.jl:549
 [12] #__solve#799
    @ ~/.julia/packages/OrdinaryDiffEq/k912u/src/solve.jl:7 [inlined]
 [13] __solve
    @ ~/.julia/packages/OrdinaryDiffEq/k912u/src/solve.jl:1 [inlined]
 [14] solve_call(_prob::ODEProblem{…}, args::CompositeAlgorithm{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:612
 [15] solve_call
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:569 [inlined]
 [16] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::ModelingToolkit.MTKParameters{…}, args::CompositeAlgorithm{…}; kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1080
 [17] solve_up
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1066 [inlined]
 [18] #solve#51
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1003 [inlined]
 [19] solve
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:993 [inlined]
 [20] #__solve#804
    @ ~/.julia/packages/OrdinaryDiffEq/k912u/src/solve.jl:542 [inlined]
 [21] __solve
    @ ~/.julia/packages/OrdinaryDiffEq/k912u/src/solve.jl:541 [inlined]
 [22] #__solve#72
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1394 [inlined]
 [23] __solve
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1386 [inlined]
 [24] #solve_call#44
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:612 [inlined]
 [25] solve_call
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:569 [inlined]
 [26] solve_up(::ODEProblem{…}, ::Nothing, ::Vector{…}, ::ModelingToolkit.MTKParameters{…}; kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1072
 [27] solve_up
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1066 [inlined]
 [28] #solve#51
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1003 [inlined]
 [29] solve(::ODEProblem{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:993
 [30] top-level scope
    @ ~/Desktop/Julia Playground/tmp/playground.jl:11
Some type information was truncated. Use `show(err)` to see complete types.

While for continuous:

continuous_events = [[X ~ 2.0] => [p ~ 2p]]
@mtkbuild osys = ODESystem(eqs, t; continuous_events)
oprob = ODEProblem(osys, [X => 1.0], (0.0, 5.0), [p => 2.0, d => 1.0])
sol = solve(oprob)

you get a 

ERROR: MethodError: no method matching get_tmp_cache(::OrdinaryDiffEq.ODEIntegrator{…}, ::CompositeAlgorithm{…}, ::OrdinaryDiffEq.DefaultCache{…})

Closest candidates are:
  get_tmp_cache(::Any, ::CompositeAlgorithm, ::OrdinaryDiffEq.CompositeCache)
   @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/k912u/src/integrators/integrator_interface.jl:179
  get_tmp_cache(::Any, ::OrdinaryDiffEqAlgorithm, ::OrdinaryDiffEq.OrdinaryDiffEqMutableCache)
   @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/k912u/src/integrators/integrator_interface.jl:125
  get_tmp_cache(::Any, ::OrdinaryDiffEqAlgorithm, ::OrdinaryDiffEq.OrdinaryDiffEqConstantCache)
   @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/k912u/src/integrators/integrator_interface.jl:118
  ...

Stacktrace:
  [1] get_tmp_cache
    @ ~/.julia/packages/OrdinaryDiffEq/k912u/src/integrators/integrator_interface.jl:111 [inlined]
  [2] (::ModelingToolkit.var"#456#464"{…})(u::Vector{…}, t::Float64, integ::OrdinaryDiffEq.ODEIntegrator{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/cAhZr/src/systems/callbacks.jl:486
  [3] determine_event_occurance
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/callbacks.jl:293 [inlined]
  [4] find_callback_time(integrator::OrdinaryDiffEq.ODEIntegrator{…}, callback::ContinuousCallback{…}, counter::Int64)
    @ DiffEqBase ~/.julia/packages/DiffEqBase/PBhFc/src/callbacks.jl:401
  [5] macro expansion
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/callbacks.jl:132 [inlined]
  [6] find_first_continuous_callback(integrator::OrdinaryDiffEq.ODEIntegrator{…}, callbacks::Tuple{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/PBhFc/src/callbacks.jl:127
  [7] find_first_continuous_callback
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/callbacks.jl:125 [inlined]
  [8] handle_callbacks!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/k912u/src/integrators/integrator_utils.jl:332
  [9] _loopfooter!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/k912u/src/integrators/integrator_utils.jl:254
 [10] loopfooter!
    @ ~/.julia/packages/OrdinaryDiffEq/k912u/src/integrators/integrator_utils.jl:207 [inlined]
 [11] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/k912u/src/solve.jl:554
 [12] #__solve#799
    @ ~/.julia/packages/OrdinaryDiffEq/k912u/src/solve.jl:7 [inlined]
 [13] __solve
    @ ~/.julia/packages/OrdinaryDiffEq/k912u/src/solve.jl:1 [inlined]
 [14] solve_call(_prob::ODEProblem{…}, args::CompositeAlgorithm{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:612
 [15] solve_call
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:569 [inlined]
 [16] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::ModelingToolkit.MTKParameters{…}, args::CompositeAlgorithm{…}; kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1080
 [17] solve_up
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1066 [inlined]
 [18] #solve#51
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1003 [inlined]
 [19] solve
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:993 [inlined]
 [20] #__solve#804
    @ ~/.julia/packages/OrdinaryDiffEq/k912u/src/solve.jl:542 [inlined]
 [21] __solve
    @ ~/.julia/packages/OrdinaryDiffEq/k912u/src/solve.jl:541 [inlined]
 [22] #__solve#72
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1394 [inlined]
 [23] __solve
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1386 [inlined]
 [24] #solve_call#44
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:612 [inlined]
 [25] solve_call
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:569 [inlined]
 [26] solve_up(::ODEProblem{…}, ::Nothing, ::Vector{…}, ::ModelingToolkit.MTKParameters{…}; kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1072
 [27] solve_up
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1066 [inlined]
 [28] #solve#51
    @ ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:1003 [inlined]
 [29] solve(::ODEProblem{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/PBhFc/src/solve.jl:993
 [30] top-level scope
    @ ~/Desktop/Julia Playground/tmp/playground.jl:17
Some type information was truncated. Use `show(err)` to see complete types.

Doing e.g. 

sol = solve(oprob, Tsit5())

and it works though.

Not selecting a solver is fine if there are no events:

using ModelingToolkit, OrdinaryDiffEq
@variables t X(t)
@parameters p d
D = Differential(t)

eqs = [D(X) ~ p - d*X]
@mtkbuild osys = ODESystem(eqs, t)

oprob = ODEProblem(osys, [X => 1.0], (0.0, 5.0), [p => 2.0, d => 1.0])
sol = solve(oprob)
ChrisRackauckas commented 3 months ago

@oscardssmith

oscardssmith commented 3 months ago

I believe this is fixed on version 6.80 of OrdinaryDiffEq. @TorkelE can you see if updating packages fixes this?

TorkelE commented 3 months ago

the mwe was carried out on v6.80 :(

oscardssmith commented 3 months ago

oh wait, @ChrisRackauckas we need to tag 6.80.1. We haven't tagged since merging https://github.com/SciML/OrdinaryDiffEq.jl/pull/2225