running average in statistics file

src/navier-stokes.c

@@ -1,5 +1,6 @@
-#include "navier-stokes.h"
 #include "io.h"
+#include "navier-stokes.h"
+#include "statistics.h"
 #include <math.h>
 #include <stdlib.h>
 #include <string.h>
@@ -132,25 +133,9 @@ int eea(
     alpha=compute_alpha(u, K1, K2, g, L);
     enstrophy=compute_enstrophy(u, K1, K2, L);
 
-    // running average
+    avg_e=average_step(energy, avg_e, t, starting_time, print_freq, first_box);
-    // reset averages
+    avg_a=average_step(alpha, avg_a, t, starting_time, print_freq, first_box);
-    if(t % print_freq == 1){
+    avg_en=average_step(enstrophy, avg_en, t, starting_time, print_freq, first_box);
-      avg_e=0;
-      avg_a=0;
-      avg_en=0;
-    }
-
-    // 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>starting_time && t%print_freq==0){
       fprintf(stderr,"%lu % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",t,t*delta, avg_a, avg_e, avg_en, alpha, energy, enstrophy);

src/statistics.c

@@ -0,0 +1,28 @@
+#include "statistics.h"
+
+// run this at each step to compute the running average
+double average_step(
+  // the value at the step
+  double val,
+  // the average computed so far
+  double avg,
+  uint64_t t,
+  uint64_t starting_time,
+  uint64_t print_freq,
+  uint64_t first_box
+){
+
+  // running average
+  // reset averages
+  if(t % print_freq == 1){
+    return(0);
+  }
+
+  // compute average
+  // different computationin first box if starting_time is not a multiple of print_freq
+  if(t < starting_time + first_box){
+    return(avg+val/first_box);
+  } else {
+    return(avg+val/print_freq);
+  }
+}

src/statistics.h

@@ -0,0 +1,9 @@
+#ifndef STATISTICS_H
+#define STATISTICS_H
+
+#include <stdint.h>
+
+// run this at each step to compute the running average
+double average_step( double val, double avg, uint64_t t, uint64_t starting_time, uint64_t print_freq, uint64_t first_box);
+
+#endif