Remove init_enstrophy and init_energy from savefile

In particular, removed the 'init_flow' parameter

README.md

@@ -40,12 +40,14 @@ The available commands are
 * `enstrophy`: to compute the enstrophy and various other observables. This
   command prints
   ```
-  step_index time average(alpha) average(energy) average(enstrophy) alpha energy enstrophy
+  step_index time average(alpha) average(energy) average(enstrophy) alpha energy enstrophy Re(u(1,1)) Re(u(1,2))
   ```
   where the averages are running averages over `print_freq`. In addition, if
   the algorithm has an adaptive step, an extra column is printed with `delta`.
   In addition, if alpha has a negative value (even in between `print_freq`
-  intervals), a line is printed with the information.
+  intervals), a line is printed with the information. The two components (1,1)
+  and (1,2) of u are included to more easily identify periodic trajectories, or
+  the presence of multiple attractors.
 
 * `lyapunov`: to compute the Lyapunov exponents. This command prints
   ```
@@ -56,7 +58,14 @@ The available commands are
 
 * `uk`: to compute the Fourier transform of the solution.
 
-* `quiet`: does not print anything, useful for debugging.
+* `quiet`: does not print anything (useful for debugging).
+
+* `enstrophy_print_init`: to compute the enstrophy and various other
+  observables for the initial condition (useful for debugging). The command
+  prints
+  ```
+  alpha energy enstrophy
+  ```
 
 
 # Parameters
@@ -148,12 +157,12 @@ should be a `;` sperated list of `key=value` pairs. The possible keys are
   is negative, its value is printed as a comment.
 
 * `lyapunov_reset` (double, default: `print_freq`): if this is set, then, when
-  computing the Lyapnuov exponents, the tangent flow will renormalize itself at
+  computing the Lyapunov exponents, the tangent flow will renormalize itself at
   times proportional to `lyapunov_reset`. This option is incompatible with
   `lyapunov_maxu`.
 
 * `lyapunov_maxu` (double, default: unset): if this is set, then, when
-  computing the Lyapnuov exponents, the tangent flow will renormalize itself
+  computing the Lyapunov exponents, the tangent flow will renormalize itself
   whenever the norm of the tangent flow exceeds  `lyapunov_maxu`.  This option
   is incompatible with `lyapunov_reset`.
 
@@ -164,14 +173,15 @@ should be a `;` sperated list of `key=value` pairs. The possible keys are
 
 # Interrupting and resuming the computation
 
-The `enstrophy` computation can be interrupted by sending Nstrophy the `SIGINT`
-signal (e.g. by pressing `Ctrl-C`.) When Nstrophy receives the `SIGINT` signal,
-it finishes its current step and writes the value of uk, either to `savefile`
-if such a file was specified on the command line (using the `-s` flag), or to
-`stderr`. In addition, when a `savefile` is specified it writes the command
-that needs to be used to resume the computation (which essentially just sets
-the appropriate `starting_time` and `init:file:<savefile>` parameters. The data
-written to the `savefile` is binary.
+The `enstrophy` and `lyapunov` computations can be interrupted by sending
+Nstrophy the `SIGINT` signal (e.g. by pressing `Ctrl-C`.) When Nstrophy
+receives the `SIGINT` signal, it finishes its current step and writes the value
+of uk, either to `savefile` if such a file was specified on the command line
+(using the `-s` flag), or to `stderr`. In addition, when a `savefile` is
+specified it writes the command that needs to be used to resume the computation
+(which essentially just sets the appropriate `starting_time` and
+`init:file:<savefile>` parameters). The data written to the `savefile` is
+binary.
 
 
 # License

src/constants.cpp

@@ -21,6 +21,7 @@ limitations under the License.
 #define COMMAND_QUIET 3
 #define COMMAND_RESUME 4
 #define COMMAND_LYAPUNOV 5
+#define COMMAND_ENSTROPHY_PRINT_INIT 6
 
 #define DRIVING_ZERO 1
 #define DRIVING_TEST 2

src/io.c

@@ -51,6 +51,30 @@ int write_vec_bin(_Complex double* vec, int K1, int K2, FILE* file){
   return 0;
 }
 
+// write complex vector (stored as 2 doubles) indexed by k1,k2 to file in binary format
+int write_vec2_bin(double* vec, int K1, int K2, FILE* file){
+  // do nothing if there is no file
+  if(file==NULL){
+    return 0;
+  }
+
+  fwrite(vec, sizeof(double), 2*(K1*(2*K2+1)+K2), file);
+
+  return 0;
+}
+
+// write complex matrix (stored as 2 doubles) indexed by k1,k2 to file in binary format
+int write_mat2_bin(double* mat, int K1, int K2, FILE* file){
+  // do nothing if there is no file
+  if(file==NULL){
+    return 0;
+  }
+
+  fwrite(mat, sizeof(double), 4*(K1*(2*K2+1)+K2)*(K1*(2*K2+1)+K2), file);
+
+  return 0;
+}
+
 // read complex vector indexed by k1,k2 from file
 int read_vec(_Complex double* out, int K1, int K2, FILE* file){
   int kx,ky;
@@ -149,15 +173,53 @@ int read_vec(_Complex double* out, int K1, int K2, FILE* file){
 
 // read complex vector indexed by k1,k2 from file in binary format
 int read_vec_bin(_Complex double* out, int K1, int K2, FILE* file){
-  char c;
-  int ret;
-
   // do nothing if there is no file
   if(file==NULL){
     return 0;
   }
 
   // seek past initial comments
+  seek_past_initial_comments(file);
+
+  fread(out, sizeof(_Complex double), K1*(2*K2+1)+K2, file);
+
+  return 0;
+}
+
+// read complex vector (represented as 2 doubles) indexed by k1,k2 from file in binary format
+int read_vec2_bin(double* out, int K1, int K2, FILE* file){
+  if(file==NULL){
+    return 0;
+  }
+
+  // seek past initial comments
+  seek_past_initial_comments(file);
+
+  fread(out, sizeof(double), 2*(K1*(2*K2+1)+K2), file);
+  return 0;
+}
+
+// read complex matrix (represented as 2 doubles) indexed by k1,k2 from file in binary format
+int read_mat2_bin(double* out, int K1, int K2, FILE* file){
+  if(file==NULL){
+    return 0;
+  }
+
+  // seek past initial comments
+  seek_past_initial_comments(file);
+
+  fread(out, sizeof(double), 4*(K1*(2*K2+1)+K2)*(K1*(2*K2+1)+K2), file);
+  return 0;
+}
+
+// ignore comments at beginning of file
+int seek_past_initial_comments(FILE* file){
+  char c;
+  int ret;
+
+  if(file==NULL){
+    return 0;
+  }
   while(true){
     ret=fscanf(file, "%c", &c);
     if (ret==1 && c=='#'){
@@ -184,8 +246,6 @@ int read_vec_bin(_Complex double* out, int K1, int K2, FILE* file){
     }
   }
 
-  fread(out, sizeof(_Complex double), K1*(2*K2+1)+K2, file);
-
   return 0;
 }
 
@@ -198,23 +258,38 @@ int remove_entry(
   char* rw_ptr;
   char* bfr;
   char* entry_ptr=entry;
+  // whether to write the entry
   int go=1;
+  // whether the pointer is at the beginning of the entry
+  int at_top=1;
 
   for(ptr=param_str, rw_ptr=ptr; *ptr!='\0'; ptr++){
+    // only match entries if one is at the beginning of an entry
+    if(at_top==1){
+      // check that the entry under ptr matches entry
       for(bfr=ptr,entry_ptr=entry; *bfr==*entry_ptr; bfr++, entry_ptr++){
 	// check if reached end of entry
 	if(*(bfr+1)=='=' && *(entry_ptr+1)=='\0'){
+	  // match: do not write entry
 	  go=0;
 	  break;
 	}
       }
+    }
+
+    // write entry
     if(go==1){
       *rw_ptr=*ptr;
       rw_ptr++;
     }
+
+    // next iterate will no longer be at the beginning of the entry
+    at_top=0;
+
     //reset
     if(*ptr==';'){
       go=1;
+      at_top=1;
     }
   }
   *rw_ptr='\0';
@@ -262,6 +337,8 @@ int save_state(
     strcpy(params, params_string);
     remove_entry(params, "starting_time");
     remove_entry(params, "init");
+    remove_entry(params, "init_enstrophy");
+    remove_entry(params, "init_energy");
     if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
       remove_entry(params, "delta");
     }

src/io.h

@@ -22,10 +22,21 @@ limitations under the License.
 // write complex vector indexed by k1,k2 to file
 int write_vec(_Complex double* u, int K1, int K2, FILE* file);
 int write_vec_bin(_Complex double* u, int K1, int K2, FILE* file);
+// write complex vector (stored as 2 doubles) indexed by k1,k2 to file in binary format
+int write_vec2_bin(double* vec, int K1, int K2, FILE* file);
+// write complex matrix (stored as 2 doubles) indexed by k1,k2 to file in binary format
+int write_mat2_bin(double* mat, int K1, int K2, FILE* file);
 
 // read complex vector indexed by k1,k2 from file
 int read_vec(_Complex double* u, int K1, int K2, FILE* file);
 int read_vec_bin(_Complex double* u, int K1, int K2, FILE* file);
+// read complex vector (represented as 2 doubles) indexed by k1,k2 from file in binary format
+int read_vec2_bin(double* out, int K1, int K2, FILE* file);
+// read complex matrix (represented as 2 doubles) indexed by k1,k2 from file in binary format
+int read_mat2_bin(double* out, int K1, int K2, FILE* file);
+
+// ignore comments at beginning of file
+int seek_past_initial_comments(FILE* file);
 
 // remove an entry from params string (inplace)
 int remove_entry(char* param_str, char* entry);

src/lyapunov.c

@@ -16,6 +16,7 @@ limitations under the License.
 
 #include "constants.cpp"
 #include "lyapunov.h"
+#include "io.h"
 #include <openblas64/cblas.h>
 #include <openblas64/lapacke.h>
 #include <math.h>
@@ -48,7 +49,18 @@ int lyapunov(
   unsigned int algorithm,
   unsigned int algorithm_lyapunov,
   double starting_time,
-  unsigned int nthreads
+  unsigned int nthreads,
+  double* flow0,
+  double* lyapunov_avg0,
+  double prevtime, // the previous time at which a QR decomposition was performed
+  double lyapunov_startingtime, // the time at which the lyapunov exponent computation was started (useful in interrupted computation)
+  FILE* savefile,
+  FILE* utfile,
+  // for interrupt recovery
+  const char* cmd_string,
+  const char* params_string,
+  const char* savefile_string,
+  const char* utfile_string
 ){
   double* lyapunov;
   double* lyapunov_avg;
@@ -82,24 +94,19 @@ int lyapunov(
 
   // copy initial condition
   copy_u(u, u0, K1, K2);
-  // initialize flow with identity matrix
+
+  // initialize flow and averages
   for (i=0;i<MATSIZE;i++){
     for (j=0;j<MATSIZE;j++){
-      if(i!=j){
-	flow[i*MATSIZE+j]=0.;
-      } else {
-	flow[i*MATSIZE+j]=1.;
-      }
+      flow[i*MATSIZE+j]=flow0[i*MATSIZE+j];
     }
+    lyapunov_avg[i]=lyapunov_avg0[i];
   }
 
   // store step (useful for adaptive step methods)
   double step=delta;
   double next_step=step;
 
-  // save times at which Lyapunov exponents are computed
-  double prevtime=starting_time;
-
   // iterate
   time=starting_time;
   while(final_time<0 || time<final_time){
@@ -118,19 +125,20 @@ int lyapunov(
       // size of flow (for reset)
       for(i=0;i<MATSIZE;i++){
 	for(j=0;j<MATSIZE;j++){
-	  norm+=fabs(flow[i*MATSIZE+j]);
-	  if(norm>lyapunov_reset){
+	  norm+=flow[i*MATSIZE+j]*flow[i*MATSIZE+j]/MATSIZE;
+	  if(sqrt(norm)>lyapunov_reset){
 	    break;
 	  }
 	}
-	if(norm>lyapunov_reset){
+	if(sqrt(norm)>lyapunov_reset){
 	  break;
 	}
       }
     }
 
     // QR decomposition
-    if((lyapunov_trigger==LYAPUNOV_TRIGGER_TIME && time>(n+1)*lyapunov_reset) || (lyapunov_trigger==LYAPUNOV_TRIGGER_SIZE && norm>lyapunov_reset)){
+    // Do it also if it is the last step
+    if((lyapunov_trigger==LYAPUNOV_TRIGGER_TIME && time>(n+1)*lyapunov_reset) || (lyapunov_trigger==LYAPUNOV_TRIGGER_SIZE && norm>lyapunov_reset) || time>=final_time){
       n++;
 
       // QR decomposition
@@ -141,14 +149,14 @@ int lyapunov(
 	lyapunov[i]=log(fabs(flow[i*MATSIZE+i]))/(time-prevtime);
       }
 
-      // sort lyapunov exponents
-      qsort(lyapunov, MATSIZE, sizeof(double), compare_double);
+      //// sort lyapunov exponents
+      //qsort(lyapunov, MATSIZE, sizeof(double), compare_double);
 
       // average lyapunov
       for(i=0; i<MATSIZE; i++){
 	// exclude inf
-	if((! isinf(lyapunov[i])) && (time>starting_time)){
-	  lyapunov_avg[i]=lyapunov_avg[i]*(prevtime-starting_time)/(time-starting_time)+lyapunov[i]*(time-prevtime)/(time-starting_time);
+	if((! isinf(lyapunov[i])) && (time>lyapunov_startingtime)){
+	  lyapunov_avg[i]=lyapunov_avg[i]*(prevtime-lyapunov_startingtime)/(time-lyapunov_startingtime)+lyapunov[i]*(time-prevtime)/(time-lyapunov_startingtime);
 	}
       }
 
@@ -157,11 +165,11 @@ int lyapunov(
 	printf("% .15e % .15e % .15e\n",time, lyapunov[i], lyapunov_avg[i]);
       }
       printf("\n");
-      fprintf(stderr,"% .15e",time);
-      // print largest and smallest lyapunov exponent to stderr
-      if(MATSIZE>0){
-	fprintf(stderr," % .15e % .15e\n", lyapunov[0], lyapunov[MATSIZE-1]);
-      }
+      fprintf(stderr,"% .15e\n",time);
+      //// print largest and smallest lyapunov exponent to stderr
+      //if(MATSIZE>0){
+      //  fprintf(stderr," % .15e % .15e\n", lyapunov[0], lyapunov[MATSIZE-1]);
+      //}
 
       // set flow to Q:
       LAPACKE_dorgqr(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, MATSIZE, flow, MATSIZE, tmp11);
@@ -169,6 +177,24 @@ int lyapunov(
       // reset prevtime
       prevtime=time;
     }
+
+    // catch abort signal
+    if (g_abort){
+      // print u to stderr if no savefile
+      if (savefile==NULL){
+	savefile=stderr;
+      }
+      break;
+    }
+  }
+
+  if(savefile!=NULL){
+    lyapunov_save_state(flow, u, lyapunov_avg, prevtime, lyapunov_startingtime, savefile, K1, K2, cmd_string, params_string, savefile_string, utfile_string, utfile, COMMAND_LYAPUNOV, algorithm, step, time, nthreads);
+  }
+
+  // save final u to utfile in txt format
+  if(utfile!=NULL){
+    write_vec(u, K1, K2, utfile);
   }
 
   lyapunov_free_tmps(lyapunov, lyapunov_avg, flow, u, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, tmp10, tmp11, fftu1, fftu2, fftu3, fftu4, fft1, ifft, algorithm, algorithm_lyapunov);
@@ -214,6 +240,7 @@ int lyapunov_D_step(
 ){
   int lx,ly;
   double alpha;
+  int n;
 
   // compute fft of u for future use
   lyapunov_fft_u_T(u,K1,K2,N1,N2,fftu1,fftu2,fftu3,fftu4);
@@ -237,17 +264,18 @@ int lyapunov_D_step(
   // loop over columns
   for(lx=0;lx<=K1;lx++){
     for(ly=(lx>0 ? -K2 : 1);ly<=K2;ly++){
-
+      for(n=0;n<=1;n++){
 	// do not use adaptive step algorithms for the tangent flow: it must be locked to the same times as u
 	if(algorithm==ALGORITHM_RK2){
-	lyapunov_D_step_rk2(flow+2*klookup_sym(lx,ly,K2)*MATSIZE, u, K1, K2, N1, N2, alpha, delta, L, g, tmp4, fftu1, fftu2, fftu3, fftu4, fft1, ifft, tmp1, tmp2, irreversible);
+	  lyapunov_D_step_rk2(flow+(2*klookup_sym(lx,ly,K2)+n)*MATSIZE, u, K1, K2, N1, N2, alpha, delta, L, g, tmp4, fftu1, fftu2, fftu3, fftu4, fft1, ifft, tmp1, tmp2, irreversible);
 	} else if (algorithm==ALGORITHM_RK4) {
-	lyapunov_D_step_rk4(flow+2*klookup_sym(lx,ly,K2)*MATSIZE, u, K1, K2, N1, N2, alpha, delta, L, g, tmp4, fftu1, fftu2, fftu3, fftu4, fft1, ifft, tmp1, tmp2, tmp3, irreversible);
+	  lyapunov_D_step_rk4(flow+(2*klookup_sym(lx,ly,K2)+n)*MATSIZE, u, K1, K2, N1, N2, alpha, delta, L, g, tmp4, fftu1, fftu2, fftu3, fftu4, fft1, ifft, tmp1, tmp2, tmp3, irreversible);
 	} else {
 	  fprintf(stderr,"bug: unknown algorithm for tangent flow: %u, contact ian.jauslin@rutgers.edu\n",algorithm);
 	}
       }
     }
+  }
 
   return(0);
 }
@@ -747,3 +775,41 @@ int lyapunov_free_tmps(
 
   return(0);
 }
+
+// save state of lyapunov computation
+int lyapunov_save_state(
+  double* flow,
+  _Complex double* u,
+  double* lyapunov_avg,
+  double prevtime,
+  double lyapunov_startingtime,
+  FILE* savefile,
+  int K1,
+  int K2,
+  const char* cmd_string,
+  const char* params_string,
+  const char* savefile_string,
+  const char* utfile_string,
+  FILE* utfile,
+  unsigned int command,
+  unsigned int algorithm,
+  double step,
+  double time,
+  unsigned int nthreads
+){
+  // save u and step
+  save_state(u, savefile, K1, K2, cmd_string, params_string, savefile_string, utfile_string, utfile, command, algorithm, step, time, nthreads);
+
+  if(savefile!=stderr && savefile!=stdout){
+    // save tangent flow
+    write_mat2_bin(flow,K1,K2,savefile);
+    // save time of QR decomposition
+    fwrite(&prevtime, sizeof(double), 1, savefile);
+    // save time at which the lyapunov computation started
+    fwrite(&lyapunov_startingtime, sizeof(double), 1, savefile);
+    // save lyapunov averages
+    write_vec2_bin(lyapunov_avg,K1,K2,savefile);
+  }
+
+  return 0;
+}

src/lyapunov.h

@@ -20,7 +20,7 @@ limitations under the License.
 #include "navier-stokes.h"
 
 // compute Lyapunov exponents
-int lyapunov( int K1, int K2, int N1, int N2, double final_time, double lyapunov_reset, unsigned int lyapunov_trigger, 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, bool keep_en_cst, double target_en, unsigned int algorithm, unsigned int algorithm_lyapunov, double starting_time, unsigned int nthreads);
+int lyapunov( int K1, int K2, int N1, int N2, double final_time, double lyapunov_reset, unsigned int lyapunov_trigger, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_cost, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, unsigned int algorithm_lyapunov, double starting_time, unsigned int nthreads, double* flow0, double* lyapunov_avg0, double prevtime, double lyapunov_startingtime, FILE* savefile, FILE* utfile, const char* cmd_string, const char* params_string, const char* savefile_string, const char* utfile_string);
 
 // comparison function for qsort
 int compare_double(const void* x, const void* y);
@@ -56,6 +56,9 @@ int lyapunov_init_tmps( double ** lyapunov, double ** lyapunov_avg, double ** fl
 // release vectors
 int lyapunov_free_tmps( double * lyapunov, double * lyapunov_avg, double * flow, _Complex double * u, double * tmp1, double * tmp2, double * tmp3, _Complex double * tmp4, _Complex double * tmp5, _Complex double * tmp6, _Complex double * tmp7, _Complex double * tmp8, _Complex double * tmp9, _Complex double * tmp10, double * tmp11, fft_vect fftu1, fft_vect fftu2, fft_vect fftu3, fft_vect fftu4, fft_vect fft1, fft_vect ifft, unsigned int algorithm, unsigned int algorithm_lyapunov);
 
+// save state of lyapunov computation
+int lyapunov_save_state( double* flow, _Complex double* u, double* lyapunov_avg, double prevtime, double lyapunov_startingtime, FILE* savefile, int K1, int K2, const char* cmd_string, const char* params_string, const char* savefile_string, const char* utfile_string, FILE* utfile, unsigned int command, unsigned int algorithm, double step, double time, unsigned int nthreads);
+
 #endif
 
 

src/main.c

@@ -28,6 +28,7 @@ limitations under the License.
 #include "dstring.h"
 #include "init.h"
 #include "int_tools.h"
+#include "io.h"
 #include "lyapunov.h"
 #include "navier-stokes.h"
 
@@ -62,7 +63,6 @@ typedef struct nstrophy_parameters {
   FILE* drivingfile;
   double lyapunov_reset;
   unsigned int lyapunov_trigger;
-  double D_epsilon;
   bool print_alpha;
 } nstrophy_parameters;
 
@@ -83,6 +83,8 @@ int args_from_file(dstring* params, unsigned int* command, unsigned int* nthread
 _Complex double* set_driving(nstrophy_parameters parameters);
 // set initial condition
 _Complex double* set_init(nstrophy_parameters* parameters);
+// set initial tangent flow for lyapunov exponents
+int set_lyapunov_flow_init( double** lyapunov_flow0, double** lyapunov_avg0, double* lyapunov_prevtime, double* lyapunov_startingtime, bool fromfile, nstrophy_parameters parameters);
 
 // signal handler
 void sig_handler (int signo);
@@ -116,6 +118,10 @@ int main (
   unsigned int nthreads=1;
   _Complex double* u0;
   _Complex double *g;
+  double* lyapunov_flow0;
+  double* lyapunov_avg0;
+  double lyapunov_prevtime;
+  double lyapunov_startingtime;
   dstring savefile_str;
   dstring utfile_str;
   dstring initfile_str;
@@ -123,6 +129,7 @@ int main (
   dstring resumefile_str;
   FILE* savefile=NULL;
   FILE* utfile=NULL;
+  bool resume=false;
 
   command=0;
 
@@ -162,6 +169,8 @@ int main (
 
   // if command is 'resume', then read args from file
   if(command==COMMAND_RESUME){
+    // remember that the original command was resume (to set values from init file)
+    resume=true;
     ret=args_from_file(&param_str, &command, &nthreads, &savefile_str, &utfile_str, dstring_to_str_noinit(&resumefile_str));
     if(ret<0){
       dstring_free(param_str);
@@ -225,6 +234,20 @@ int main (
   g=set_driving(parameters);
   // set initial condition
   u0=set_init(&parameters);
+  // read extra values from init file when resuming a computation
+  if(resume){
+    // read start time
+    fread(&(parameters.starting_time), sizeof(double), 1, parameters.initfile);
+    // if adaptive step algorithm
+    if(parameters.algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
+      // read delta
+      fread(&(parameters.delta), sizeof(double), 1, parameters.initfile);
+    }
+  }
+  // set initial condition for the lyapunov flow
+  if (command==COMMAND_LYAPUNOV){
+    set_lyapunov_flow_init(&lyapunov_flow0, &lyapunov_avg0, &lyapunov_prevtime, &lyapunov_startingtime, resume, parameters);
+  }
 
   // if init_enstrophy is not set in the parameters, then compute it from the initial condition
   if (parameters.init_enstrophy_or_energy!=FIX_ENSTROPHY){
@@ -256,6 +279,10 @@ int main (
       dstring_free(drivingfile_str);
       free(g);
       free(u0);
+      if (command==COMMAND_LYAPUNOV){
+	free(lyapunov_flow0);
+	free(lyapunov_avg0);
+      }
       return(-1);
     }
   }
@@ -272,6 +299,10 @@ int main (
       dstring_free(drivingfile_str);
       free(g);
       free(u0);
+      if (command==COMMAND_LYAPUNOV){
+	free(lyapunov_flow0);
+	free(lyapunov_avg0);
+      }
       return(-1);
     }
   }
@@ -286,7 +317,7 @@ 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, parameters.max_delta, parameters.adaptive_cost, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_enstrophy, 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_cost, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_enstrophy, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile, utfile, (char*)argv[0], dstring_to_str_noinit(&param_str), dstring_to_str_noinit(&savefile_str), dstring_to_str_noinit(&utfile_str));
   }
   else if(command==COMMAND_ENSTROPHY){
     // register signal handler to handle aborts
@@ -294,11 +325,17 @@ int main (
     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, parameters.max_delta, parameters.adaptive_cost, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_enstrophy, parameters.algorithm, parameters.print_freq, parameters.starting_time, parameters.print_alpha, nthreads, savefile, utfile, (char*)argv[0], dstring_to_str_noinit(&param_str), dstring_to_str_noinit(&savefile_str), dstring_to_str_noinit(&utfile_str));
   }
+  else if(command==COMMAND_ENSTROPHY_PRINT_INIT){
+    enstrophy_print_init(parameters.K1, parameters.K2, parameters.L, u0, g);
+  }
   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.max_delta, parameters.adaptive_cost, parameters.starting_time, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_enstrophy, parameters.algorithm, nthreads, savefile);
   }
   else if(command==COMMAND_LYAPUNOV){
-    lyapunov(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.lyapunov_reset, parameters.lyapunov_trigger, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_cost, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_enstrophy, parameters.algorithm, parameters.algorithm_lyapunov, parameters.starting_time, nthreads);
+    // register signal handler to handle aborts
+    signal(SIGINT, sig_handler);
+    signal(SIGTERM, sig_handler);
+    lyapunov(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.lyapunov_reset, parameters.lyapunov_trigger, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_cost, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_enstrophy, parameters.algorithm, parameters.algorithm_lyapunov, parameters.starting_time, nthreads, lyapunov_flow0, lyapunov_avg0, lyapunov_prevtime, lyapunov_startingtime, savefile, utfile, (char*)argv[0], dstring_to_str_noinit(&param_str), dstring_to_str_noinit(&savefile_str), dstring_to_str_noinit(&utfile_str));
   }
   else if(command==0){
     fprintf(stderr, "error: no command specified\n");
@@ -307,6 +344,10 @@ int main (
 
   free(g);
   free(u0);
+  if (command==COMMAND_LYAPUNOV){
+    free(lyapunov_flow0);
+    free(lyapunov_avg0);
+  }
 
   // free strings
   dstring_free(param_str);
@@ -523,6 +564,9 @@ int read_args(
       else if(strcmp(argv[i], "enstrophy")==0){
 	*command=COMMAND_ENSTROPHY;
       }
+      else if(strcmp(argv[i], "enstrophy_print_init")==0){
+	*command=COMMAND_ENSTROPHY_PRINT_INIT;
+      }
       else if(strcmp(argv[i], "quiet")==0){
 	*command=COMMAND_QUIET;
       }
@@ -541,7 +585,7 @@ int read_args(
   }
 
   // check that if the command is 'resume', then resumefile has been set
-  if(*command==COMMAND_RESUME && resumefile_str->length==0){
+  if(*command==COMMAND_RESUME && (resumefile_str==NULL || resumefile_str->length==0)){
     fprintf(stderr, "error: 'resume' command used, but no resume file\n");
     print_usage();
     return(-1);
@@ -852,13 +896,6 @@ int set_parameter(
       return(-1);
     }
   }
-  else if (strcmp(lhs,"D_epsilon")==0){
-    ret=sscanf(rhs,"%lf",&(parameters->D_epsilon));
-    if(ret!=1){
-      fprintf(stderr, "error: parameter 'D_epsilon' should be a double\n       got '%s'\n",rhs);
-      return(-1);
-    }
-  }
   else if (strcmp(lhs,"driving")==0){
     if (strcmp(rhs,"zero")==0){
       parameters->driving=DRIVING_ZERO;
@@ -1100,12 +1137,6 @@ _Complex double* set_init(
 
   case INIT_FILE:
     init_file_bin(u0, parameters->K1, parameters->K2, parameters->initfile);
-    // read start time
-    fread(&(parameters->starting_time), sizeof(double), 1, parameters->initfile);
-    if(parameters->algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
-      // read delta
-      fread(&(parameters->delta), sizeof(double), 1, parameters->initfile);
-    }
     break;
 
   case INIT_FILE_TXT:
@@ -1140,3 +1171,49 @@ _Complex double* set_init(
   
   return u0;
 }
+
+// set initial tangent flow for lyapunov exponents
+int set_lyapunov_flow_init(
+  double** lyapunov_flow0,
+  double** lyapunov_avg0,
+  double* lyapunov_prevtime,
+  double* lyapunov_startingtime,
+  bool fromfile, // whether to initialize from file
+  nstrophy_parameters parameters
+){
+  *lyapunov_flow0=calloc(4*(parameters.K1*(2*parameters.K2+1)+parameters.K2)*(parameters.K1*(2*parameters.K2+1)+parameters.K2),sizeof(double));
+  *lyapunov_avg0=calloc(2*(parameters.K1*(2*parameters.K2+1)+parameters.K2),sizeof(double));
+
+  // read from file or init from identity matrix
+  if(fromfile){
+    // read flow
+    read_mat2_bin(*lyapunov_flow0, parameters.K1, parameters.K2, parameters.initfile);
+    // read time of last QR decomposition
+    fread(lyapunov_prevtime, sizeof(double), 1, parameters.initfile);
+    // read time at which lyapunov computation was started
+    fread(lyapunov_startingtime, sizeof(double), 1, parameters.initfile);
+    // read averages
+    read_vec2_bin(*lyapunov_avg0, parameters.K1, parameters.K2, parameters.initfile);
+  } else {
+    // init with identity
+    int i,j;
+    for (i=0;i<2*(parameters.K1*(2*parameters.K2+1)+parameters.K2);i++){
+      for (j=0;j<2*(parameters.K1*(2*parameters.K2+1)+parameters.K2);j++){
+	if(i!=j){
+	  (*lyapunov_flow0)[i*2*(parameters.K1*(2*parameters.K2+1)+parameters.K2)+j]=0.;
+	} else {
+	  (*lyapunov_flow0)[i*2*(parameters.K1*(2*parameters.K2+1)+parameters.K2)+j]=1.;
+	}
+      }
+    }
+    // init prevtime
+    *lyapunov_prevtime=parameters.starting_time;
+    // init starting_time
+    *lyapunov_startingtime=parameters.starting_time;
+    // init averages
+    for(i=0;i<2*(parameters.K1*(2*parameters.K2+1)+parameters.K2);i++){
+      (*lyapunov_avg0)[i]=0;
+    }
+  }
+  return 0;
+}

src/navier-stokes.c

@@ -23,6 +23,8 @@ limitations under the License.
 #include <stdlib.h>
 #include <string.h>
 
+#define USIZE (K1*(2*K2+1)+K2)
+
 // compute solution as a function of time
 int uk(
   int K1,
@@ -46,7 +48,13 @@ int uk(
   double print_freq,
   double starting_time,
   unsigned int nthreads,
-  FILE* savefile
+  FILE* savefile,
+  FILE* utfile,
+  // for interrupt recovery
+  const char* cmd_string,
+  const char* params_string,
+  const char* savefile_string,
+  const char* utfile_string
 ){
   _Complex double* u;
   _Complex double* tmp1;
@@ -93,7 +101,6 @@ int uk(
     ns_step(algorithm, u, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_cost, L, g, time, starting_time, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en);
 
     time+=step;
-    step=next_step;
     
     if(time>(n+1)*print_freq){
       n++;
@@ -113,11 +120,25 @@ int uk(
       fprintf(stderr,"\n");
       printf("\n");
     }
+
+    // get ready for next step
+    step=next_step;
+
+    // catch abort signal
+    if (g_abort){
+      break;
+    }
   }
 
   // TODO: update handling of savefile
-  // save final entry to savefile
-  write_vec_bin(u, K1, K2, savefile);
+  if(savefile!=NULL){
+    save_state(u, savefile, K1, K2, cmd_string, params_string, savefile_string, utfile_string, utfile, COMMAND_UK, algorithm, step, time, nthreads);
+  }
+
+  // save final u to utfile in txt format
+  if(utfile!=NULL){
+    write_vec(u, K1, K2, utfile);
+  }
 
   ns_free_tmps(u, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, fft1, fft2, ifft, algorithm);
   return(0);
@@ -223,11 +244,19 @@ int enstrophy(
       // print to stderr so user can follow along
       if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
 	fprintf(stderr,"% .8e % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",time, avg_a, avg_energy, avg_enstrophy, alpha, energy, enstrophy, step);
+	if(K1>=1 && K2>=2){
+	  printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_energy, avg_enstrophy, alpha, energy, enstrophy, step, __real__ u[klookup_sym(1,1,K2)], __real__ u[klookup_sym(1,2,K2)]);
+	}else{
 	  printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_energy, avg_enstrophy, alpha, energy, enstrophy, step);
+	}
       } else {
 	fprintf(stderr,"% .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",time, avg_a, avg_energy, avg_enstrophy, alpha, energy, enstrophy);
+	if(K1>=1 && K2>=2){
+	  printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_energy, avg_enstrophy, alpha, energy, enstrophy, __real__ u[klookup_sym(1,1,K2)], __real__ u[klookup_sym(1,2,K2)]);
+	}else{
 	  printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_energy, avg_enstrophy, alpha, energy, enstrophy);
 	}
+      }
 
       // reset averages
       avg_a=0;
@@ -268,6 +297,25 @@ int enstrophy(
   return(0);
 }
 
+// compute enstrophy, alpha for the initial condition (useful for debugging)
+int enstrophy_print_init(
+  int K1,
+  int K2,
+  double L,
+  _Complex double* u0,
+  _Complex double* g
+){
+  double alpha, enstrophy, energy;
+
+  alpha=compute_alpha(u0, K1, K2, g, L);
+  enstrophy=compute_enstrophy(u0, K1, K2, L);
+  energy=compute_energy(u0, K1, K2);
+
+  // print
+  printf("% .15e % .15e % .15e\n", alpha, energy, enstrophy);
+  return(0);
+}
+
 // compute solution as a function of time, but do not print anything (useful for debugging)
 int quiet(
   int K1,
@@ -353,30 +401,30 @@ int ns_init_tmps(
   unsigned int algorithm
 ){
   // velocity field
-  *u=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
+  *u=calloc(USIZE,sizeof(_Complex double));
 
   // allocate tmp vectors for computation
   if(algorithm==ALGORITHM_RK2){
-    *tmp1=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp2=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
+    *tmp1=calloc(USIZE,sizeof(_Complex double));
+    *tmp2=calloc(USIZE,sizeof(_Complex double));
   } else if (algorithm==ALGORITHM_RK4){
-    *tmp1=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp2=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp3=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
+    *tmp1=calloc(USIZE,sizeof(_Complex double));
+    *tmp2=calloc(USIZE,sizeof(_Complex double));
+    *tmp3=calloc(USIZE,sizeof(_Complex double));
   } else if (algorithm==ALGORITHM_RKF45 || algorithm==ALGORITHM_RKDP54){
-    *tmp1=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp2=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp3=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp4=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp5=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp6=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp7=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
+    *tmp1=calloc(USIZE,sizeof(_Complex double));
+    *tmp2=calloc(USIZE,sizeof(_Complex double));
+    *tmp3=calloc(USIZE,sizeof(_Complex double));
+    *tmp4=calloc(USIZE,sizeof(_Complex double));
+    *tmp5=calloc(USIZE,sizeof(_Complex double));
+    *tmp6=calloc(USIZE,sizeof(_Complex double));
+    *tmp7=calloc(USIZE,sizeof(_Complex double));
   } else if (algorithm==ALGORITHM_RKBS32){
-    *tmp1=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp2=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp3=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp4=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp5=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
+    *tmp1=calloc(USIZE,sizeof(_Complex double));
+    *tmp2=calloc(USIZE,sizeof(_Complex double));
+    *tmp3=calloc(USIZE,sizeof(_Complex double));
+    *tmp4=calloc(USIZE,sizeof(_Complex double));
+    *tmp5=calloc(USIZE,sizeof(_Complex double));
   } else {
     fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
   };
@@ -462,7 +510,7 @@ int copy_u(
 ){
   int i;
 
-  for(i=0;i<K1*(2*K2+1)+K2;i++){
+  for(i=0;i<USIZE;i++){
     u[i]=u0[i];
   }
 
@@ -524,7 +572,6 @@ int ns_step(
   return (0);
 }
 
-// TODO: do not need to use klookup in any of the rk routines
 // RK 4 algorithm
 int ns_step_rk4(
   _Complex double* u,
@@ -546,69 +593,53 @@ int ns_step_rk4(
   bool keep_en_cst,
   double target_en
 ){
-  int kx,ky;
+  int i;
 
   // k1
   ns_rhs(tmp1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
   // add to output
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp3[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta/6*tmp1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp3[i]=u[i]+delta/6*tmp1[i];
   }
 
   // u+h*k1/2
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp2[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta/2*tmp1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp2[i]=u[i]+delta/2*tmp1[i];
   }
   // k2
   ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
   // add to output
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp3[klookup_sym(kx,ky,K2)]+=delta/3*tmp1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp3[i]+=delta/3*tmp1[i];
   }
 
   // u+h*k2/2
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp2[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta/2*tmp1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp2[i]=u[i]+delta/2*tmp1[i];
   }
   // k3
   ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
   // add to output
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp3[klookup_sym(kx,ky,K2)]+=delta/3*tmp1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp3[i]+=delta/3*tmp1[i];
   }
 
   // u+h*k3
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp2[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta*tmp1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp2[i]=u[i]+delta*tmp1[i];
   }
   // k4
   ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
   // add to output
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      u[klookup_sym(kx,ky,K2)]=tmp3[klookup_sym(kx,ky,K2)]+delta/6*tmp1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    u[i]=tmp3[i]+delta/6*tmp1[i];
   }
 
   // keep enstrophy constant
   if(keep_en_cst){
     double en=compute_enstrophy(u, K1, K2, L);
-    for(kx=0;kx<=K1;kx++){
-      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-	u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en);
-      }
+    for(i=0;i<USIZE;i++){
+      u[i]*=sqrt(target_en/en);
     }
   }
 
@@ -635,33 +666,27 @@ int ns_step_rk2(
   bool keep_en_cst,
   double target_en
 ){
-  int kx,ky;
+  int i;
 
   // k1
   ns_rhs(tmp1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
 
   // u+h*k1/2
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp2[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta/2*tmp1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp2[i]=u[i]+delta/2*tmp1[i];
   }
   // k2
   ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
   // add to output
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      u[klookup_sym(kx,ky,K2)]+=delta*tmp1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    u[i]+=delta*tmp1[i];
   }
 
   // keep enstrophy constant
   if(keep_en_cst){
     double en=compute_enstrophy(u, K1, K2, L);
-    for(kx=0;kx<=K1;kx++){
-      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-	u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en);
-      }
+    for(i=0;i<USIZE;i++){
+      u[i]*=sqrt(target_en/en);
     }
   }
 
@@ -701,7 +726,7 @@ int ns_step_rkf45(
   // whether to compute k1 or leave it as is
   bool compute_k1
 ){
-  int kx,ky;
+  int i;
   double cost;
 
   // k1: u(t)
@@ -710,53 +735,41 @@ int ns_step_rkf45(
   }
 
   // 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)/4*k1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)/4*k1[i];
   }
   ns_rhs(k2, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
 
   // k3 : u(t+3/8*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./32*k1[klookup_sym(kx,ky,K2)]+9./32*k2[klookup_sym(kx,ky,K2)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(3./32*k1[i]+9./32*k2[i]);
   }
   ns_rhs(k3, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
 
   // k4 : u(t+12./13*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)*(1932./2197*k1[klookup_sym(kx,ky,K2)]-7200./2197*k2[klookup_sym(kx,ky,K2)]+7296./2197*k3[klookup_sym(kx,ky,K2)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(1932./2197*k1[i]-7200./2197*k2[i]+7296./2197*k3[i]);
   }
   ns_rhs(k4, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
 
   // k5 : u(t+1*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)*(439./216*k1[klookup_sym(kx,ky,K2)]-8*k2[klookup_sym(kx,ky,K2)]+3680./513*k3[klookup_sym(kx,ky,K2)]-845./4104*k4[klookup_sym(kx,ky,K2)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(439./216*k1[i]-8*k2[i]+3680./513*k3[i]-845./4104*k4[i]);
   }
   ns_rhs(k5, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
 
   // k6 : u(t+1./2*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)*(-8./27*k1[klookup_sym(kx,ky,K2)]+2*k2[klookup_sym(kx,ky,K2)]-3544./2565*k3[klookup_sym(kx,ky,K2)]+1859./4104*k4[klookup_sym(kx,ky,K2)]-11./40*k5[klookup_sym(kx,ky,K2)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(-8./27*k1[i]+2*k2[i]-3544./2565*k3[i]+1859./4104*k4[i]-11./40*k5[i]);
   }
   ns_rhs(k6, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
   
   // next step
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+  for(i=0;i<USIZE;i++){
     // u
-      tmp[klookup_sym(kx,ky,K2)]=u[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)]);
+    tmp[i]=u[i]+(*delta)*(25./216*k1[i]+1408./2565*k3[i]+2197./4104*k4[i]-1./5*k5[i]);
     // U: save to k6, which is no longer needed
-      k6[klookup_sym(kx,ky,K2)]=u[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)]);
-    }
+    k6[i]=u[i]+(*delta)*(16./135*k1[i]+6656./12825*k3[i]+28561./56430*k4[i]-9./50*k5[i]+2./55*k6[i]);
   }
 
   // cost function
@@ -765,10 +778,8 @@ int ns_step_rkf45(
   // compare relative error with tolerance
   if(cost<tolerance){
     // copy 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)];
-      }
+    for(i=0;i<USIZE;i++){
+      u[i]=tmp[i];
     }
     // next delta to use in future steps (do not exceed max_delta)
     *next_delta=fmin(max_delta, (*delta)*pow(tolerance/cost,0.2));
@@ -776,10 +787,8 @@ int ns_step_rkf45(
     // keep enstrophy constant
     if(keep_en_cst){
       double en=compute_enstrophy(u, K1, K2, L);
-      for(kx=0;kx<=K1;kx++){
-	for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-	  u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en);
-	}
+      for(i=0;i<USIZE;i++){
+	u[i]*=sqrt(target_en/en);
       }
     }
   }
@@ -826,7 +835,7 @@ int ns_step_rkbs32(
   // whether to compute k1
   bool compute_k1
 ){
-  int kx,ky;
+  int i;
   double cost;
 
   // k1: u(t)
@@ -836,36 +845,28 @@ int ns_step_rkbs32(
   }
 
   // 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)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)/2*(*k1)[i];
   }
   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)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(3./4*k2[i]);
   }
   ns_rhs(k3, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
 
   // k4 : u(t+delta)
   // tmp computed 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)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(2./9*(*k1)[i]+1./3*k2[i]+4./9*k3[i]);
   }
   ns_rhs(*k4, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
 
   // next step
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+  for(i=0;i<USIZE;i++){
     // U: store in k3, which is no longer needed
-      k3[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*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)]);
-    }
+    k3[i]=u[i]+(*delta)*(7./24*(*k1)[i]+1./4*k2[i]+1./3*k3[i]+1./8*(*k4)[i]);
   }
   
   // compute cost
@@ -874,10 +875,8 @@ int ns_step_rkbs32(
   // compare cost with tolerance
   if(cost<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)];
-      }
+    for(i=0;i<USIZE;i++){
+      u[i]=tmp[i];
     }
     // next delta to use in future steps (do not exceed max_delta)
     *next_delta=fmin(max_delta, (*delta)*pow(tolerance/cost,1./3));
@@ -890,10 +889,8 @@ int ns_step_rkbs32(
     // keep enstrophy constant
     if(keep_en_cst){
       double en=compute_enstrophy(u, K1, K2, L);
-      for(kx=0;kx<=K1;kx++){
-	for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-	  u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en);
-	}
+      for(i=0;i<USIZE;i++){
+	u[i]*=sqrt(target_en/en);
       }
     }
   }
@@ -941,7 +938,7 @@ int ns_step_rkdp54(
   // whether to compute k1
   bool compute_k1
 ){
-  int kx,ky;
+  int i;
   double cost;
 
   // k1: u(t)
@@ -951,61 +948,47 @@ int ns_step_rkdp54(
   }
 
   // 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)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)/5*(*k1)[i];
   }
   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)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(3./40*(*k1)[i]+9./40*(*k2)[i]);
   }
   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)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(44./45*(*k1)[i]-56./15*(*k2)[i]+32./9*k3[i]);
   }
   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)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(19372./6561*(*k1)[i]-25360./2187*(*k2)[i]+64448./6561*k3[i]-212./729*k4[i]);
   }
   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)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(9017./3168*(*k1)[i]-355./33*(*k2)[i]+46732./5247*k3[i]+49./176*k4[i]-5103./18656*k5[i]);
   }
   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++){
+  for(i=0;i<USIZE;i++){
     // 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)]);
-    }
+    tmp[i]=u[i]+(*delta)*(35./384*(*k1)[i]+500./1113*k3[i]+125./192*k4[i]-2187./6784*k5[i]+11./84*k6[i]);
   }
   // store in k2, which is not needed anymore
   ns_rhs(*k2, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
 
   //next step
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+  for(i=0;i<USIZE;i++){
     // U: store in k6, which is not needed anymore
-      k6[klookup_sym(kx,ky,K2)]=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)]);
-    }
+    k6[i]=u[i]+(*delta)*(5179./57600*(*k1)[i]+7571./16695*k3[i]+393./640*k4[i]-92097./339200*k5[i]+187./2100*k6[i]+1./40*(*k2)[i]);
   }
   
   // compute cost
@@ -1014,10 +997,8 @@ int ns_step_rkdp54(
   // compare relative error with tolerance
   if(cost<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)];
-      }
+    for(i=0;i<USIZE;i++){
+      u[i]=tmp[i];
     }
     // next delta to use in future steps (do not exceed max_delta)
     *next_delta=fmin(max_delta, (*delta)*pow(tolerance/cost,0.2));
@@ -1030,10 +1011,8 @@ int ns_step_rkdp54(
     // keep enstrophy constant
     if(keep_en_cst){
       double en=compute_enstrophy(u, K1, K2, L);
-      for(kx=0;kx<=K1;kx++){
-	for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-	  u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en);
-	}
+      for(i=0;i<USIZE;i++){
+	u[i]*=sqrt(target_en/en);
       }
     }
   }
@@ -1138,7 +1117,7 @@ int ns_rhs(
     alpha=compute_alpha(u,K1,K2,g,L);
   }
 
-  for(i=0; i<K1*(2*K2+1)+K2; i++){
+  for(i=0; i<USIZE; i++){
     out[i]=0;
   }
   for(kx=0;kx<=K1;kx++){

src/navier-stokes.h

@@ -31,10 +31,12 @@ 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, double max_delta, unsigned int adaptive_cost, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, 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_cost, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, double print_freq, double starting_time, unsigned int nthreadsl, FILE* savefile, FILE* utfile, const char* cmd_string, const char* params_string, const char* savefile_string, const char* utfile_string);
 
 // 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, double max_delta, unsigned int adaptive_cost, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, double print_freq, double starting_time, bool print_alpha, unsigned int nthreads, FILE* savefile, FILE* utfile, const char* cmd_string, const char* params_string, const char* savefile_string, const char* utfile_string);
+// compute enstrophy, alpha for the initial condition (useful for debugging)
+int enstrophy_print_init( int K1, int K2, double L, _Complex double* u0, _Complex double* g);
 
 // 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 max_delta, unsigned int adaptive_cost, double starting_time, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, unsigned int nthreads, FILE* savefile);