Saving only part of the solution in DifferentialEquations.jl

General Usage
,
1

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?

2

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

https://docs.juliadiffeq.org/latest/features/callback_functions.html#Using-Callbacks-1

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.

3

Jea there is a saving callback
http://docs.juliadiffeq.org/latest/features/callback_library.html#SavingCallback-1
which is a easy to use wrapper for exactly this purpose.

EDIT:
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
4

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
5

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
6

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).

7

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.