DifferentialEquations.jl change_t_via_interpolation! without save_everystep for large, expensive ODE

Hi all,

I am writing a domain-specific simulator built on top of DifferentialEquations.jl; I am solving a mathematically simple, but large / expensive system of ODEs (>10,000), and wish to minimise memory usage by setting save_everystep=false and periodically exporting some data to disk.

With no exporting to disk and using the SciML integrator interface, a simulation can finish in 27 function evaluations (as shown by integrator.stats); however, adding the periodic saving increases that to 603 function evaluations, both using DiffEqCallbacks or explicitly step!(integrator, next_export_time, true). I also set u_modified!(integrator, false) to avoid unnecessary saves as there are no discontinuities added.

Would it make sense to do each step!(integrator) - i.e. take the maximal successful time step - then use change_t_via_interpolation! to “go back” to each data-exporting timepoint between integrator.tprev and integrator.t? Do the ODE higher-order interpolants apply when save_everystep=false, but just between t and tprev?

I implemented this “maximal stepping, then back-interpolating” (though it’s quite a bit of code; I’m hoping this post is enough to spot any problems with the general approach, but if an MWE is needed please let me know) - however, the calls to change_t_via_interpolation! seem to destroy the integration accuracy in future timesteps. In other words, incrementally interpolating between tprev and t using change_t_via_interpolation! does not return the integrator to t. Is that expected, and is there a smarter way to let the integrator take maximal steps (and minimise function evaluations) while interpolating intra-step to export some data? I know interpolating can take some time, but I expect it to be much less than running more function evaluations.

Thanks,
Leonard

If you are using the interpolation to build the save values, there are 0 new function evaluations for any method that’s not a lazy interpolation (i.e. a Verner method).

No, just use integrator(t) and get it for free.

Yes, even if the post-solution interpolation is not turned on, the current solution interpolation is always built because it’s free.

It would increase the integration accuracy because it’s effectively just decreasing dt. If that difference is too large, then your tolerance is too large :sweat_smile:.

That’s not expected. An ODE that does this? Do you have an MWE?