Saving only part of the solution in DifferentialEquations.jl

I want to use DifferentialEquations.jl to solve a system of stochastic time dependent PDEs, with a complicated geometry (multiple extended variables occupying different geometries, 1D+2D or 2D+3D, in the same simulation).
Given the nature of the system, I plan on using a fixed step-size Integrator (Euler-Maruyama with a dt that satisfies the stability condition for a diffusion equation, which I have in my model) .
For the moment I am interested in testing the deterministic simulations (no noise), using a simple Euler integrator.
I have written down the simulation function in ODE form (simulation(du, u, p, t)).
So far, the simulation runs, but I would be interested in the following:

  • Save only part of the variable array u (for example only u[1:10], and not the whole u) with a specific sampling frequency.
    So far I use saveat = 1. / sampling_freq (I am interested in saving only part of u because the number of variables of the resulting ODE system is quite large and I need to integrate the system for a long time)

  • Additionally I would also like to save the status of the system by means of a helper function (a cost function, or an energy): E = g(u, t)
    This pseudo-variable does not have an associated differential equation, so it cannot enter the model as an explicit variable.
    I would like to save this value whenever the solution (or part of it) is saved.

The reason I need to save this variable during the simualtion and I can’t just calculate it on the final solution object as post processing is because I count on not saving part of the u array that I need to calculate E.

Is there a way I could customize the resulting solution object as I desire?

1 Like

I’m no expert but I believe you should be able to do all of the above with a discrete Callback.

Maybe there’s a better way but with this you should be able to make the callback function compute the values you need and push! them to some prepared arrays.

Jea there is a saving callback
which is a easy to use wrapper for exactly this purpose.

You can use it like

using DiffEqCallbacks, OrdinaryDiffEq, LinearAlgebra

saved_values = SavedValues(typeof(u[1:10]),typeof(g(u,0)))
cb = SavingCallback((u,t,integrator)->(u[1:10], 
                    g(u,t), saved_values; saveat= 1. / sampling_freq)
sol = solve(prob, callback=cb)

Actually you can when specifying a mass matrix (how ever then the solver is restricted to implicit methods)

1 Like

Those are both good solutions, and you can also use the save_idxs keyword argument to choose to save only select arguments. Which one is best depends on what you need.

1 Like

Thank you very much for the suggestions.
I applied the SavingCallback and the save_idxs keyword, and it is all working for now.

@MatFi Btw, in the callback, it does not accept a simple Float64 for saveat, but instead it wants the explicit range. So I end up using the following:

saved_values = SavedValues(Real, Real)
cb = SavingCallback((u, p, t, integrator) -> g(u, p, t), saved_values, saveat= 0.0:1. / sampling_freq:1.0)
sol = (prob, callback=cb, saveat = 1. / sampling_freq, save_idxs = 1:10)

1 Like

I am continuing to develop the simulation I described and I was wondering if I could do the following:

  1. Save u[idxs] with saveat = start:step:end
  2. Save u[:] at the end of the integration

I was thinking of saving u[idxs] with the SavingCallback, which stores it in the saved_values object, and save u[:] at the end via as an element of sol = solve(prob, callback = saving_cb, save_on = false).

Trying it I noticed that when I use save_on = false the saving callback will also save only initial and final state of u .

Additionally, the interpolation of intermediate solutions between time-steps is not possible when the solution is stored in the saved_values object. Although this is not an important feature, it would be nice to have (for the deterministic simulations).

Jep true. but IMHO this lies in the nature that the type of the elements contained in saved_values can be arbitrary and so I don’t think that there is a general way to do interpolation here.

I haven’t observed such a behavior, but possible. You can probably use a work around:
In stead of feeding u[idxs] to the saving callback you can try t < end ? u[idxs] : u . so that at the last step the u to be saved switches to the full range. Haven’t tested it but should work.