Using callbacks degrades ComponentArray solutions to Arrays. #284
Comments
|
Hm.. I'll give this a look. Thanks. |
|
So it looks it happens with other array types too. You get the same thing if you use LabelledArrays: using Sundials
using DifferentialEquations
using LabelledArrays
using Parameters: @unpack
function f(out,du,u,p,t)
@unpack u_ = u
out[1] = - 0.04u_[1] + 1e4*u_[2]*u_[3] - du[1]
out[2] = + 0.04u_[1] - 3e7*u_[2]^2 - 1e4*u_[2]*u_[3] - du[2]
out[3] = u_[1] + u_[2] + u_[3] - 1.0
end
u₀ = @LArray [1.0, 0, 0] (u_=1:3,)
du₀ = @LArray rand(3) (u_=1:3,)
tspan = (0.0,100000.0)
differential_vars = [true,true,false]
prob = DAEProblem{true}(f,du₀,u₀,tspan,differential_vars=differential_vars)
test_condition(u, t, integrator) = isapprox(u[1], 0.01, atol=0.01)
test_cb = DiscreteCallback(test_condition, terminate!)
# test_cb = nothing
sol = solve(prob, IDA(), callback=test_cb, jac_prototype=ones(3,3))
# sol = solve(prob, DABDF2(), callback=test_cb) |
|
It's interesting that a ComponentArray worked in the first place! Sundials as a C++ code does not compose with Julia abstractions. ComponentArrays seems to have slipped by and worked because Sundials grabbed the function pointer and used that, since it has an underlying array, but when it shows up in the DiscreteCallback it needs to be reinterpreted for the user. We should probably find a nice general way to handle this. |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
I’ve been solving a
DAEProblemwithIDA, and use aComponentArray(from ComponentArrays.jl) to represent the solution. This works, until I add aDiscreteCallback. Then, theComponentArrayis degraded to an vanillaArray. AFAIK, this degradation breaks the intended type promises of the common solve interface, and is an undocumented issue.The only potential workarounds I see are:
Arrays. This would be annoying, especially since it seems to be specific toSundials.jl.ComponentArrayin the appropriate places. If anyone knows how to do this without causing fresh array allocations, please let me know.This MWE fails on unpacking the (dummy) component
u_fromu, becauseuis anArray. If you set the callback tonothing, it works (uremains aComponentArray). If you change the solver toDImplicitEulerorDABDF2, it also works. I have not tested other callback types.Code from adapted from the DAE Example in the SciML docs.
FYI @jonniedie
The text was updated successfully, but these errors were encountered: