DifferentialEquations.jl does not save fields of custom types after the callback function

[ronisbr]

Hi guys!

I am using DifferentialEquations.jl to simulate a satellite with continuous and discrete algorithms. @ChrisRackauckas helped me and the full discussion about it can be found here:

https://discourse.julialang.org/t/get-values-of-symbols-inside-a-function/848

However, I had a problem because I needed to modify custom fields of my new type inside the callback function and use those values in the next integration step. It turns out that DifferentialEquation.jldoes not save those values as I am expecting. Consider the following simplified example:

https://gist.github.com/ronisbr/1cac82902f78cf33ad732df79a5bd369

The dynamic function and callbacks are:

function f(t,u,du)
    du[1] = -0.5*u[1] + u.f1
    du[2] = -0.5*u[2]
end

callback = @ode_callback begin
    if t > 5
        u.f1 = 1.5
    end
    @ode_savevalues
end

Hence, I was expecting that u[1] becomes different of u[2] for t > 5. However, it does not happen as we can see in the following figure:

@ChrisRackauckas pointed out that in order to make this work, I need modify the callback function as described next:

callback = @ode_callback begin
    if t > 5
        for c in cache
            c.f1 = 1.5
        end
    end
    @ode_savevalues
end

The full code can be seen here:

https://gist.github.com/ronisbr/2883ffd81abcc379a259f13faf2b92ea

Now, the result is:

It does not seem a bug, as @ChrisRackauckas said, but it really should be documented.

[ChrisRackauckas]

The plot reveals you hit a bug in the interpolation. Hold up and I'll have a solution soon.

[ChrisRackauckas]

I don't know if that plot if from a different problem (the time goes from 0 to 30 and the discontinuity is at a different spot?), but I ran your code and saw the same wobblyness in the plot. This problem occured when there's a value which was both in tstops and saveat. It had an easy fix and so master for OrdinaryDiffEq.jl should be good. I'll release that patch soon enough. Below are the graphs using the patch.

Another thing that I had to handle though is the fact that discontinuities in the derivative mess with the interpolation. The fix is to save both before and after applying a discontinuous change. I fixed up the code to do this here:

https://gist.github.com/ChrisRackauckas/fc2b439e3e35b22b145a836ab57d42ce

which gives the solution:

(with the patch applied, using saveat instead of tstops works and will lower the number of computations that have to happen). However, you see that I had to do something extra with fsal...

---

Event Handling

That makes me think I led you astray. I'm sorry! These are details that the event handling framework takes care of for you, so maybe I should have had you define an event instead! Here's the event-handling way.

For events, you define a function which is <=0 when an event occurs. Here, lets do a discrete event each value in tstop:

const tstop = [5]

function event_f(t,u) # Event when event_f(t,u) == 0
  t ∉ tstop
end

Then you say what you want. This is the discontinuous change you defined before:

function apply_event!(u,cache)
  for c in cache
      c.f1 = 1.5
  end
end

The @ode_event macro has a few optional arguments. The first two are of interest to this problem. We define the event callback like:

const rootfind_event_loc = false
const interp_points = 0

callback = @ode_callback begin
  @ode_event event_f apply_event! rootfind_event_loc interp_points
end

The first part @ode_event event_f apply_event! is clear: you say what the event is and what to do when an event happens. rootfind_event_loc is whether to "use rootfinding to hone in on the event". This makes sense for the ball example where you let adaptive timestepping go past the point where the ball hits the wall, and then this "pulls the time back" to right where the ball hits the wall. But sense we did tstops = tstop, we will hit the point exactly so there's no reason to do this. Lastly, interp_points is the number interpolated points to check for the event. For example, if you took a timestep from t=10 to t=20, 9 interpolated points would mean to check for an event at t=11, t=12, t=13, ..., but since we know we are hitting our event exactly, we set that to zero.

The code using events in its entirety can be found here:

https://gist.github.com/ChrisRackauckas/44bd64d39255fd2a31ed0731c66401a1

and produces the following plot:

I'll add this example to DiffEqTutorials.jl.

[ronisbr]

Hi @ChrisRackauckas,

Thank you very much! Don't need to apologize, it was my fault since I did not describe the problem deeply enough.

I will port the code to use those events. My only question is: inside the apply_event! function, how can I get access to the time variable? Is it possible?

[ChrisRackauckas]

Not yet, but that's in the plan:

SciML/Roadmap#10 (comment)

(Feel free to comment there on the kinds of things you need in events and callbacks!)

If you need it right now, I can help you get it going. It should be as simple as making apply_event! to take in t in the event macro definition, which I see you have worked through.

[ronisbr]

hi @ChrisRackauckas,

When I finish this simulator, I can make a PR. My workaround was very simple: I created a global variable t_ctr that is set inside event_f.

[ChrisRackauckas]

Fixed by the new callbacks + cache iterators.