interrupting lyapunov mentioned in README

README.md

@@ -164,14 +164,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`
+The `enstrophy` and `lyapunov` computations can be interrupted by sending
-signal (e.g. by pressing `Ctrl-C`.) When Nstrophy receives the `SIGINT` signal,
+Nstrophy the `SIGINT` signal (e.g. by pressing `Ctrl-C`.) When Nstrophy
-it finishes its current step and writes the value of uk, either to `savefile`
+receives the `SIGINT` signal, it finishes its current step and writes the value
-if such a file was specified on the command line (using the `-s` flag), or to
+of uk, either to `savefile` if such a file was specified on the command line
-`stderr`. In addition, when a `savefile` is specified it writes the command
+(using the `-s` flag), or to `stderr`. In addition, when a `savefile` is
-that needs to be used to resume the computation (which essentially just sets
+specified it writes the command that needs to be used to resume the computation
-the appropriate `starting_time` and `init:file:<savefile>` parameters. The data
+(which essentially just sets the appropriate `starting_time` and
-written to the `savefile` is binary.
+`init:file:<savefile>` parameters). The data written to the `savefile` is
+binary.
 
 
 # License

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;
 }
 
@@ -275,6 +335,10 @@ int save_state(
     if(savefile==stderr || savefile==stdout){
       fprintf(savefile,";starting_time=%.15e", time);
     }
+    // instruction to read init flow from file if computation is lyapunov
+    if(command==COMMAND_LYAPUNOV){
+      fprintf(savefile,";init_flow=file");
+    }
     fprintf(savefile,"\"");
   }
 

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,16 @@ int lyapunov(
   unsigned int algorithm,
   unsigned int algorithm_lyapunov,
   double starting_time,
-  unsigned int nthreads
+  unsigned int nthreads,
+  double* flow0,
+  double* lyapunov_avg0,
+  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,15 +92,13 @@ 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]=flow0[i*MATSIZE+j];
-	flow[i*MATSIZE+j]=0.;
-      } else {
-	flow[i*MATSIZE+j]=1.;
-      }
     }
+    lyapunov_avg[i]=lyapunov_avg0[i];
   }
 
   // store step (useful for adaptive step methods)
@@ -168,7 +176,25 @@ 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, 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,35 @@ int lyapunov_free_tmps(
 
   return(0);
 }
+
+// save state of lyapunov computation
+int lyapunov_save_state(
+  double* flow,
+  _Complex double* u,
+  double* lyapunov_avg,
+  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 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, 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, 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,7 @@ typedef struct nstrophy_parameters {
   FILE* drivingfile;
   double lyapunov_reset;
   unsigned int lyapunov_trigger;
-  double D_epsilon;
+  bool init_flow_file;
   bool print_alpha;
 } nstrophy_parameters;
 
@@ -83,6 +84,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, nstrophy_parameters parameters);
 
 // signal handler
 void sig_handler (int signo);
@@ -116,6 +119,8 @@ int main (
   unsigned int nthreads=1;
   _Complex double* u0;
   _Complex double *g;
+  double* lyapunov_flow0;
+  double* lyapunov_avg0;
   dstring savefile_str;
   dstring utfile_str;
   dstring initfile_str;
@@ -225,6 +230,10 @@ int main (
   g=set_driving(parameters);
   // set initial condition
   u0=set_init(&parameters);
+  // set initial condition for the lyapunov flow
+  if (command==COMMAND_LYAPUNOV){
+    set_lyapunov_flow_init(&lyapunov_flow0, &lyapunov_avg0, 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 +265,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 +285,10 @@ int main (
       dstring_free(drivingfile_str);
       free(g);
       free(u0);
+      if (command==COMMAND_LYAPUNOV){
+	free(lyapunov_flow0);
+	free(lyapunov_avg0);
+      }
       return(-1);
     }
   }
@@ -298,7 +315,7 @@ int main (
     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);
+    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, 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 +324,10 @@ int main (
 
   free(g);
   free(u0);
+  if (command==COMMAND_LYAPUNOV){
+    free(lyapunov_flow0);
+    free(lyapunov_avg0);
+  }
 
   // free strings
   dstring_free(param_str);
@@ -541,7 +562,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);
@@ -581,6 +602,7 @@ int set_default_params(
   parameters->algorithm_lyapunov=ALGORITHM_RK4;
   // default lyapunov_reset will be print_time (set later) for now set target to 0 to indicate it hasn't been set explicitly
   parameters->lyapunov_trigger=0;
+  parameters->init_flow_file=false;
 
   parameters->print_alpha=false;
   
@@ -852,13 +874,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;
@@ -985,6 +1000,16 @@ int set_parameter(
       return(-1);
     }
   }
+  else if (strcmp(lhs,"init_flow")==0){
+    if(strcmp(rhs,"file")==0){
+      parameters->init_flow_file=true;
+    } else if (strcmp(rhs,"identity")==0){
+      parameters->init_flow_file=false;
+    } else {
+      fprintf(stderr, "error: parameter 'init_flow' should be 'file' or 'identity'\n       got '%s'\n",rhs);
+      return(-1);
+    }
+  }
   else{
     fprintf(stderr, "error: unrecognized parameter '%s'\n",lhs);
     return(-1);
@@ -1140,3 +1165,36 @@ _Complex double* set_init(
   
   return u0;
 }
+
+// set initial tangent flow for lyapunov exponents
+int set_lyapunov_flow_init(
+  double** lyapunov_flow0,
+  double** lyapunov_avg0,
+  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));
+
+  if(parameters.init_flow_file){
+    // read flow
+    read_mat2_bin(*lyapunov_flow0, parameters.K1, parameters.K2, 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.;
+	}
+      }
+    }
+    for(i=0;i<2*(parameters.K1*(2*parameters.K2+1)+parameters.K2);i++){
+      (*lyapunov_avg0)[i]=0;
+    }
+  }
+  return 0;
+}