Temporal integration in stdlib (_Integrate)

Projects
1

Motivation: There has been some rumours about a big rework of a widely used ocean model here where I work, and knowing how painful and intricate is to rework the time scheme of an ocean model, I wanted to see if I can come up with a design to ease the process.

Disclaimer: I am not the most acknowledged person when it comes to time integration, but I am studying them now. I have designed this addition with Claude (especially the examples), so it might be sub-optimal with respect to something coded by a real expert of the field, but still better than what I could have achieved alone.

Prior art that I plagiarised: I have taken inspiration from @jacobwilliams’s rklib (with permission), especially for the generation of the schemes. rklib is kept here as a golden standard, as I compare all the methods and results to what rklib does. The ODE model design I have chosen to refer to is that of @dubos’s GFmodels. This model has allowed me to feed more complicated models (e.g. a Shallow water model and a double gyre) to the same time scheme I was testing. No wonders I’ve read this discussion on an almost daily basis.

This being said, I would like to introduce the first iteration for the Integrate library

image
image1790×285 21 KB

stdlib_Integrate

Time integration of initial-value problems, dy/dt = f(t, y), y(t_0) = y_0, for Fortran, with explicit Runge-Kutta methods. Written as a candidate for the Fortran standard library: proposal.md is the full proposal, and the place to start.

The program describes what is integrated; the library provides how. There are two ways to describe the problem, and they run the same methods with the same steps:

  • plain arrays: the state is one rank-1 array, and the program writes only a right-hand side f(t, y, dydt, stage);
  • model classes: the program extends four abstract types (a model, its state, its tendencies, a scratch space), so that the state can be any data structure (several fields, multi-dimensional grids, halos, tiles), which the library never copies or reorganises. This is what models built on the method of lines need, and what existing Fortran ODE libraries do not offer.
use stdlib_Integrate, only: rk_scheme, rk_array_scheme_t, stage_t

type(rk_array_scheme_t) :: scheme

call rk_scheme(scheme, "rk4", size(y), dt)     ! the method, by name
do it = 1, n_steps
    call scheme%advance(f, time, y)             ! y and time, one step forward
end do

What is there

  • three methods at fixed steps: Wicker-Skamarock RK3, classic RK4, Dormand-Prince 5(4);
  • classic RK4 as a hardcoded step, and a generic engine that runs any explicit Butcher tableau, including one built by hand: every other method runs on it, and adding a method is adding a tableau;
  • a hook to fill halos and apply boundary conditions before every evaluation.

Not yet: adaptive steps, low-storage forms, and other families of methods (multistep, implicit); the proposal says how they would fit.

Double precision only for now. Working and tested with fpm and gfortran.

Looking forward to hear what you think about this API, and what you would need to fit it to your needs. The choice of the methods, the extension to different ranks and precisions can come later.

Cheers,
Francesco

4 Likes