Print negative values of alpha

src/io.c

@@ -17,6 +17,7 @@ limitations under the License.
 #include <stdbool.h>
 #include <stdlib.h>
 #include <string.h>
+#include "constants.cpp"
 #include "io.h"
 #include "navier-stokes.h"
 
@@ -70,7 +71,7 @@ int read_vec(_Complex double* out, int K1, int K2, FILE* file){
   }
 
   // allocate line buffer
-  line=calloc(sizeof(char), len);
+  line=calloc(len,sizeof(char));
 
   while(1){
     c=fgetc(file);
@@ -125,7 +126,7 @@ int read_vec(_Complex double* out, int K1, int K2, FILE* file){
       // check that there is room in buffer
       if(pos==len){
 	// too short: reallocate
-	line_realloc=calloc(sizeof(char), 2*len);
+	line_realloc=calloc(2*len,sizeof(char));
 	for(pos=0;pos<len;pos++){
 	  line_realloc[pos]=line[pos];
 	}
@@ -225,3 +226,119 @@ int remove_entry(
 
   return 0;
 }
+
+
+// save state to savefile
+int save_state(
+  _Complex double* u,
+  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
+){
+  if(savefile==NULL){
+    fprintf(stderr, "error while writing savefile: savefile not open\n");
+    return -1;
+  }
+
+  fprintf(savefile,"# Continue computation with\n");
+
+  // command to resume
+  fprintf(savefile,"#! ");
+  fprintf(savefile, cmd_string);
+  // params
+  // allocate buffer for params
+  if(params_string!=NULL) {
+    char* params=calloc(strlen(params_string)+1,sizeof(char));
+    strcpy(params, params_string);
+    remove_entry(params, "starting_time");
+    remove_entry(params, "init");
+    if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
+      remove_entry(params, "delta");
+    }
+    fprintf(savefile," -p \"%s;init=file:%s", params, savefile_string);
+    free(params);
+    // write delta if adaptive, and not writing binary
+    if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD && (savefile==stderr || savefile==stdout)){
+      fprintf(savefile,";delta=%.15e", step);
+    }
+    // write starting_time if not writing binary
+    if(savefile==stderr || savefile==stdout){
+      fprintf(savefile,";starting_time=%.15e", time);
+    }
+    fprintf(savefile,"\"");
+  }
+
+  // utfile
+  if(utfile!=NULL){
+    fprintf(savefile," -u \"%s\"", utfile_string);
+  }
+  // threads
+  if(nthreads!=1){
+    fprintf(savefile," -t %d", nthreads);
+  }
+
+  switch (command){
+  case COMMAND_UK:
+    fprintf(savefile," uk\n");
+    break;
+  case COMMAND_ENSTROPHY:
+    fprintf(savefile," enstrophy\n");
+    break;
+  case COMMAND_QUIET:
+    fprintf(savefile," quiet\n");
+    break;
+  case COMMAND_LYAPUNOV:
+    fprintf(savefile," lyapunov\n");
+    break;
+  case COMMAND_RESUME:
+    fprintf(savefile," resume\n");
+    break;
+
+  default:
+    fprintf(stderr,"error while writing savefile: unrecognized command %d\n", command);
+    return(-1);
+    break;
+  }
+
+  // save final u to savefile
+  if(savefile==stderr || savefile==stdout){
+    write_vec(u, K1, K2, savefile);
+  } else {
+    write_vec_bin(u, K1, K2, savefile);
+    // last binary entry: starting time
+    fwrite(&time, sizeof(double), 1, savefile);
+    // extra binary data for adaptive algorithm
+    if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
+      fwrite(&step, sizeof(double), 1, savefile);
+    }
+  }
+
+  return 0;
+}
+
+// write u to utfile in plain text
+int write_utfile(
+  _Complex double* u,
+  FILE* utfile,
+  int K1,
+  int K2
+){
+
+  if(utfile==NULL){
+    fprintf(stderr, "error while writing utfile: utfile is not open\n");
+    return -1;
+  }
+
+  write_vec(u, K1, K2, utfile);
+  return 0;
+}

src/io.h

@@ -30,4 +30,9 @@ int read_vec_bin(_Complex double* u, int K1, int K2, FILE* file);
 // remove an entry from params string (inplace)
 int remove_entry(char* param_str, char* entry);
 
+// save state to savefile
+int save_state(_Complex double* u, 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);
+
+// write u to utfile in plain text
+int write_utfile(_Complex double* u, FILE* utfile, int K1, int K2);
 #endif

src/lyapunov.c

@@ -301,14 +301,14 @@ int lyapunov_init_tmps(
 ){
   ns_init_tmps(u, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, fft1, fft2, ifft, K1, K2, N1, N2, nthreads, algorithm);
 
-  *lyapunov=calloc(sizeof(double),MATSIZE);
-  *lyapunov_avg=calloc(sizeof(double),MATSIZE);
-  *Du_prod=calloc(sizeof(double),MATSIZE*MATSIZE);
-  *Du=calloc(sizeof(double),MATSIZE*MATSIZE);
-  *prevu=calloc(sizeof(_Complex double),USIZE);
-  *tmp1=calloc(sizeof(_Complex double),USIZE);
-  *tmp2=calloc(sizeof(_Complex double),USIZE);
-  *tmp10=calloc(sizeof(double),MATSIZE*MATSIZE);
+  *lyapunov=calloc(MATSIZE,sizeof(double));
+  *lyapunov_avg=calloc(MATSIZE,sizeof(double));
+  *Du_prod=calloc(MATSIZE*MATSIZE,sizeof(double));
+  *Du=calloc(MATSIZE*MATSIZE,sizeof(double));
+  *prevu=calloc(USIZE,sizeof(_Complex double));
+  *tmp1=calloc(USIZE,sizeof(_Complex double));
+  *tmp2=calloc(USIZE,sizeof(_Complex double));
+  *tmp10=calloc(MATSIZE*MATSIZE,sizeof(double));
 
   return(0);
 }

src/main.c

@@ -958,7 +958,7 @@ int args_from_file(
 _Complex double* set_driving(
   nstrophy_parameters parameters
 ){
-  _Complex double* g=calloc(sizeof(_Complex double),parameters.K1*(2*parameters.K2+1)+parameters.K2);
+  _Complex double* g=calloc(parameters.K1*(2*parameters.K2+1)+parameters.K2,sizeof(_Complex double));
 
   switch(parameters.driving){
   case DRIVING_ZERO:
@@ -987,7 +987,7 @@ _Complex double* set_driving(
 _Complex double* set_init(
   nstrophy_parameters* parameters
 ){
-  _Complex double* u0=calloc(sizeof(_Complex double),parameters->K1*(2*parameters->K2+1)+parameters->K2);
+  _Complex double* u0=calloc(parameters->K1*(2*parameters->K2+1)+parameters->K2,sizeof(_Complex double));
 
   switch(parameters->init){
   case INIT_RANDOM:

src/navier-stokes.c

@@ -149,10 +149,10 @@ int enstrophy(
   FILE* savefile,
   FILE* utfile,
   // for interrupt recovery
-  char* cmd_string,
-  char* params_string,
-  char* savefile_string,
-  char* utfile_string
+  const char* cmd_string,
+  const char* params_string,
+  const char* savefile_string,
+  const char* utfile_string
 ){
   _Complex double* u;
   _Complex double* tmp1;
@@ -163,9 +163,9 @@ int enstrophy(
   _Complex double* tmp6;
   _Complex double* tmp7;
   double time;
-  double alpha, enstrophy;
+  double alpha, enstrophy, energy;
   double prevtime;
-  double avg_a,avg_en,avg_en_x_a;
+  double avg_a,avg_enstrophy,avg_energy;
   // index
   fft_vect fft1;
   fft_vect fft2;
@@ -181,8 +181,8 @@ int enstrophy(
 
   // init running average
   avg_a=0;
-  avg_en=0;
-  avg_en_x_a=0;
+  avg_enstrophy=0;
+  avg_energy=0;
   prevtime=starting_time;
 
   // period
@@ -198,6 +198,7 @@ int enstrophy(
 
     alpha=compute_alpha(u, K1, K2, g, L);
     enstrophy=compute_enstrophy(u, K1, K2, L);
+    energy=compute_energy(u, K1, K2);
 
     // check that anything is nan
     if(isnan(alpha) || isnan(enstrophy)){
@@ -207,33 +208,39 @@ int enstrophy(
 
     // add to running averages (estimate the total duration of interval as print_freq, will be adjusted later)
     avg_a+=alpha*(step/print_freq);
-    avg_en+=enstrophy*(step/print_freq);
-    avg_en_x_a+=enstrophy*alpha*(step/print_freq);
+    avg_enstrophy+=enstrophy*(step/print_freq);
+    avg_energy+=energy*(step/print_freq);
 
     if(time>(n+1)*print_freq){
       n++;
 
       // adjust duration of interval to its actual value
       avg_a*=print_freq/(time-prevtime);
-      avg_en*=print_freq/(time-prevtime);
-      avg_en_x_a*=print_freq/(time-prevtime);
+      avg_enstrophy*=print_freq/(time-prevtime);
+      avg_energy*=print_freq/(time-prevtime);
 
       // 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_en, avg_en_x_a, alpha, enstrophy, alpha*enstrophy, step);
-	printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_en, avg_en_x_a, alpha, enstrophy, alpha*enstrophy, step);
+	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);
       } else {
-	fprintf(stderr,"% .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",time, avg_a, avg_en, avg_en_x_a, alpha, enstrophy, alpha*enstrophy);
-	printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_en, avg_en_x_a, alpha, enstrophy, alpha*enstrophy);
+	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);
       }
 
       // reset averages
       avg_a=0;
-      avg_en=0;
-      avg_en_x_a=0;
+      avg_enstrophy=0;
+      avg_energy=0;
       prevtime=time;
     }
 
+    // print alpha when it gets negative
+    if(alpha<0){
+      fprintf(stderr,"## negative alpha: % .8e at time % .8e\n",alpha, time);
+      printf("## negative alpha: % .15e at time % .15e\n",alpha, time);
+    }
+
     // get ready for next step
     step=next_step;
 
@@ -248,57 +255,7 @@ int enstrophy(
   }
 
   if(savefile!=NULL){
-    fprintf(savefile,"# Continue computation with\n");
-
-    // command to resume
-    fprintf(savefile,"#! ");
-    fprintf(savefile, cmd_string);
-    // params
-    // allocate buffer for params
-    if(params_string!=NULL) {
-      char* params=calloc(sizeof(char), strlen(params_string)+1);
-      strcpy(params, params_string);
-      remove_entry(params, "starting_time");
-      remove_entry(params, "init");
-      if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
-	remove_entry(params, "delta");
-      }
-      fprintf(savefile," -p \"%s;init=file:%s", params, savefile_string);
-      free(params);
-      // write delta if adaptive, and not writing binary
-      if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD && (savefile==stderr || savefile==stdout)){
-	fprintf(savefile,";delta=%.15e", step);
-      }
-      // write starting_time if not writing binary
-      if(savefile==stderr || savefile==stdout){
-	fprintf(savefile,";starting_time=%.15e", time);
-      }
-      fprintf(savefile,"\"");
-    }
-
-    // utfile
-    if(utfile!=NULL){
-      fprintf(savefile," -u \"%s\"", utfile_string);
-    }
-    // threads
-    if(nthreads!=1){
-      fprintf(savefile," -t %d", nthreads);
-    }
-
-    fprintf(savefile," enstrophy\n");
-
-    // save final u to savefile
-    if(savefile==stderr || savefile==stdout){
-      write_vec(u, K1, K2, savefile);
-    } else {
-      write_vec_bin(u, K1, K2, savefile);
-      // last binary entry: starting time
-      fwrite(&time, sizeof(double), 1, savefile);
-      // extra binary data for adaptive algorithm
-      if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
-	fwrite(&step, sizeof(double), 1, savefile);
-      }
-    }
+    save_state(u, savefile, K1, K2, cmd_string, params_string, savefile_string, utfile_string, utfile, COMMAND_ENSTROPHY, algorithm, step, time, nthreads);
   }
 
   // save final u to utfile in txt format
@@ -395,30 +352,30 @@ int ns_init_tmps(
   unsigned int algorithm
 ){
   // velocity field
-  *u=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
+  *u=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
 
   // allocate tmp vectors for computation
   if(algorithm==ALGORITHM_RK2){
-    *tmp1=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
-    *tmp2=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
+    *tmp1=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
+    *tmp2=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
   } else if (algorithm==ALGORITHM_RK4){
-    *tmp1=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
-    *tmp2=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
-    *tmp3=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
+    *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));
   } else if (algorithm==ALGORITHM_RKF45 || algorithm==ALGORITHM_RKDP54){
-    *tmp1=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
-    *tmp2=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
-    *tmp3=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
-    *tmp4=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
-    *tmp5=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
-    *tmp6=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
-    *tmp7=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
+    *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));
   } else if (algorithm==ALGORITHM_RKBS32){
-    *tmp1=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
-    *tmp2=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
-    *tmp3=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
-    *tmp4=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
-    *tmp5=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
+    *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));
   } else {
     fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
   };

src/navier-stokes.h

@@ -34,7 +34,7 @@ typedef struct fft_vects {
 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, bool keep_en_cst, double target_en, 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, 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, double print_freq, double starting_time, unsigned int nthreads, FILE* savefile, FILE* utfile, char* cmd_string, char* params_string, char* savefile_string, char* utfile_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, bool keep_en_cst, double target_en, unsigned int algorithm, double print_freq, double starting_time, unsigned int nthreads, FILE* savefile, FILE* utfile, const char* cmd_string, const char* params_string, const char* savefile_string, const char* utfile_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 max_delta, unsigned int adaptive_norm, 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);