standardize statistics

I'm not using the ideal estimate of the integral length since I'm doing
the integral with the 1D spectrum, but I think this should be close
enough.

bfps/NavierStokes.py

@@ -530,6 +530,7 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
             ii1 = iter1 // self.parameters['niter_stat']
             self.statistics['kshell'] = data_file['kspace/kshell'].value
             self.statistics['kM'] = data_file['kspace/kM'].value
+            self.statistics['dk'] = data_file['kspace/dk'].value
             if self.particle_species > 0:
                 self.trajectories = [
                         data_file['particles/' + key + '/state'][
@@ -575,8 +576,18 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
         return None
     def compute_time_averages(self):
         for key in ['energy', 'enstrophy']:
-            self.statistics[key + '(t)'] = np.sum(self.statistics[key + '(t, k)'], axis = 1)
-        for key in ['energy', 'enstrophy', 'vel_max', 'mean_trS2']:
+            self.statistics[key + '(t)'] = (self.statistics['dk'] *
+                                            np.sum(self.statistics[key + '(t, k)'], axis = 1))
+        self.statistics['Uint(t)'] = np.sqrt(2*self.statistics['energy(t)'] / 3)
+        self.statistics['Lint(t)'] = ((self.statistics['dk']*np.pi / (2*self.statistics['Uint(t)']**2)) *
+                                      np.nansum(self.statistics['energy(t, k)'] /
+                                                self.statistics['kshell'][None, :], axis = 1))
+        for key in ['energy',
+                    'enstrophy',
+                    'vel_max',
+                    'mean_trS2',
+                    'Uint',
+                    'Lint']:
             if key + '(t)' in self.statistics.keys():
                 self.statistics[key] = np.average(self.statistics[key + '(t)'], axis = 0)
         for suffix in ['', '(t)']:
@@ -584,14 +595,19 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
                                                    self.statistics['enstrophy' + suffix]*2)
             self.statistics['etaK'    + suffix] = (self.parameters['nu']**3 /
                                                    self.statistics['diss' + suffix])**.25
-            self.statistics['Rlambda' + suffix] = (2*np.sqrt(5./3) *
-                                                   (self.statistics['energy' + suffix] /
-                                                   (self.parameters['nu']*self.statistics['diss' + suffix])**.5))
             self.statistics['tauK'    + suffix] =  (self.parameters['nu'] /
                                                     self.statistics['diss' + suffix])**.5
-        self.statistics['Tint'] = 2*self.statistics['energy'] / self.statistics['diss']
-        self.statistics['Lint'] = (2*self.statistics['energy'])**1.5 / self.statistics['diss']
-        self.statistics['Taylor_microscale'] = (10 * self.parameters['nu'] * self.statistics['energy'] / self.statistics['diss'])**.5
+            self.statistics['Re' + suffix] = (self.statistics['Uint' + suffix] *
+                                              self.statistics['Lint' + suffix] /
+                                              self.parameters['nu'])
+            self.statistics['lambda' + suffix] = (15 * self.parameters['nu'] *
+                                                  self.statistics['Uint' + suffix]**2 /
+                                                  self.statistics['diss' + suffix])**.5
+            self.statistics['Rlambda' + suffix] = (self.statistics['Uint' + suffix] *
+                                                   self.statistics['lambda' + suffix] /
+                                                   self.parameters['nu'])
+        self.statistics['Tint'] = self.statistics['Lint'] / self.statistics['Uint']
+        self.statistics['Taylor_microscale'] = self.statistics['lambda']
         return None
     def set_plt_style(
             self,

tests/test_base.py

@@ -96,7 +96,10 @@ def launch(
     c.parameters['famplitude'] = 0.2
     c.fill_up_fluid_code()
     if c.parameters['nparticles'] > 0:
-        c.add_particle_fields(name = 'regular', neighbours = opt.neighbours, smoothness = opt.smoothness)
+        c.add_particle_fields(
+                name = 'regular',
+                neighbours = opt.neighbours,
+                smoothness = opt.smoothness)
         c.add_particle_fields(kcut = 'fs->kM/2', name = 'filtered', neighbours = opt.neighbours)
         c.add_particles(
                 kcut = 'fs->kM/2',
@@ -108,9 +111,8 @@ def launch(
         #            neighbours = opt.neighbours,
         #            smoothness = opt.smoothness,
         #            fields_name = 'regular')
-        for info in [(2, 'Heun'),
-                     (2, 'AdamsBashforth'),
-                     (4, 'cRK4'),
+        for info in [(2, 'AdamsBashforth'),
+                     (3, 'AdamsBashforth'),
                      (4, 'AdamsBashforth'),
                      (6, 'AdamsBashforth')]:
             c.add_particles(

tests/test_plain.py

@@ -38,7 +38,11 @@ def plain(opt):
     opt.work_dir = wd + '/N{0:0>3x}_1'.format(opt.n)
     c0 = launch(opt, dt = 0.2/opt.n)
     c0.compute_statistics()
-    df = c0.get_data_file()
+    print ('Re = {0:.0f}'.format(c0.statistics['Re']))
+    print ('Rlambda = {0:.0f}'.format(c0.statistics['Rlambda']))
+    print ('Lint = {0:.4e}, etaK = {1:.4e}'.format(c0.statistics['Lint'], c0.statistics['etaK']))
+    print ('Tint = {0:.4e}, tauK = {1:.4e}'.format(c0.statistics['Tint'], c0.statistics['tauK']))
+    print ('kMetaK = {0:.4e}'.format(c0.statistics['kM']*c0.statistics['etaK']))
     for s in range(c0.particle_species):
         acceleration_test(c0, species = s, m = 1)
     if not opt.multiplejob: