RKDP54

README.md

@@ -93,9 +93,11 @@ should be a `;` sperated list of `key=value` pairs. The possible keys are
 
 * `random_seed` (int, default ): seed for random initialization.
 
-* `algorithm`: `RK4` for Runge-Kutta 4, `RK2` for Runge-Kutta 2, `RKF45` for
+* `algorithm`: fixed step methods: `RK4` for Runge-Kutta 4, `RK2` for
-  the Runge-Kutta-Fehlberg adaptive step method, `RKBS23` for the
+  Runge-Kutta 2.
-  Runge-Kutta-Bogacki-Shampine method.
+  adaptive step methods: `RKF45` for Runge-Kutta-Fehlberg (order
+  4), `RKDP54` for Runge-Kutta-Dormand-Prince (order 5), `RKBS32` for
+  Runge-Kutta-Bogacki-Shampine (order 3).
 
 * `adaptive_tolerance` (double, default 1e-11): when using an adaptive step
   method, this is the maximal allowed relative error.

src/constants.cpp

@@ -34,5 +34,6 @@ limitations under the License.
 #define ALGORITHM_ADAPTIVE_THRESHOLD 1000
 // adaptive algorithms: index > ALGORITHM_ADAPTIVE_THRESHOLD
 #define ALGORITHM_RKF45 1001
-#define ALGORITHM_RKBS32 1002
+#define ALGORITHM_RKDP54 1002
+#define ALGORITHM_RKBS32 1003
 

src/main.c

@@ -265,6 +265,9 @@ int print_params(
   case ALGORITHM_RKF45:
     fprintf(file,", algorithm=RKF45, tolerance=%.15e, factor=%.15e",parameters.adaptive_tolerance, parameters.adaptive_factor);
     break;
+  case ALGORITHM_RKDP54:
+    fprintf(file,", algorithm=RKDP54, tolerance=%.15e, factor=%.15e",parameters.adaptive_tolerance, parameters.adaptive_factor);
+    break;
   case ALGORITHM_RKBS32:
     fprintf(file,", algorithm=RKBS32, tolerance=%.15e, factor=%.15e",parameters.adaptive_tolerance, parameters.adaptive_factor);
     break;
@@ -674,6 +677,9 @@ int set_parameter(
     else if (strcmp(rhs,"RKF45")==0){
       parameters->algorithm=ALGORITHM_RKF45;
     }
+    else if (strcmp(rhs,"RKDP54")==0){
+      parameters->algorithm=ALGORITHM_RKDP54;
+    }
     else if (strcmp(rhs,"RKBS32")==0){
       parameters->algorithm=ALGORITHM_RKBS32;
     }

src/navier-stokes.c

@@ -90,6 +90,9 @@ int uk(
       ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
     } else if (algorithm==ALGORITHM_RKF45) {
       ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, true);
+    } else if (algorithm==ALGORITHM_RKDP54) {
+      // only compute k1 on the first step
+      ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, time==starting_time);
     } else if (algorithm==ALGORITHM_RKBS32) {
       // only compute k1 on the first step
       ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, time==starting_time);
@@ -196,6 +199,9 @@ int enstrophy(
       ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
     } else if (algorithm==ALGORITHM_RKF45) {
       ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, true);
+    } else if (algorithm==ALGORITHM_RKDP54) {
+      // only compute k1 on the first step
+      ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, time==starting_time);
     } else if (algorithm==ALGORITHM_RKBS32) {
       // only compute k1 on the first step
       ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, time==starting_time);
@@ -351,6 +357,9 @@ int quiet(
       ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
     } else if (algorithm==ALGORITHM_RKF45) {
       ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, true);
+    } else if (algorithm==ALGORITHM_RKDP54) {
+      // only compute k1 on the first step
+      ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, time==starting_time);
     } else if (algorithm==ALGORITHM_RKBS32) {
       // only compute k1 on the first step
       ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, time==starting_time);
@@ -401,7 +410,7 @@ int ns_init_tmps(
     *tmp1=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
     *tmp2=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
     *tmp3=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
-  } else if (algorithm==ALGORITHM_RKF45){
+  } else if (algorithm==ALGORITHM_RKF45 || algorithm==ALGORITHM_RKDP54){
     *tmp1=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
     *tmp2=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
     *tmp3=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
@@ -466,7 +475,7 @@ int ns_free_tmps(
     free(tmp1);
     free(tmp2);
     free(tmp3);
-  } else if (algorithm==ALGORITHM_RKF45){
+  } else if (algorithm==ALGORITHM_RKF45 || algorithm==ALGORITHM_RKDP54){
     free(tmp1);
     free(tmp2);
     free(tmp3);
@@ -838,6 +847,132 @@ int ns_step_rkbs32(
   return 0;
 }
 
+// next time step
+// adaptive RK algorithm (Runge-Kutta-Dormand-Prince method)
+int ns_step_rkdp54(
+  _Complex double* u,
+  double tolerance,
+  double factor,
+  int K1,
+  int K2,
+  int N1,
+  int N2,
+  double nu,
+  double* delta,
+  double* next_delta,
+  double L,
+  _Complex double* g,
+  fft_vect fft1,
+  fft_vect fft2,
+  fft_vect ifft,
+  // the pointers k1 and k2 will be exchanged at the end of the routine
+  _Complex double** k1,
+  _Complex double** k2,
+  _Complex double* k3,
+  _Complex double* k4,
+  _Complex double* k5,
+  _Complex double* k6,
+  _Complex double* tmp,
+  bool irreversible,
+  // whether to compute k1
+  bool compute_k1
+){
+  int kx,ky;
+  double err,relative;
+
+  // k1: u(t)
+  // only compute it if it is the first step (otherwise, it has already been computed due to the First Same As Last property)
+  if(compute_k1){
+    ns_rhs(*k1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
+  }
+
+  // k2 : u(t+1/5*delta)
+  for(kx=0;kx<=K1;kx++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)/5*(*k1)[klookup_sym(kx,ky,K2)];
+    }
+  }
+  ns_rhs(*k2, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
+
+  // k3 : u(t+3/10*delta)
+  for(kx=0;kx<=K1;kx++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(3./40*(*k1)[klookup_sym(kx,ky,K2)]+9./40*(*k2)[klookup_sym(kx,ky,K2)]);
+    }
+  }
+  ns_rhs(k3, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
+
+  // k4 : u(t+4/5*delta)
+  for(kx=0;kx<=K1;kx++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(44./45*(*k1)[klookup_sym(kx,ky,K2)]-56./15*(*k2)[klookup_sym(kx,ky,K2)]+32./9*k3[klookup_sym(kx,ky,K2)]);
+    }
+  }
+  ns_rhs(k4, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
+
+  // k5 : u(t+8/9*delta)
+  for(kx=0;kx<=K1;kx++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(19372./6561*(*k1)[klookup_sym(kx,ky,K2)]-25360./2187*(*k2)[klookup_sym(kx,ky,K2)]+64448./6561*k3[klookup_sym(kx,ky,K2)]-212./729*k4[klookup_sym(kx,ky,K2)]);
+    }
+  }
+  ns_rhs(k5, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
+
+  // k6 : u(t+delta)
+  for(kx=0;kx<=K1;kx++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(9017./3168*(*k1)[klookup_sym(kx,ky,K2)]-355./33*(*k2)[klookup_sym(kx,ky,K2)]+46732./5247*k3[klookup_sym(kx,ky,K2)]+49./176*k4[klookup_sym(kx,ky,K2)]-5103./18656*k5[klookup_sym(kx,ky,K2)]);
+    }
+  }
+  ns_rhs(k6, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
+
+  // k7 : u(t+delta)
+  for(kx=0;kx<=K1;kx++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+      // tmp computed here is the step
+      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(35./384*(*k1)[klookup_sym(kx,ky,K2)]+500./1113*k3[klookup_sym(kx,ky,K2)]+125./192*k4[klookup_sym(kx,ky,K2)]-2187./6784*k5[klookup_sym(kx,ky,K2)]+11./84*k6[klookup_sym(kx,ky,K2)]);
+    }
+  }
+  // store in k2, which is not needed anymore
+  ns_rhs(*k2, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
+  
+  // compute error
+  err=0;
+  relative=0;
+  for(kx=0;kx<=K1;kx++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+      // difference between 5th order and 4th order
+      err+=cabs((*delta)*(-71./57600*(*k1)[klookup_sym(kx,ky,K2)]+71./16695*k3[klookup_sym(kx,ky,K2)]-71./1920*k4[klookup_sym(kx,ky,K2)]+17253./339200*k5[klookup_sym(kx,ky,K2)]-22./525*k6[klookup_sym(kx,ky,K2)]+1./40*(*k2)[klookup_sym(kx,ky,K2)]));
+      relative+=cabs(tmp[klookup_sym(kx,ky,K2)]-u[klookup_sym(kx,ky,K2)]);
+    }
+  }
+
+  // compare relative error with tolerance
+  if(err<relative*tolerance){
+    // add to output
+    for(kx=0;kx<=K1;kx++){
+      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+	u[klookup_sym(kx,ky,K2)]=tmp[klookup_sym(kx,ky,K2)];
+      }
+    }
+    // next delta to use in future steps
+    *next_delta=(*delta)*pow(relative*tolerance/err,1./5);
+
+    // k1 in the next step will be this k4 (first same as last)
+    tmp=*k1;
+    *k1=*k2;
+    *k2=tmp;
+  }
+  // error too big: repeat with smaller step
+  else{
+    *delta=factor*(*delta)*pow(relative*tolerance/err,1./5);
+    // this will reuse the same k1 without re-computing it
+    ns_step_rkdp54(u,tolerance,factor,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible,false);
+  }
+
+  return 0;
+}
+
 // right side of Irreversible/Reversible Navier-Stokes equation
 int ns_rhs(
   _Complex double* out,

src/navier-stokes.h

@@ -55,8 +55,10 @@ int copy_u( _Complex double* u, _Complex double* u0, int K1, int K2);
 int ns_step_rk4( _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2,fft_vect ifft, _Complex double* tmp1, _Complex double *tmp2, _Complex double *tmp3, bool irreversible);
 // RK2
 int ns_step_rk2( _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2,fft_vect ifft, _Complex double* tmp1, _Complex double *tmp2, bool irreversible);
-// adaptive RK algorithm (Runge-Kutta-Fehlberg)
+// Runge-Kutta-Fehlberg
 int ns_step_rkf45( _Complex double* u, double tolerance, double factor, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double* k1, _Complex double* k2, _Complex double* k3, _Complex double* k4, _Complex double* k5, _Complex double* k6, _Complex double* tmp, bool irreversible, bool compute_k1);
+// Runge-Kutta-Dromand-Prince
+int ns_step_rkdp54( _Complex double* u, double tolerance, double factor, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double** k1, _Complex double** k2, _Complex double* k3, _Complex double* k4, _Complex double* k5, _Complex double* k6, _Complex double* tmp, bool irreversible, bool compute_k1);
 // Runge-Kutta-Bogacki-Shampine
 int ns_step_rkbs32( _Complex double* u, double tolerance, double factor, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double** k1, _Complex double* k2, _Complex double* k3, _Complex double** k4, _Complex double* tmp, bool irreversible, bool compute_k1);