Let me just start by saying, I’m not convinced this array-based approach is a good idea at all, and have been against it for about the past year now. I think we need to do something drastically differently and just remove this from the docs, and hopefully this post will explain why.
First of all,
What you did is not the recommended approach, precisely because of issues with discontinuities like this. What’s recommended is that you use a DiscreteCallback that updates the value. That will save two values at the discontinuity, one pre and one post discontinuity. The reason is because that makes sure the two exact (non-unique) values are known at that point, which ensures that you can always know the left-continuous and right-continuous behavior (assuming you keep the solution array ordered). So if you followed the documentation, you would not have this issue.
However, given that you seem very interested in this topic and possibly want to try and reduce memory, let me dig into why trying to assume that a point in time with a discontinuity has a unique value that can be saved, and also why DEDataArray is somewhat problematic.
The cache is just literally the cache. When you do f(du,u,p,t), what’s du? The cache variable. You’re just updating the values in the future cache variable. Of course, the user both interacts with the cache variables and doesn’t. The user never “explicitly” wants to interact with the cache, but we do hand it to the user as the du. But then when we do the broadcasts in the next step, those values need to be correct or the next u will not have the correct discrete values. That’s why they need to be updated.
It’s basically random due to the broadcast. Let me explain a bit. The interpolation is just calculated here:
Inside of that broadcast expression
@inbounds @.. out = y₀ + dt*(k[1]*b1Θ + k[2]*b2Θ + k[3]*b3Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ + k[7]*b7Θ)
you have DEDataArrays of both the forward and backward values. The array values are easy: you perform broadcast. But… what do you do with the discrete values?
https://github.com/SciML/DiffEqBase.jl/blob/master/src/data_array.jl#L148-L149
Let’s unpack that. Notice that there is no in-place broadcast overload, which means that the values of the DEDataArray are updated in-place but the discrete values are unchanged. The values are thus whatever was already put into the DEDataArray. Which comes from:
So technically this issue can be fixed by generating the in-place part with the right-hand side, but… that’s somewhat odd? I guess we can do that if all tests pass.
That description of broadcast gives a full scope of the issue with handling these objects. When you broadcast in-place, which is what all of the internal parts of the algorithm do, you can control what the value of the discrete variables are because they don’t change in broadcast, so you just need to set them correctly before hitting the algorithm which is what the cache loop is doing. However, for out-of-place broadcast, it’s going to just grab it from some array in the broadcast expression (at least broadcast always lowers the same way). However, if you dig in you’ll see:
y₀.ctrl_u = 0.4938506200402129
(k[1]).ctrl_u = 1.2099094768961067
at that point, so this issue of why you have to update the cache directly comes into play here.
So we could have it similar the right-most point, but that means that DEDataArray will work in DiffEq the way you’d like, which we can throw a test on and it should be okay-ish.
That seems a little snarky. I tried tens of explicit RKs and they seem to work. What don’t work are methods that use linear algebra because DEDataArray doesn’t have BLAS overloads. In theory doing linear algebra is just telling BLAS how to lower it, like https://github.com/SciML/LabelledArrays.jl/blob/master/src/larray.jl#L88 . If you’re willing to test it out, try adding that one line for DEDataArray and see if that fixes linear algebra. If it does, send a PR to DiffEqBase.
So…
In conclusion, if you donate some test cases we can probably get this fixed up for you with a really tricky fix as I described there. But again, trying to save one unique point at a discontinuity is just a problematic idea, and DEDataArray will always need a bit of manual intervention because what to do with the discrete variables when doing the array operations is ill-defined. And if you test out linear algebra a bit we can get that overload in there too. But I generally don’t think that DEDataArray is the best way to do explicit algebraic variables, but don’t have a much better idea right now.