get_du cause accessing undefined reference

[shipengcheng1230]

Hello, I am trying to retrieve the derivatives and save them to elsewhere and find a problem. Consider the following naive example:

using DifferentialEquations

function foo(du, u, p, t)
    du.x[1] .= u.x[1] / 10
    du.x[2] .= u.x[2] / 100
    du.x[3] .= u.x[3] / 1000
end
u0 = ArrayPartition(rand(2, 3), rand(5), rand(3, 2))
tspan = (0.0, 1000.0)
prob = ODEProblem(foo, u0, tspan)
cb1 = FunctionCallingCallback((u, t, integrator) -> integrator.fsallast)
cb2 = FunctionCallingCallback((u, t, integrator) -> integrator.cache.tmp)
cb3 = FunctionCallingCallback((u, t, integrator) -> get_du(integrator))

solve(prob, Tsit5(); callback=cb1) # ERROR: UndefRefError: access to undefined reference
solve(prob, Tsit5(); callback=cb2) # OK
solve(prob, Tsit5(); callback=cb3) # ERROR: UndefRefError: access to undefined reference

Accessing field fsallast of integrator would trigger this error. However, the following seems works without a glitch:

integrator = init(prob, Tsit5())
for i = 1: 100
    step!(integrator)
    x = get_du(integrator)
    y = integrator.cache.tmp
    @show x.x[1][1, 1], y.x[1][1, 1]
end # seems no error

But the derivatives are different (which one should I use?). How do I reliably retrieve the derivatives in such case? I am all on the latest Julia and packages.

Thanks!

[ChrisRackauckas]

I'd recommend integrator(t,Val{1}) for now, but we should definitely do a getproperty overload for du.

[shipengcheng1230]

I tried your solution and for the given example above:

cb4 = FunctionCallingCallback((u, t, integrator) -> integrator(integrator.t, Val{1}))
solve(prob, Tsit5(); callback=cb4) # perfect

integrator = init(prob, Tsit5())
for i = 1: 10
    step!(integrator)
    x = get_du(integrator).x[1][1]
    y = integrator.cache.tmp.x[1][1]
    z = integrator(integrator.t, Val{1}).x[1][1]
    @show x, y, z
end

# the ouputs
(x, y, z) = (0.032038298097310874, 0.3203829488581116, 0.03203829809731048)
(x, y, z) = (0.03397110155476714, 0.3397209584709375, 0.033971101554766715)
(x, y, z) = (0.0610254155996639, 0.6226214991916795, 0.06102541559966307)
(x, y, z) = (0.1365346808286418, 1.4179867554173868, 0.13653468082863987)
(x, y, z) = (0.3133795300172519, 3.261577991818207, 0.3133795300172463)
(x, y, z) = (0.8675610917180112, 9.164868397121804, 0.8675610917179984)
(x, y, z) = (2.4280376159197767, 25.66997888206144, 2.428037615919732)
(x, y, z) = (7.758461565274317, 82.73784855736018, 7.75846156527423)
(x, y, z) = (25.0530816002286, 267.3307717738807, 25.05308160022844)
(x, y, z) = (87.85608609552757, 941.3764940703618, 87.85608609552658)

There is slight difference between get_du and integrator(t,Val{1}). BTW, I wasn't able to find the corresponding code that has integrator(t, ::Val{N}) signature. Could you please give me a hint? Thanks!

[ChrisRackauckas]

BTW, I wasn't able to find the corresponding code that has integrator(t, ::Val{N}) signature. Could you please give me a hint? Thanks!

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/integrator_utils.jl#L374

which goes into the dense output handling

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/dense/generic_dense.jl

which then uses the analytical derivative of the interpolant, or just the point value if at an endpoint (with a bit of autodiff handling). That would be the numerical difference.

integrator.cache.tmp isn't guaranteed to have any value. It's a cache that can be written to, and so other things will overwrite it as they please. But you can use it as a cache as well. That explains the behavior there.

I think that settles the issue? Let me know if you have anything else.