RKDP54

README.md

@@ -93,8 +93,11 @@ should be a `;` sperated list of `key=value` pairs. The possible keys are
 
 * `random_seed` (int, default ): seed for random initialization.
 
-* `algorithm`: either `RK4` for Runge-Kutta 4, `RK2` for Runge-Kutta 2, or
-  `RKF45` for the Runge-Kutta-Fehlberg adaptive step method.
+* `algorithm`: fixed step methods: `RK4` for Runge-Kutta 4, `RK2` for
+  Runge-Kutta 2.
+  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,4 +34,6 @@ limitations under the License.
 #define ALGORITHM_ADAPTIVE_THRESHOLD 1000
 // adaptive algorithms: index > ALGORITHM_ADAPTIVE_THRESHOLD
 #define ALGORITHM_RKF45 1001
+#define ALGORITHM_RKDP54 1002
+#define ALGORITHM_RKBS32 1003
 

src/main.c

@@ -265,6 +265,12 @@ 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;
   default:
     fprintf(file,", algorithm=RK4");
     break;
@@ -671,6 +677,12 @@ 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;
+    }
     else{
       fprintf(stderr, "error: unrecognized algorithm '%s'\n",rhs);
       return(-1);

src/navier-stokes.c

@@ -89,7 +89,13 @@ int uk(
     } else if (algorithm==ALGORITHM_RK4) {
       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);
+      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);
     } else {
       fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
     }
@@ -192,7 +198,13 @@ int enstrophy(
     } else if (algorithm==ALGORITHM_RK4) {
       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);
+      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);
     } else {
       fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
     }
@@ -216,7 +228,7 @@ int enstrophy(
       avg_en_x_a*=print_freq/(time-prevtime);
 
       // print to stderr so user can follow along
-      if(algorithm==ALGORITHM_RKF45){
+      if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
 	fprintf(stderr,"% .8e % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",time, avg_a, avg_en, avg_en_x_a, alpha, enstrophy, alpha*enstrophy, step);
 	printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy, step);
       } else {
@@ -344,7 +356,13 @@ int quiet(
     } else if (algorithm==ALGORITHM_RK4) {
       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);
+      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);
     } else {
       fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
     }
@@ -392,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);
@@ -400,6 +418,12 @@ int ns_init_tmps(
     *tmp5=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
     *tmp6=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
     *tmp7=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
+  } else if (algorithm==ALGORITHM_RKBS32){
+    *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);
+    *tmp4=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
+    *tmp5=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
   } else {
     fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
   };
@@ -451,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);
@@ -459,6 +483,12 @@ int ns_free_tmps(
     free(tmp5);
     free(tmp6);
     free(tmp7);
+  } else if (algorithm==ALGORITHM_RKBS32){
+    free(tmp1);
+    free(tmp2);
+    free(tmp3);
+    free(tmp4);
+    free(tmp5);
   } else {
     fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
   };
@@ -631,13 +661,17 @@ int ns_step_rkf45(
   _Complex double* k5,
   _Complex double* k6,
   _Complex double* tmp,
-  bool irreversible
+  bool irreversible,
+  // whether to compute k1 or leave it as is
+  bool compute_k1
 ){
   int kx,ky;
   double err,relative;
 
   // k1: u(t)
+  if(compute_k1){
     ns_rhs(k1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
+  }
 
   // k2 : u(t+1/4*delta)
   for(kx=0;kx<=K1;kx++){
@@ -687,7 +721,7 @@ int ns_step_rkf45(
       // difference between 5th order and 4th order
       err+=cabs((*delta)*(1./360*k1[klookup_sym(kx,ky,K2)]-128./4275*k3[klookup_sym(kx,ky,K2)]-2197./75240*k4[klookup_sym(kx,ky,K2)]+1./50*k5[klookup_sym(kx,ky,K2)]+2./55*k6[klookup_sym(kx,ky,K2)]));
       // next step
-      tmp[klookup_sym(kx,ky,K2)]=(*delta)*(16./135*k1[klookup_sym(kx,ky,K2)]+6656./12825*k3[klookup_sym(kx,ky,K2)]+28561./56430*k4[klookup_sym(kx,ky,K2)]-9./50*k5[klookup_sym(kx,ky,K2)]+2./55*k6[klookup_sym(kx,ky,K2)]);
+      tmp[klookup_sym(kx,ky,K2)]=(*delta)*(25./216*k1[klookup_sym(kx,ky,K2)]+1408./2565*k3[klookup_sym(kx,ky,K2)]+2197./4104*k4[klookup_sym(kx,ky,K2)]-1./5*k5[klookup_sym(kx,ky,K2)]);
       relative+=cabs(tmp[klookup_sym(kx,ky,K2)]);
     }
   }
@@ -701,12 +735,239 @@ int ns_step_rkf45(
       }
     }
     // next delta to use in future steps
-    *next_delta=factor*(*delta)*pow(relative*tolerance/err,0.2);
+    *next_delta=(*delta)*pow(relative*tolerance/err,0.2);
   }
   // error too big: repeat with smaller step
   else{
     *delta=factor*(*delta)*pow(relative*tolerance/err,0.2);
-    ns_step_rkf45(u,tolerance,factor,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible);
+    // do not recompute k1
+    ns_step_rkf45(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;
+}
+
+
+// next time step
+// adaptive RK algorithm (Runge-Kutta-Bogacki-Shampine method)
+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,
+  // the pointers k1 and k4 will be exchanged at the end of the routine
+  _Complex double** k1,
+  _Complex double* k2,
+  _Complex double* k3,
+  _Complex double** k4,
+  _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/4*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)/2*(*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/4*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./4*k2[klookup_sym(kx,ky,K2)]);
+    }
+  }
+  ns_rhs(k3, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
+
+  // k4 : u(t+delta)
+  // tmp cpmputed here is the next step
+  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)*(2./9*(*k1)[klookup_sym(kx,ky,K2)]+1./3*k2[klookup_sym(kx,ky,K2)]+4./9*k3[klookup_sym(kx,ky,K2)]);
+    }
+  }
+  ns_rhs(*k4, 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)*(5./72*(*k1)[klookup_sym(kx,ky,K2)]-1./12*k2[klookup_sym(kx,ky,K2)]-1./9*k3[klookup_sym(kx,ky,K2)]+1./8*(*k4)[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./3);
+
+    // k1 in the next step will be this k4 (first same as last)
+    tmp=*k1;
+    *k1=*k4;
+    *k4=tmp;
+  }
+  // error too big: repeat with smaller step
+  else{
+    *delta=factor*(*delta)*pow(relative*tolerance/err,1./3);
+    // this will reuse the same k1 without re-computing it
+    ns_step_rkbs32(u,tolerance,factor,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,tmp,irreversible,false);
+  }
+
+  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;

src/navier-stokes.h

@@ -55,8 +55,12 @@ 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)
-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);
+// 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);
 
 // right side of Irreversible/reversible Navier-Stokes equation
 int ns_rhs( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, bool irreversible);