Remove avg_window option: always use print_freq

1 files changed, 18 insertions, 103 deletions

diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 40fb9aa..b9a3979 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -91,7 +91,6 @@ int eea(
   _Complex double* g,
   bool irreversible,
   unsigned int print_freq,
-  unsigned int running_avg_window,
   unsigned int starting_time,
   unsigned int nthreads,
   FILE* savefile,
@@ -106,11 +105,7 @@ int eea(
   _Complex double* tmp3;
   double alpha, energy, enstrophy;
   double avg_e,avg_a,avg_en;
-  // quantities needed to compute averages
-  double *save_print_e, *save_print_a, *save_print_en;
-  double *save_print_short_e, *save_print_short_a, *save_print_short_en;
   // index
-  unsigned int i;
   unsigned int t;
   fft_vect fft1;
   fft_vect fft2;
@@ -126,27 +121,8 @@ int eea(
   avg_a=0;
   avg_en=0;
 
-  // need this for running average
-  unsigned int short_len=running_avg_window % print_freq;
-  // number of full print intervals in running average window
-  unsigned int nr_print_in_window = (running_avg_window-short_len)/print_freq;
-  double avg_print_e=0;
-  double avg_print_a=0;
-  double avg_print_en=0;
-  double avg_print_short_e=0;
-  double avg_print_short_a=0;
-  double avg_print_short_en=0;
-
-  // save e a en for running average
-  // only need to save if window is non zero
-  if(running_avg_window!=0){
-    save_print_e=calloc(sizeof(double), nr_print_in_window);
-    save_print_a=calloc(sizeof(double), nr_print_in_window);
-    save_print_en=calloc(sizeof(double), nr_print_in_window);
-    save_print_short_e=calloc(sizeof(double), nr_print_in_window+1);
-    save_print_short_a=calloc(sizeof(double), nr_print_in_window+1);
-    save_print_short_en=calloc(sizeof(double), nr_print_in_window+1);
-  }
+  // special first case when starting_time is not a multiple of print_freq
+  unsigned int first_box = print_freq - (starting_time % print_freq);
 
   // iterate
   for(t=starting_time;t<starting_time+nsteps;t++){
@@ -160,78 +136,26 @@ int eea(
     enstrophy=compute_enstrophy(u, K1, K2, L);
     
     // running average
-    // if window=0, then take average over all times
-    if(running_avg_window==0){
-      if(t!=0){
-	avg_e=avg_e-(avg_e-energy)/t;
-	avg_a=avg_a-(avg_a-alpha)/t;
-	avg_en=avg_en-(avg_en-enstrophy)/t;
-      }
+    // reset averages
+    if(t % print_freq == 1){
+      avg_e=0;
+      avg_a=0;
+      avg_en=0;
     }
-    else{
-      // reset print interval averages
-      if(t % print_freq == 1){
-	avg_print_e=0;
-	avg_print_a=0;
-	avg_print_en=0;
-	avg_print_short_e=0;
-	avg_print_short_a=0;
-	avg_print_short_en=0;
-      }
-
-      // compute average over print interval
-      avg_print_e=avg_print_e+energy/print_freq;
-      avg_print_a=avg_print_a+alpha/print_freq;
-      avg_print_en=avg_print_en+enstrophy/print_freq;
 
-      // compute average over the last short_len elements of a print interval
-      if (short_len != 0 && (t % print_freq > print_freq-short_len || t % print_freq == 0)){
-	avg_print_short_e+=energy/short_len;
-	avg_print_short_a+=alpha/short_len;
-	avg_print_short_en+=enstrophy/short_len;
-      }
-
-      if(t % print_freq ==0 && t>starting_time){
-	// compute averages
-	if (t > running_avg_window + starting_time) {
-	  avg_e=save_print_short_e[nr_print_in_window]*((double)short_len/running_avg_window);
-	  avg_a=save_print_short_a[nr_print_in_window]*((double)short_len/running_avg_window);
-	  avg_en=save_print_short_en[nr_print_in_window]*((double)short_len/running_avg_window);
-	  for(i=0;i<nr_print_in_window;i++){
-	    // prevent warnings: these will be initialized as long as t > running_avg_window
-	    #pragma GCC diagnostic push
-	    #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
-	    avg_e+=save_print_e[i]*((double)print_freq/running_avg_window);
-	    avg_a+=save_print_a[i]*((double)print_freq/running_avg_window);
-	    avg_en+=save_print_en[i]*((double)print_freq/running_avg_window);
-	    #pragma GCC diagnostic pop
-	  }
-	}
-
-	// memorize averages for future use
-	// need to keep track of case ==0, otherwise nr_print_in_window-1 is not an unsigned int
-	if(nr_print_in_window>0){
-	  for(i=0;i<nr_print_in_window-1;i++){
-	    save_print_e[i+1]=save_print_e[i];
-	    save_print_a[i+1]=save_print_a[i];
-	    save_print_en[i+1]=save_print_en[i];
-	  }
-	}
-	for(i=0;i<nr_print_in_window;i++){
-	  save_print_short_e[i+1]=save_print_short_e[i];
-	  save_print_short_a[i+1]=save_print_short_a[i];
-	  save_print_short_en[i+1]=save_print_short_en[i];
-	}
-	save_print_e[0]=avg_print_e;
-	save_print_a[0]=avg_print_a;
-	save_print_en[0]=avg_print_en;
-	save_print_short_e[0]=avg_print_short_e;
-	save_print_short_a[0]=avg_print_short_a;
-	save_print_short_en[0]=avg_print_short_en;
-      }
+    // compute average
+    // different computationin first box if starting_time is not a multiple of print_freq
+    if(t < starting_time + first_box){
+      avg_e+=energy/first_box;
+      avg_a+=alpha/first_box;
+      avg_en+=enstrophy/first_box;
+    } else {
+      avg_e+=energy/print_freq;
+      avg_a+=alpha/print_freq;
+      avg_en+=enstrophy/print_freq;
     }
 
-    if(t>running_avg_window + starting_time && t%print_freq==0){
+    if(t>starting_time && t%print_freq==0){
       fprintf(stderr,"%d % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",t,t*delta, avg_a, avg_e, avg_en, alpha, energy, enstrophy);
       printf("%8d % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",t,t*delta, avg_a, avg_e, avg_en, alpha, energy, enstrophy);
     }
@@ -272,15 +196,6 @@ int eea(
     write_u_bin(u, K1, K2, savefile);
   }
 
-  if(running_avg_window!=0){
-    free(save_print_e);
-    free(save_print_a);
-    free(save_print_en);
-    free(save_print_short_e);
-    free(save_print_short_a);
-    free(save_print_short_en);
-  }
-
   ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft);
   return(0);
 }