diff --git a/bfps/DNS.py b/bfps/DNS.py
index 7c5fdfa7dfd97f9ded6ca387d32af9b7eec4e579..c407a955c3d7e9ffd50bd69850b3b64296c65aa1 100644
--- a/bfps/DNS.py
+++ b/bfps/DNS.py
@@ -208,8 +208,12 @@ class DNS(_code):
         return os.path.join(self.work_dir, self.simname + '_particles.h5')
     def get_particle_file(self):
         return h5py.File(self.get_particle_file_name(), 'r')
+    def get_cache_file_name(self):
+        return os.path.join(self.work_dir, self.simname + '_cache.h5')
+    def get_cache_file(self):
+        return h5py.File(self.get_cache_file_name(), 'r')
     def get_postprocess_file_name(self):
-        return os.path.join(self.work_dir, self.simname + '_postprocess.h5')
+        return self.get_cache_file_name()
     def get_postprocess_file(self):
         return h5py.File(self.get_postprocess_file_name(), 'r')
     def compute_statistics(self, iter0 = 0, iter1 = None):
@@ -225,7 +229,7 @@ class DNS(_code):
         tensors, and the enstrophy spectrum is also used to
         compute the dissipation :math:`\\varepsilon(t)`.
         These basic quantities are stored in a newly created HDF5 file,
-        ``simname_postprocess.h5``.
+        ``simname_cache.h5``.
         """
         if len(list(self.statistics.keys())) > 0:
             return None
@@ -254,8 +258,15 @@ class DNS(_code):
                 computation_needed =  not (ii0 == pp_file['ii0'].value and
                                            ii1 == pp_file['ii1'].value)
                 if computation_needed:
-                    for k in pp_file.keys():
-                        del pp_file[k]
+                    for k in ['t', 'vel_max(t)', 'renergy(t)',
+                              'energy(t)', 'enstrophy(t)',
+                              'energy(k)', 'enstrophy(k)',
+                              'energy(t, k)',
+                              'enstrophy(t, k)',
+                              'R_ij(t)',
+                              'ii0', 'ii1', 'iter0', 'iter1']:
+                        if k in pp_file.keys():
+                            del pp_file[k]
             if computation_needed:
                 pp_file['iter0'] = iter0
                 pp_file['iter1'] = iter1
@@ -264,39 +275,63 @@ class DNS(_code):
                 pp_file['t'] = (self.parameters['dt']*
                                 self.parameters['niter_stat']*
                                 (np.arange(ii0, ii1+1).astype(np.float)))
-                pp_file['energy(t, k)'] = (
-                    data_file['statistics/spectra/velocity_velocity'][ii0:ii1+1, :, 0, 0] +
-                    data_file['statistics/spectra/velocity_velocity'][ii0:ii1+1, :, 1, 1] +
-                    data_file['statistics/spectra/velocity_velocity'][ii0:ii1+1, :, 2, 2])/2
-                pp_file['enstrophy(t, k)'] = (
+                phi_ij = data_file['statistics/spectra/velocity_velocity'][ii0:ii1+1]
+                discrete_Fourier_prefactor = 1. / (self.parameters['dkx']*
+                                                   self.parameters['dky']*
+                                                   self.parameters['dkz'])
+                pp_file['R_ij(t)'] = self.statistics['dk']*np.sum(phi_ij, axis = 1)*discrete_Fourier_prefactor
+                energy_tk = discrete_Fourier_prefactor*(
+                    phi_ij[:, :, 0, 0] +
+                    phi_ij[:, :, 1, 1] +
+                    phi_ij[:, :, 2, 2])/2
+                pp_file['energy(t)'] = (self.statistics['dk'] *
+                                        np.sum(energy_tk, axis = 1))
+                pp_file['energy(k)'] = np.mean(energy_tk, axis = 0)
+                enstrophy_tk = discrete_Fourier_prefactor*(
                     data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 0, 0] +
                     data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 1, 1] +
                     data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 2, 2])/2
-                pp_file['vel_max(t)'] = data_file['statistics/moments/velocity']  [ii0:ii1+1, 9, 3]
+                pp_file['enstrophy(t)'] = (self.statistics['dk'] *
+                                           np.sum(enstrophy_tk, axis = 1))
+                pp_file['enstrophy(k)'] = np.mean(enstrophy_tk, axis = 0)
+                pp_file['vel_max(t)'] = data_file['statistics/moments/velocity'][ii0:ii1+1, 9, 3]
                 pp_file['renergy(t)'] = data_file['statistics/moments/velocity'][ii0:ii1+1, 2, 3]/2
             for k in ['t',
-                      'energy(t, k)',
-                      'enstrophy(t, k)',
+                      'energy(t)',
+                      'energy(k)',
+                      'enstrophy(t)',
+                      'enstrophy(k)',
+                      'R_ij(t)',
                       'vel_max(t)',
                       'renergy(t)']:
                 if k in pp_file.keys():
                     self.statistics[k] = pp_file[k].value
             self.compute_time_averages()
         return None
+    def compute_Reynolds_stress_invariants(
+            self):
+        Rij = self.statistics['R_ij(t)']
+        Rij /= (2*self.statistics['energy(t)'][:, None, None])
+        Rij[:, 0, 0] -= 1./3
+        Rij[:, 1, 1] -= 1./3
+        Rij[:, 2, 2] -= 1./3
+        self.statistics['I2(t)'] = np.sqrt(np.einsum('...ij,...ij', Rij, Rij, optimize = True) / 6)
+        self.statistics['I3(t)'] = np.cbrt(np.einsum('...ij,...jk,...ki', Rij, Rij, Rij, optimize = True) / 6)
+        return None
     def compute_time_averages(self):
         """Compute easy stats.
 
         Further computation of statistics based on the contents of
-        ``simname_postprocess.h5``.
+        ``simname_cache.h5``.
         Standard quantities are as follows
         (consistent with [Ishihara]_):
 
         .. math::
 
             U_{\\textrm{int}}(t) = \\sqrt{\\frac{2E(t)}{3}}, \\hskip .5cm
-            L_{\\textrm{int}}(t) = \\frac{\pi}{2U_{int}^2(t)} \\int \\frac{dk}{k} E(t, k), \\hskip .5cm
-            T_{\\textrm{int}}(t) =
-            \\frac{L_{\\textrm{int}}(t)}{U_{\\textrm{int}}(t)}
+            L_{\\textrm{int}} = \\frac{\pi}{2U_{int}^2} \\int \\frac{dk}{k} E(k), \\hskip .5cm
+            T_{\\textrm{int}} =
+            \\frac{L_{\\textrm{int}}}{U_{\\textrm{int}}}
 
             \\eta_K = \\left(\\frac{\\nu^3}{\\varepsilon}\\right)^{1/4}, \\hskip .5cm
             \\tau_K = \\left(\\frac{\\nu}{\\varepsilon}\\right)^{1/2}, \\hskip .5cm
@@ -313,21 +348,14 @@ class DNS(_code):
                       J. Fluid Mech.,
                       **592**, 335-366, 2007
         """
-        for key in ['energy', 'enstrophy']:
-            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',
-                    'Uint',
-                    'Lint']:
+                    'mean_trS2',
+                    'Uint']:
             if key + '(t)' in self.statistics.keys():
                 self.statistics[key] = np.average(self.statistics[key + '(t)'], axis = 0)
+        self.statistics['vel_max'] = np.max(self.statistics['vel_max(t)'])
         for suffix in ['', '(t)']:
             self.statistics['diss'    + suffix] = (self.parameters['nu'] *
                                                    self.statistics['enstrophy' + suffix]*2)
@@ -335,9 +363,6 @@ class DNS(_code):
                                                    self.statistics['diss' + suffix])**.25
             self.statistics['tauK'    + suffix] =  (self.parameters['nu'] /
                                                     self.statistics['diss' + suffix])**.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
@@ -348,6 +373,13 @@ class DNS(_code):
                                                  self.statistics['etaK' + suffix])
             if self.parameters['dealias_type'] == 1:
                 self.statistics['kMeta' + suffix] *= 0.8
+        self.statistics['Lint'] = ((self.statistics['dk']*np.pi /
+                                    (2*self.statistics['Uint']**2)) *
+                                   np.nansum(self.statistics['energy(k)'] /
+                                                self.statistics['kshell']))
+        self.statistics['Re'] = (self.statistics['Uint'] *
+                                 self.statistics['Lint'] /
+                                 self.parameters['nu'])
         self.statistics['Tint'] = self.statistics['Lint'] / self.statistics['Uint']
         self.statistics['Taylor_microscale'] = self.statistics['lambda']
         return None
@@ -836,21 +868,21 @@ class DNS(_code):
         """
         np.random.seed(rseed)
         Kdata00 = scalar_generator(
-                self.parameters['nz']//2,
-                self.parameters['ny']//2,
-                self.parameters['nx']//2,
+                self.parameters['nz'],
+                self.parameters['ny'],
+                self.parameters['nx'],
                 p = spectra_slope,
                 amplitude = amplitude).astype(self.ctype)
         Kdata01 = scalar_generator(
-                self.parameters['nz']//2,
-                self.parameters['ny']//2,
-                self.parameters['nx']//2,
+                self.parameters['nz'],
+                self.parameters['ny'],
+                self.parameters['nx'],
                 p = spectra_slope,
                 amplitude = amplitude).astype(self.ctype)
         Kdata02 = scalar_generator(
-                self.parameters['nz']//2,
-                self.parameters['ny']//2,
-                self.parameters['nx']//2,
+                self.parameters['nz'],
+                self.parameters['ny'],
+                self.parameters['nx'],
                 p = spectra_slope,
                 amplitude = amplitude).astype(self.ctype)
         Kdata0 = np.zeros(
diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index 98169ab00b12981f4fbbc59ca417306a7f3226b6..62a2263e9ed65e5dbc7ab32738a0430ad6436a5d 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -627,15 +627,18 @@ class NavierStokes(_fluid_particle_base):
                                 self.parameters['niter_stat']*
                                 (np.arange(ii0, ii1+1).astype(np.float)))
                 phi_ij = data_file['statistics/spectra/velocity_velocity'][ii0:ii1+1]
-                pp_file['R_ij(t)'] = self.statistics['dk']*np.sum(phi_ij, axis = 1)
-                energy_tk = (
+                discrete_Fourier_prefactor = 1. / (self.parameters['dkx']*
+                                                   self.parameters['dky']*
+                                                   self.parameters['dkz'])
+                pp_file['R_ij(t)'] = self.statistics['dk']*np.sum(phi_ij, axis = 1)*discrete_Fourier_prefactor
+                energy_tk = discrete_Fourier_prefactor*(
                     phi_ij[:, :, 0, 0] +
                     phi_ij[:, :, 1, 1] +
                     phi_ij[:, :, 2, 2])/2
                 pp_file['energy(t)'] = (self.statistics['dk'] *
                                         np.sum(energy_tk, axis = 1))
                 pp_file['energy(k)'] = np.mean(energy_tk, axis = 0)
-                enstrophy_tk = (
+                enstrophy_tk = discrete_Fourier_prefactor*(
                     data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 0, 0] +
                     data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 1, 1] +
                     data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 2, 2])/2