search

SciML / Sundials.jl

Julia interface to Sundials, including a nonlinear solver (KINSOL), ODE's (CVODE and ARKODE), and DAE's (IDA) in a SciML scientific machine learning enabled manner
https://diffeq.sciml.ai
BSD 2-Clause "Simplified" License
209 stars 77 forks source link

Error with SavingCallback defaulting to Array{Float64,1} with CVODE? #275

Open HelgavonLichtenstein opened 4 years ago

HelgavonLichtenstein commented 4 years ago

I'm not very good at making a MWE so I took the example from and replaced Tsit5() with CVODE_BDF(). I assumed it would still work since it's just a different solver, albeit not very efficient since it is a stiff solver. 

prob = ODEProblem((du,u,p,t) -> du .= u, rand(4,4), (0.0,1.0))
saved_values = SavedValues(Float64, Tuple{Float64,Float64})
cb = SavingCallback((u,t,integrator)->(tr(u),norm(u)), saved_values)
sol = solve(prob, CVODE_BDF(), callback=cb)

It errors with 

Closest candidates are:
  tr(!Matched::Array{T,2}) where T at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\LinearAlgebra\src\dense.jl:331
  tr(!Matched::Hermitian) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\LinearAlgebra\src\symmetric.jl:363
  tr(!Matched::Diagonal) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\LinearAlgebra\src\diagonal.jl:547
  ...
in top-level scope at untitled:6
in  at DiffEqBase\ytJuW\src\solve.jl:100
in #solve#458 at DiffEqBase\ytJuW\src\solve.jl:102 
in solve_up##kw at DiffEqBase\ytJuW\src\solve.jl:107 
in #solve_up#459 at DiffEqBase\ytJuW\src\solve.jl:114 
in solve_call##kw at DiffEqBase\ytJuW\src\solve.jl:65 
in #solve_call#455 at DiffEqBase\ytJuW\src\solve.jl:92
in  at Sundials\UNvwn\src\common_interface\solve.jl:10
in __solve##kw at Sundials\UNvwn\src\common_interface\solve.jl:10 
in __solve##kw at Sundials\UNvwn\src\common_interface\solve.jl:10 
in __solve##kw at Sundials\UNvwn\src\common_interface\solve.jl:10 
in __solve##kw at Sundials\UNvwn\src\common_interface\solve.jl:10 
in #__solve#22 at Sundials\UNvwn\src\common_interface\solve.jl:10
in  at Sundials\UNvwn\src\common_interface\solve.jl:42
in #__init#23 at Sundials\UNvwn\src\common_interface\solve.jl:322
in initialize_callbacks! at Sundials\UNvwn\src\common_interface\solve.jl:1318
in initialize_callbacks! at Sundials\UNvwn\src\common_interface\solve.jl:1323
in initialize! at DiffEqBase\ytJuW\src\callbacks.jl:315 
in initialize! at DiffEqBase\ytJuW\src\callbacks.jl:325 
in saving_initialize at DiffEqCallbacks\yz7kf\src\saving.jl:91 
in SavingAffect at DiffEqCallbacks\yz7kf\src\saving.jl:46 
in  at DiffEqCallbacks\yz7kf\src\saving.jl:77
in  at untitled:5

My code creates this error, which is almost identical apart from in #solve_up#459 at DiffEqBase\ytJuW\src\solve.jl:114 from example above in #solve_up#459 at DiffEqBase\ytJuW\src\solve.jl:110 in error from below

MethodError: Cannot `convert` an object of type Array{Float64,1} to an object of type Float64
Closest candidates are:
  convert(::Type{T}, !Matched::Unitful.Quantity) where T<:Real at C:\Users\Helga\.julia\packages\Unitful\MOEUx\src\conversion.jl:145
  convert(::Type{T}, !Matched::Unitful.Level) where T<:Real at C:\Users\Helga\.julia\packages\Unitful\MOEUx\src\logarithm.jl:22
  convert(::Type{T}, !Matched::Unitful.Gain) where T<:Real at C:\Users\Helga\.julia\packages\Unitful\MOEUx\src\logarithm.jl:62
  ...
in top-level scope at untitled:51
in  at DiffEqBase\ytJuW\src\solve.jl:100
in #solve#458 at DiffEqBase\ytJuW\src\solve.jl:102 
in solve_up##kw at DiffEqBase\ytJuW\src\solve.jl:107 
in #solve_up#459 at DiffEqBase\ytJuW\src\solve.jl:110
in solve_call##kw at DiffEqBase\ytJuW\src\solve.jl:65 
in #solve_call#455 at DiffEqBase\ytJuW\src\solve.jl:92
in  at Sundials\UNvwn\src\common_interface\solve.jl:10
in __solve##kw at Sundials\UNvwn\src\common_interface\solve.jl:10 
in __solve##kw at Sundials\UNvwn\src\common_interface\solve.jl:10 
in __solve##kw at Sundials\UNvwn\src\common_interface\solve.jl:10 
in __solve##kw at Sundials\UNvwn\src\common_interface\solve.jl:10 
in #__solve#22 at Sundials\UNvwn\src\common_interface\solve.jl:10
in  at Sundials\UNvwn\src\common_interface\solve.jl:42
in #__init#23 at Sundials\UNvwn\src\common_interface\solve.jl:322
in initialize_callbacks! at Sundials\UNvwn\src\common_interface\solve.jl:1318
in initialize_callbacks! at Sundials\UNvwn\src\common_interface\solve.jl:1323
in initialize! at DiffEqBase\ytJuW\src\callbacks.jl:315 
in initialize! at DiffEqBase\ytJuW\src\callbacks.jl:321 
in initialize! at DiffEqBase\ytJuW\src\callbacks.jl:325 
in saving_initialize at DiffEqCallbacks\yz7kf\src\saving.jl:91 
in SavingAffect at DiffEqCallbacks\yz7kf\src\saving.jl:46 
in  at DiffEqCallbacks\yz7kf\src\saving.jl:77
in copyat_or_push! at RecursiveArrayTools\DBKJE\src\utils.jl
in copyto! at base\multidimensional.jl:962
in setindex! at base\array.jl:826
ChrisRackauckas commented 4 years ago 
using Sundials, DiffEqCallbacks, LinearAlgebra
prob = ODEProblem((du,u,p,t) -> du .= u, rand(4,4), (0.0,1.0))
saved_values = SavedValues(Float64, Tuple{Float64,Float64})
cb = SavingCallback((u,t,integrator)->(tr(reshape(u,4,4)),norm(reshape(u,4,4))), saved_values)
sol = solve(prob, CVODE_BDF(), callback=cb)

is a workaround. We have to internally flatten for Sundials' API, and we normally reshape before f to hide that from the user. But it looks like it's exposing itself here. This one is one of those issues that is particularly difficult to solve, but extremely easy to workaround, so I'm just going to leave it hanging until after JuliaCon.