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`.
 
@@ -161,11 +170,6 @@ should be a `;` sperated list of `key=value` pairs. The possible keys are
   be `RK4` for Runge-Kutta 4 (default) or `RK2` for Runge-Kutta 2. Adaptive
   step algorithms cannot be used for the tangent flow.
 
-* `init_flow` (`identity` or `file` (default)): if set to `file`, then read the
-  initial condition for the tangent flow (used for the Lyapunov exponent
-  computation) from the init file (the same as for `init`, which needs to be
-  specified). Otherwise, the flow is initialized as the identity matrix.
-
 
 # Interrupting and resuming the computation
 

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

@@ -258,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++){
-    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'){
-	go=0;
-	break;
+    // 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';
@@ -322,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");
     }
@@ -335,10 +352,6 @@ 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/main.c

@@ -63,7 +63,6 @@ typedef struct nstrophy_parameters {
   FILE* drivingfile;
   double lyapunov_reset;
   unsigned int lyapunov_trigger;
-  bool init_flow_file;
   bool print_alpha;
 } nstrophy_parameters;
 
@@ -85,7 +84,7 @@ _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, nstrophy_parameters parameters);
+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);
@@ -130,6 +129,7 @@ int main (
   dstring resumefile_str;
   FILE* savefile=NULL;
   FILE* utfile=NULL;
+  bool resume=false;
 
   command=0;
 
@@ -169,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);
@@ -232,9 +234,19 @@ 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, parameters);
+    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
@@ -305,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
@@ -313,6 +325,9 @@ 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);
   }
@@ -549,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;
       }
@@ -607,7 +625,6 @@ 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;
   
@@ -686,12 +703,6 @@ int read_params(
     parameters->lyapunov_reset=parameters->print_freq;
   }
 
-  // check that if flow_init is used, then so is init
-  if(parameters->init_flow_file && parameters->init!=INIT_FILE){
-    fprintf(stderr, "error: cannot use 'init_flow:file' if 'init' is not a binary file\n");
-    return(-1);
-  }
-
   return(0);
 }
 
@@ -1011,16 +1022,6 @@ 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);
@@ -1136,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:
@@ -1183,12 +1178,14 @@ int set_lyapunov_flow_init(
   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));
 
-  if(parameters.init_flow_file){
+  // 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

src/navier-stokes.c

@@ -48,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;
@@ -95,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++;
@@ -115,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);
@@ -225,10 +244,18 @@ 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);
-	printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e\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);
-	printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e\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
@@ -270,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,

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);