sandbox/Antoonvh/lsrk.h
Low-Storage RK-3 and RK-4 schemes
The default is 3-stage 3rd order method, if defined otherwise, a 5-stage 4th order method is used.
As suggested in,
van Heerwaarden, et al. MicroHH 1.0: … Geoscientific Model Development (2017).
#include "run.h"
#ifndef RKORDER
#define RKORDER (3)
#endif
#if (RKORDER == 3)
// Williamson, J. H.: Low-Storage Runge-Kutta schemes, J.
// Comput.Phys., 35, 48–56, 1980.
#define STAGES (3)
double An[STAGES] = {0., -5./9., -153./128.};
double Bn[STAGES] = {1./3., 15./16., 8./15.};
#else
// Carpenter, M. H. and Kennedy, C. A.: Fourth-order
// 2N-storageRunge-Kutta schemes, Tech. Rep. TM-109112, NASA
// LangleyResearch Center, 1994
#define STAGES (5)
double An[STAGES] = {0.,
-567301805773. /1357537059087.,
-2404267990393./2016746695238.,
-3550918686646./2091501179385.,
-1275806237668./842570457699.};
double Bn[STAGES] = {1432997174477./9575080441755. ,
5161836677717./13612068292357.,
1720146321549./2090206949498. ,
3134564353537./4481467310338. ,
2277821191437./14882151754819.};
#endif
scalar * dfl = NULL;
void A_Time_Step (scalar * fl, double dt,
void (* Lu) (scalar * ul, scalar * dul)) {
if (dfl == NULL)
dfl = list_clone (fl);
scalar * dfltmp = list_clone (dfl);
for (int Stp = 0; Stp < STAGES; Stp++) {
Lu (fl, dfltmp);
scalar f, df, dftmp;
foreach() {
for (f, df, dftmp in fl, dfl, dfltmp) {
df[] = An[Stp]*df[] + dftmp[];
f[] += Bn[Stp]*df[]*dt;
}
}
}
delete (dfltmp); free (dfltmp); dfltmp = NULL;
}
event rm_dfl (t = end) {
delete (dfl); free (dfl); dfl = NULL;
}