max_delta and adaptive norm in rkf45 and rkbs32

README.md

@@ -105,6 +105,14 @@ should be a `;` sperated list of `key=value` pairs. The possible keys are
 * `adaptive_factor` (double, default 0.9): when using the RKF45 method, the
   step gets adjusted by `factor*delta*(tolerance/error)^(1/5)`.
 
+* `max_delta` (double, default 1e-2): when using the adaptive step, do not
+  exceet `delta_max`.
+
+* `adaptive_norm`: norm to use to estimate the error of the adaptive method:
+  `L1` (default) for the normalized L1 norm, `k3` for the normalized L1 norm of
+  f_k/|k|^3, `k32` for the normalized L1 norm, `enstrophy` for the enstrophy
+  norm.
+
 
 # Interrupting/resuming the computation
 

docs/nstrophy_doc.tex

@@ -103,6 +103,103 @@ We truncate the Fourier modes and assume that $\hat u_k=0$ if $|k_1|>K_1$ or $|k
 \end{equation}
 \bigskip
 
+\point{\bf Runge-Kutta methods}.
+To solve the equation numerically, we will use Runge-Kutta methods, which compute an approximate value $\hat u_k^{(n)}$ for $\hat u_k(t_n)$.
+{\tt nstrophy} supports the 4th order Runge-Kutta ({\tt RK4}) and 2nd order Runge-Kutta ({\tt RK2}) algorithms.
+In addition, several variable step methods are implemented:
+\begin{itemize}
+  \item the Runge-Kutta-Dormand-Prince method ({\tt RKDP54}), which is of 5th order, and adjusts the step by comparing to a 4th order method;
+  \item the Runge-Kutta-Fehlberg method ({\tt RKF45}), which is of 4th order, and adjusts the step by comparing to a 5th order method;
+  \item the Runge-Kutta-Bogacki-Shampine method ({\tt RKBS32}), which is of 3d order, and adjusts the step by comparing to a 2nd order method.
+\end{itemize}
+In these adaptive step methods, two steps are computed at different orders: $\hat u_k^{(n)}$ and $\hat U_k^{(n)}$, the step size is adjusted at every step in such a way that the error is small enough:
+\begin{equation}
+  \|\hat u^{(n)}-\hat U^{(n)}\|
+  <\epsilon_{\mathrm{target}}
+\end{equation}
+for some given $\epsilon_{\mathrm{target}}$, set using the {\tt adaptive\_tolerance} parameter.
+The choice of the norm matters, and will be discussed below.
+If the error is larger than the target, then the step size is decreased.
+How this is done depends on the order of algorithm.
+If the order is $q$ (here we mean the smaller of the two orders, so 4 for {\tt RKDP54} and {\tt RKF45} and 2 for {\tt RKBS32}), then we expect
+\begin{equation}
+  \|\hat u^{(n)}-\hat U^{(n)}\|=\delta_n^qC_n
+  .
+\end{equation}
+We wish to set $\delta_{n+1}$ so that
+\begin{equation}
+  \delta_{n+1}^qC_n=\epsilon_{\mathrm{target}}
+\end{equation}
+so
+\begin{equation}
+  \delta_{n+1}
+  =\left(\frac{\epsilon_{\mathrm{target}}}{C_n}\right)^{\frac1q}
+  =\delta_n\left(\frac{\epsilon_{\mathrm{target}}}{\|\hat u^{(n)}-\hat U^{(n)}\|}\right)^{\frac1q}
+  .
+  \label{adaptive_delta}
+\end{equation}
+(Actually, to be safe and ensure that $\delta$ decreases sufficiently, we multiply this by a safety factor that can be set using the {\tt adaptive\_factor} parameter.)
+If the error is smaller than the target, we increase $\delta$ using\-~(\ref{adaptive_delta}) (without the safety factor).
+To be safe, we also set a maximal value for $\delta$ via the {\tt max\_delta} parameter.
+\bigskip
+
+\indent
+The choice of the norm $\|\cdot\|$ matters.
+It can be made by specifying the parameter {\tt adaptive\_norm}.
+\begin{itemize}
+  \item A naive choice is to take $\|\cdot\|$ to be the normalized $L_1$ norm:
+  \begin{equation}
+    \|f\|:=
+    \frac1{\mathcal N}\sum_k|f_k|
+    ,\quad
+    \mathcal N:=\sum_k|\hat u_k^{(n)}-\hat u_k^{(n-1)}|
+    .
+  \end{equation}
+  This norm is selected by choosing {\tt adaptive\_norm=L1}.
+
+  \item Empirically, we have found that $|\hat u-\hat U|$ behaves like $k^{-3}$ for {\tt RKDP54} and {\tt RKF45}, and like $k^{-\frac32}$ for {\tt RKBS32}, so a norm of the form
+  \begin{equation}
+    \|f\|:=\frac1{\mathcal N}\sum_k|f_k|k^{-3}
+    ,\quad
+    \mathcal N:=\sum_k|\hat u_k^{(n)}-\hat u_k^{(n-1)}|k^{-3}
+  \end{equation}
+  or
+  \begin{equation}
+    \|f\|:=\frac1{\mathcal N}\sum_k|f_k|k^{-\frac32}
+    ,\quad
+    \mathcal N:=\sum_k|\hat u_k^{(n)}-\hat u_k^{(n-1)}|k^{-\frac32}
+  \end{equation}
+  are sensible choices.
+  These norms are selected by choosing {\tt adaptive\_norm=k3} and {\tt adaptive\_norm=k32} respectively.
+
+  \item
+  Another option is to define a norm based on the expression of the enstrophy\-~(\ref{enstrophy}):
+  \begin{equation}
+    \|f\|:=\frac1{\mathcal N}\sqrt{\sum_k k^2|f_k|^2}
+    \quad
+    \mathcal N:=\left(\frac{\sqrt{\sum_k k^2|\hat u_k^{(n)}|^2}+\sqrt{\sum_k k^2|\hat U_k^{(n)}|^2}}{\sum_k k^2|\hat u_k^{(n)}|^2}\right)^{\frac13}
+    .
+  \end{equation}
+  Doing so controls the error of the enstrophy through
+  \begin{equation}
+    \mathcal N^2|\mathcal En(\hat u)-\mathcal En(\hat U)|\equiv|\|\hat u\|^2-\|\hat U\|^2|\leqslant \|\hat u-\hat U\|(\|\hat u\|+\|\hat U\|)
+  \end{equation}
+  so
+  \begin{equation}
+    \mathcal N^2
+    |\mathcal En(\hat u)-\mathcal En(\hat U)|\leqslant
+    \|\hat u-\hat U\|\frac1{\mathcal N}\left(\sqrt{\sum_k k^2|\hat u_k|^2}+\sqrt{\sum_k k^2|\hat U_k|^2}\right)
+  \end{equation}
+  and thus
+  \begin{equation}
+    \frac{|\mathcal En(\hat u)-\mathcal En(\hat U)|}{\mathcal En(\hat u)}\leqslant
+    \|\hat u-\hat U\|
+    .
+  \end{equation}
+  This norm is selected by choosing {\tt adaptive\_norm=enstrophy}.
+\end{itemize}
+\bigskip
+
 \point{\bf Reality}.
 Since $U$ is real, $\hat U_{-k}=\hat U_k^*$, and so
 \begin{equation}
@@ -271,6 +368,7 @@ and
     +e^{-\frac{4\pi^2}{L^2}\nu t}\sqrt{E(0)}
   \right)^2
   .
+  \label{enstrophy}
 \end{equation}
 \bigskip
 

src/complex_tools.c

@@ -0,0 +1,22 @@
+/*
+Copyright 2017-2023 Ian Jauslin
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include "complex_tools.h"
+
+// magnitude squared
+double cabs2(_Complex double x){
+  return (__real__ x)*(__real__ x)+(__imag__ x)*(__imag__ x);
+}

src/complex_tools.h

@@ -0,0 +1,23 @@
+/*
+Copyright 2017-2023 Ian Jauslin
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#ifndef COMPLEXTOOLS_H
+#define COMPLEXTOOLS_H
+
+// magnitude squared
+double cabs2(_Complex double x);
+
+#endif

src/constants.cpp

@@ -37,3 +37,7 @@ limitations under the License.
 #define ALGORITHM_RKDP54 1002
 #define ALGORITHM_RKBS32 1003
 
+#define NORM_L1 1
+#define NORM_k3 2
+#define NORM_k32 3
+#define NORM_ENSTROPHY 4

src/main.c

@@ -47,6 +47,8 @@ typedef struct nstrophy_parameters {
   double L;
   double adaptive_tolerance;
   double adaptive_factor;
+  double max_delta;
+  unsigned int adaptive_norm;
   double print_freq;
   int seed;
   double starting_time;
@@ -170,16 +172,16 @@ int main (
 
   // run command
   if (command==COMMAND_UK){
-    uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, u0, g, parameters.irreversible, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile);
+    uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_norm, u0, g, parameters.irreversible, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile);
   }
   else if(command==COMMAND_ENSTROPHY){
     // register signal handler to handle aborts
     signal(SIGINT, sig_handler);
     signal(SIGTERM, sig_handler);
-    enstrophy(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, u0, g, parameters.irreversible, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile, (char*)argv[0], param_str, savefile_str);
+    enstrophy(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_norm, u0, g, parameters.irreversible, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile, (char*)argv[0], param_str, savefile_str);
   }
   else if(command==COMMAND_QUIET){
-    quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.starting_time, u0, g, parameters.irreversible, parameters.algorithm, nthreads, savefile);
+    quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_norm, parameters.starting_time, u0, g, parameters.irreversible, parameters.algorithm, nthreads, savefile);
   }
   else if(command==0){
     fprintf(stderr, "error: no command specified\n");
@@ -264,19 +266,36 @@ int print_params(
     fprintf(file,", algorithm=RK2");
     break;
   case ALGORITHM_RKF45:
-    fprintf(file,", algorithm=RKF45, tolerance=%.15e, factor=%.15e",parameters.adaptive_tolerance, parameters.adaptive_factor);
+    fprintf(file,", algorithm=RKF45, tolerance=%.15e",parameters.adaptive_tolerance);
     break;
   case ALGORITHM_RKDP54:
-    fprintf(file,", algorithm=RKDP54, tolerance=%.15e, factor=%.15e",parameters.adaptive_tolerance, parameters.adaptive_factor);
+    fprintf(file,", algorithm=RKDP54, tolerance=%.15e",parameters.adaptive_tolerance);
     break;
   case ALGORITHM_RKBS32:
-    fprintf(file,", algorithm=RKBS32, tolerance=%.15e, factor=%.15e",parameters.adaptive_tolerance, parameters.adaptive_factor);
+    fprintf(file,", algorithm=RKBS32, tolerance=%.15e",parameters.adaptive_tolerance);
     break;
   default:
     fprintf(file,", algorithm=RK4");
     break;
   }
 
+  if(parameters.algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
+    switch(parameters.adaptive_norm){
+    case NORM_L1:
+      fprintf(file,", norm=L1");
+      break;
+    case NORM_k3:
+      fprintf(file,", norm=k3");
+      break;
+    case NORM_k32:
+      fprintf(file,", norm=k32");
+      break;
+    case NORM_ENSTROPHY:
+      fprintf(file,", norm=enstrophy");
+      break;
+    }
+  }
+
   fprintf(file,"\n");
 
   return 0;
@@ -397,6 +416,8 @@ int read_params(
   parameters->L=2*M_PI;
   parameters->adaptive_tolerance=1e-11;
   parameters->adaptive_factor=0.9;
+  parameters->max_delta=1e-2;
+  parameters->adaptive_norm=NORM_L1;
   parameters->final_time=100000;
   parameters->print_freq=1;
   parameters->starting_time=0;
@@ -597,6 +618,31 @@ int set_parameter(
       return(-1);
     }
   }
+  else if (strcmp(lhs,"max_delta")==0){
+    ret=sscanf(rhs,"%lf",&(parameters->max_delta));
+    if(ret!=1){
+      fprintf(stderr, "error: parameter 'max_delta' should be a double\n       got '%s'\n",rhs);
+      return(-1);
+    }
+  }
+  else if (strcmp(lhs,"adaptive_norm")==0){
+    if (strcmp(rhs,"L1")==0){
+      parameters->adaptive_norm=NORM_L1;
+    }
+    else if (strcmp(rhs,"k3")==0){
+      parameters->adaptive_norm=NORM_k3;
+    }
+    else if (strcmp(rhs,"k32")==0){
+      parameters->adaptive_norm=NORM_k32;
+    }
+    else if (strcmp(rhs,"enstrophy")==0){
+      parameters->adaptive_norm=NORM_ENSTROPHY;
+    }
+    else{
+      fprintf(stderr, "error: unrecognized adaptive_norm '%s'\n",rhs);
+      return(-1);
+    }
+  }
   else if (strcmp(lhs,"print_freq")==0){
     ret=sscanf(rhs,"%lf",&(parameters->print_freq));
     if(ret!=1){

src/navier-stokes.c

@@ -17,6 +17,7 @@ limitations under the License.
 #include "constants.cpp"
 #include "io.h"
 #include "navier-stokes.h"
+#include "complex_tools.h"
 #include <math.h>
 #include <stdint.h>
 #include <stdlib.h>
@@ -34,6 +35,8 @@ int uk(
   double L,
   double adaptive_tolerance,
   double adaptive_factor,
+  double max_delta,
+  unsigned int adaptive_norm,
   _Complex double* u0,
   _Complex double* g,
   bool irreversible,
@@ -89,13 +92,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, true);
+      ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, 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);
+      ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, 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);
+      ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, 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);
     }
@@ -142,6 +145,8 @@ int enstrophy(
   double L,
   double adaptive_tolerance,
   double adaptive_factor,
+  double max_delta,
+  unsigned int adaptive_norm,
   _Complex double* u0,
   _Complex double* g,
   bool irreversible,
@@ -198,13 +203,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, true);
+      ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, 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);
+      ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, 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);
+      ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, 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);
     }
@@ -318,6 +323,8 @@ int quiet(
   double L,
   double adaptive_tolerance,
   double adaptive_factor,
+  double max_delta,
+  unsigned int adaptive_norm,
   double starting_time,
   _Complex double* u0,
   _Complex double* g,
@@ -356,13 +363,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, true);
+      ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, 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);
+      ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, 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);
+      ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, 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);
     }
@@ -642,6 +649,8 @@ int ns_step_rkf45(
   _Complex double* u,
   double tolerance,
   double factor,
+  double max_delta,
+  unsigned int adaptive_norm,
   int K1,
   int K2,
   int N1,
@@ -715,17 +724,59 @@ int ns_step_rkf45(
   
   // 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
-      // use the norm |u_i|/k^3
-      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)]))/pow(kx*kx+ky*ky,1.5);
-      // next step
-      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)])/pow(kx*kx+ky*ky,1.5);
+  if(adaptive_norm==NORM_L1){
+    relative=0;
+    for(kx=0;kx<=K1;kx++){
+      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+	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)*(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)]);
+      }
     }
   }
+  else if(adaptive_norm==NORM_k3){
+    relative=0;
+    for(kx=0;kx<=K1;kx++){
+      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+	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)]))/pow(kx*kx+ky*ky,1.5);
+	// next step
+	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)])/pow(kx*kx+ky*ky,1.5);
+      }
+    }
+  }
+  else if(adaptive_norm==NORM_k32){
+    relative=0;
+    for(kx=0;kx<=K1;kx++){
+      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+	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)]))/pow(kx*kx+ky*ky,0.75);
+	// next step
+	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)])/pow(kx*kx+ky*ky,0.75);
+      }
+    }
+  }
+  else if(adaptive_norm==NORM_ENSTROPHY){
+    double sumu, sumU;
+    sumu=0;
+    sumU=0;
+    for(kx=0;kx<=K1;kx++){
+      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+	err+=(kx*kx+ky*ky)*cabs2((*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)*(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)]);
+	sumU+=(kx*kx+ky*ky)*cabs2((*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)]));
+	sumu+=(kx*kx+ky*ky)*cabs2(u[klookup_sym(kx,ky,K2)]+tmp[klookup_sym(kx,ky,K2)]);
+      }
+    }
+    err=sqrt(err);
+    relative=pow((sqrt(sumu)+sqrt(sumU))/sumu, 1./3);
+  }
+  else{
+    fprintf(stderr,"bug: unknown norm: %u, contact ian.jauslin@rutgers.edu\n", adaptive_norm);
+    exit(-1);
+  }
 
   // compare relative error with tolerance
   if(err<relative*tolerance){
@@ -735,14 +786,14 @@ int ns_step_rkf45(
 	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,0.2);
+    // next delta to use in future steps (do not exceed max_delta)
+    *next_delta=fmin(max_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);
     // 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);
+    ns_step_rkf45(u,tolerance,factor,max_delta,adaptive_norm,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible,false);
   }
 
   return 0;
@@ -755,6 +806,8 @@ int ns_step_rkbs32(
   _Complex double* u,
   double tolerance,
   double factor,
+  double max_delta,
+  unsigned int adaptive_norm,
   int K1,
   int K2,
   int N1,
@@ -813,14 +866,51 @@ int ns_step_rkbs32(
   
   // 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)]))/pow(kx*kx+ky*ky,0.75);
-      relative+=cabs(tmp[klookup_sym(kx,ky,K2)]-u[klookup_sym(kx,ky,K2)])/pow(kx*kx+ky*ky,0.75);
+  if(adaptive_norm==NORM_L1){
+    relative=0;
+    for(kx=0;kx<=K1;kx++){
+      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+	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)]);
+      }
     }
   }
+  else if(adaptive_norm==NORM_k3){
+    relative=0;
+    for(kx=0;kx<=K1;kx++){
+      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+	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)]))/pow(kx*kx+ky*ky,1.5);
+	relative+=cabs(tmp[klookup_sym(kx,ky,K2)]-u[klookup_sym(kx,ky,K2)])/pow(kx*kx+ky*ky,1.5);
+      }
+    }
+  }
+  else if(adaptive_norm==NORM_k32){
+    relative=0;
+    for(kx=0;kx<=K1;kx++){
+      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+	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)]))/pow(kx*kx+ky*ky,0.75);
+	relative+=cabs(tmp[klookup_sym(kx,ky,K2)]-u[klookup_sym(kx,ky,K2)])/pow(kx*kx+ky*ky,0.75);
+      }
+    }
+  }
+  else if(adaptive_norm==NORM_ENSTROPHY){
+    double sumu, sumU;
+    sumu=0;
+    sumU=0;
+    for(kx=0;kx<=K1;kx++){
+      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+	err+=(kx*kx+ky*ky)*cabs2((*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)]));
+	sumU+=(kx*kx+ky*ky)*cabs2((*delta)*(7./24*(*k1)[klookup_sym(kx,ky,K2)]+1./4*k2[klookup_sym(kx,ky,K2)]+1./3*k3[klookup_sym(kx,ky,K2)]+1./8*(*k4)[klookup_sym(kx,ky,K2)]));
+	sumu+=(kx*kx+ky*ky)*cabs2(tmp[klookup_sym(kx,ky,K2)]);
+      }
+    }
+    err=sqrt(err);
+    relative=pow((sqrt(sumu)+sqrt(sumU))/sumu, 1./3);
+  }
+  else{
+    fprintf(stderr,"bug: unknown norm: %u, contact ian.jauslin@rutgers,edu\n", adaptive_norm);
+    exit(-1);
+  }
 
   // compare relative error with tolerance
   if(err<relative*tolerance){
@@ -830,8 +920,8 @@ int ns_step_rkbs32(
 	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);
+    // next delta to use in future steps (do not exceed max_delta)
+    *next_delta=fmin(max_delta, (*delta)*pow(relative*tolerance/err,1./3));
 
     // k1 in the next step will be this k4 (first same as last)
     tmp=*k1;
@@ -842,7 +932,7 @@ int ns_step_rkbs32(
   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);
+    ns_step_rkbs32(u,tolerance,factor,max_delta,adaptive_norm,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,tmp,irreversible,false);
   }
 
   return 0;
@@ -854,6 +944,8 @@ int ns_step_rkdp54(
   _Complex double* u,
   double tolerance,
   double factor,
+  double max_delta,
+  unsigned int adaptive_norm,
   int K1,
   int K2,
   int N1,
@@ -939,15 +1031,51 @@ int ns_step_rkdp54(
   
   // 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
-      // use the norm |u_i|/k^3
-      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)]))/pow(kx*kx+ky*ky,1.5);
-      relative+=cabs(tmp[klookup_sym(kx,ky,K2)]-u[klookup_sym(kx,ky,K2)])/pow(kx*kx+ky*ky,1.5);
+  if(adaptive_norm==NORM_L1){
+    relative=0;
+    for(kx=0;kx<=K1;kx++){
+      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+	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)]);
+      }
     }
   }
+  else if(adaptive_norm==NORM_k3){
+    relative=0;
+    for(kx=0;kx<=K1;kx++){
+      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+	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)]))/pow(kx*kx+ky*ky,1.5);
+	relative+=cabs(tmp[klookup_sym(kx,ky,K2)]-u[klookup_sym(kx,ky,K2)])/pow(kx*kx+ky*ky,1.5);
+      }
+    }
+  }
+  else if(adaptive_norm==NORM_k32){
+    relative=0;
+    for(kx=0;kx<=K1;kx++){
+      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+	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)]))/pow(kx*kx+ky*ky,0.75);
+	relative+=cabs(tmp[klookup_sym(kx,ky,K2)]-u[klookup_sym(kx,ky,K2)])/pow(kx*kx+ky*ky,0.75);
+      }
+    }
+  }
+  else if(adaptive_norm==NORM_ENSTROPHY){
+    double sumu, sumU;
+    sumu=0;
+    sumU=0;
+    for(kx=0;kx<=K1;kx++){
+      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+	err+=(kx*kx+ky*ky)*cabs2((*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)]));
+	sumU+=(kx*kx+ky*ky)*cabs2(u[klookup_sym(kx,ky,K2)]+(*delta)*(5179./57600*(*k1)[klookup_sym(kx,ky,K2)]+7571./16695*k3[klookup_sym(kx,ky,K2)]+393./640*k4[klookup_sym(kx,ky,K2)]-92097./339200*k5[klookup_sym(kx,ky,K2)]+187./2100*k6[klookup_sym(kx,ky,K2)]+1./40*(*k2)[klookup_sym(kx,ky,K2)]));
+	sumu+=(kx*kx+ky*ky)*cabs2(tmp[klookup_sym(kx,ky,K2)]);
+      }
+    }
+    err=sqrt(err);
+    relative=pow((sqrt(sumu)+sqrt(sumU))/sumu, 1./3);
+  }
+  else{
+    fprintf(stderr,"bug: unknown norm: %u, contact ian.jauslin@rutgers,edu\n", adaptive_norm);
+    exit(-1);
+  }
 
   // compare relative error with tolerance
   if(err<relative*tolerance){
@@ -957,8 +1085,8 @@ int ns_step_rkdp54(
 	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);
+    // next delta to use in future steps (do not exceed max_delta)
+    *next_delta=fmin(max_delta, (*delta)*pow(relative*tolerance/err,1./5));
 
     // k1 in the next step will be this k4 (first same as last)
     tmp=*k1;
@@ -969,7 +1097,7 @@ int ns_step_rkdp54(
   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);
+    ns_step_rkdp54(u,tolerance,factor,max_delta,adaptive_norm,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

@@ -33,13 +33,13 @@ typedef struct fft_vects {
 } fft_vect;
 
 // compute u_k
-int uk( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, double print_freq, double starting_time, unsigned int nthreadsl, FILE* savefile);
+int uk( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, double print_freq, double starting_time, unsigned int nthreadsl, FILE* savefile);
 
 // compute enstrophy and alpha
-int enstrophy( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, double print_freq, double starting_time, unsigned int nthreads, FILE* savefile, char* cmd_string, char* params_string, char* savefile_string);
+int enstrophy( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, double print_freq, double starting_time, unsigned int nthreads, FILE* savefile, char* cmd_string, char* params_string, char* savefile_string);
 
 // compute solution as a function of time, but do not print anything (useful for debugging)
-int quiet( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double starting_time, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, unsigned int nthreads, FILE* savefile);
+int quiet( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, double starting_time, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, unsigned int nthreads, FILE* savefile);
 
 
 // initialize vectors for computation
@@ -56,11 +56,11 @@ int ns_step_rk4( _Complex double* u, int K1, int K2, int N1, int N2, double nu,
 // 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);
 // 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);
+int ns_step_rkf45( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_norm, 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);
+int ns_step_rkdp54( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_norm, 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);
+int ns_step_rkbs32( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_norm, 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);