Adaptive timestepping

[johnomotani]

Add several adaptive timestep methods based on [Fekete, Conde and Shadid, "Embedded pairs for optimal explicit strong stability preserving Runge-Kutta methods", Journal of Computational and Applied Mathematics 421 (2022) 114325, https://doi.org/10.1016/j.cam.2022.114325], as well as the Runge-Kutta-Fehlberg method.

Timestep limits are set by the lowest of: an estimate based on the truncation error estimated by the embedded pair of methods in the RK schemes; CFL limits from z-advection and 'vpa-advection' (i.e. parallel accelleration) terms in the kinetic equation; a maximum-increase-per-step factor. It is also possible to set a minimum timestep (which can help to get through initial transients).

Several diagnostics are written and can be plotted to help tune the adaptive timestep parameters.

For more details see the new section of the docs: https://mabarnes.github.io/moment_kinetics/previews/PR194/timestepping/#Timestepping

Also fixes several small bugs/problems in makie_post_processing.

[mrhardman]

Is the default time-stepping behaviour still intended to be the RK4 with a fixed timestep size, or do you think that these optimisations are robust enough to replace the current default behaviour? Do I understand correctly that your new optimised options are available by choosing the type input parameter?

[johnomotani]

I think we need to test out the adaptive timestepping methods on more cases before we think about making (one of) them the default. I've only tried out one case so far. Even once we are happy that it works for most cases, there are more parameters that need setting (not just dt!) so we would need to decide (based on experience) on some sensible defaults for them.

I'd suggest merging this PR now(ish), so we can experiment with the methods, and having a future PR to update default values for the settings, and possibly (if we decide it's a good idea) the 'type' option that sets the method used.

[mrhardman]

I am not against merging, but before we do I would want to discuss how this may or may not affect the readability of the RK update steps. I think that the changes in features-gyroaverages (#187) to the time stepping are also useful (isolating the moments and diagnostics calculation into a single location) and would potentially complement (but also potentially conflict with) any changes aimed at making the time stepper even more encapsulated. I haven't had a chance to understand what you have changed here, but I would like to!

[johnomotani]

I am not against merging, but before we do I would want to discuss how this may or may not affect the readability of the RK update steps. I think that the changes in features-gyroaverages (#187) to the time stepping are also useful (isolating the moments and diagnostics calculation into a single location) and would potentially complement (but also potentially conflict with) any changes aimed at making the time stepper even more encapsulated. I haven't had a chance to understand what you have changed here, but I would like to!

I think the changes here and in #187 are mostly non-overlapping, so hopefully both simplify the RK update a bit, and shouldn't be too hard to merge (although I've been too optimistic about merging more often than not!).

[mrhardman]

I think the changes here and in #187 are mostly non-overlapping, so hopefully both simplify the RK update a bit, and shouldn't be too hard to merge (although I've been too optimistic about merging more often than not!).

If you indicate when you think you are fairly converged, I'll merge this PR into a test branch off #187. #187 should only have further changes in moment_kinetics/src/gyroaverages.jl, since introducing distributed-memory MPI to the gyroaverage operator should only require changes in the gyroaverages.jl module.

[mrhardman]

Is it really true that we need Float128 types to have adaptive timestepping? Does having Float128 as a type slow down the speed of the calculation significantly?

I am referencing this comment in master, I assume it came from this PR? https://github.com/mabarnes/moment_kinetics/blob/1800f807308827a39f7e0bf0a3f5928a6990a0e1/moment_kinetics/src/runge_kutta.jl#L537C1-L537C55

Quadmath is imported in a module other than runge_kutta.jl

[johnomotani]

Is it really true that we need Float128 types to have adaptive timestepping? Does having Float128 as a type slow down the speed of the calculation significantly?

I am referencing this comment in master, I assume it came from this PR? https://github.com/mabarnes/moment_kinetics/blob/1800f807308827a39f7e0bf0a3f5928a6990a0e1/moment_kinetics/src/runge_kutta.jl#L537C1-L537C55

Quadmath is imported in a module other than runge_kutta.jl

Using Float128 is not required, and not used by default. It can be helpful for debugging, to try and make simulation results exactly identical when running on different numbers of processes - see the description for high_precision_error_sum here: https://mabarnes.github.io/moment_kinetics/dev/timestepping/#timestepping-input-parameters.

It does look like using Quadmath should be in time_advance.jl rather than moment_kinetics_input.jl. I must've forgotten to move it when I changed where some of the input was set up - I'm actually not sure how the tests are passing at the moment, possibly only because they don't use the high_precision_error_sum option... Is this causing you a problem?

[mrhardman]

I haven't noticed a problem, except that I had to manually install Quadmath in the package manager after a git pull. When I used higher precision floats in the collision operator I noticed a big speed hit, hence my concern.

[johnomotani]

had to manually install Quadmath in the package manager after a git pull.

Ah, that makes sense. It should be enough to do pkg> resolve to pull in the new dependency.