add Parseval test for vorticity

bfps/DNS.py

@@ -226,7 +226,8 @@ class DNS(_code):
                           'enstrophy(k)',
                           'R_ij(t)',
                           'vel_max(t)',
-                          'renergy(t)']:
+                          'renergy(t)',
+                          'renstrophy(t)']:
                     if k in pp_file.keys():
                         self.statistics[k] = pp_file[k][...]
                 self.statistics['kM'] = pp_file['kspace/kM'][...]
@@ -263,7 +264,9 @@ class DNS(_code):
                     computation_needed =  not (ii0 == pp_file['ii0'][...] and
                                                ii1 == pp_file['ii1'][...])
                     if computation_needed:
-                        for k in ['t', 'vel_max(t)', 'renergy(t)',
+                        for k in ['t', 'vel_max(t)',
+                                  'renergy(t)',
+                                  'renstrophy(t)',
                                   'energy(t)', 'enstrophy(t)',
                                   'energy(k)', 'enstrophy(k)',
                                   'energy(t, k)',
@@ -288,14 +291,16 @@ class DNS(_code):
                         phi_ij[:, :, 2, 2])/2
                     pp_file['energy(t)'] = np.sum(energy_tk, axis = 1)
                     pp_file['energy(k)'] = np.mean(energy_tk, axis = 0)*(4*np.pi*self.statistics['kshell']**2) / (self.statistics['dk']*self.statistics['nshell'])
+                    phi_vorticity_ij = data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1]
                     enstrophy_tk = (
-                        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
+                        phi_vorticity_ij[:, :, 0, 0] +
+                        phi_vorticity_ij[:, :, 1, 1] +
+                        phi_vorticity_ij[:, :, 2, 2])/2
                     pp_file['enstrophy(t)'] = np.sum(enstrophy_tk, axis = 1)
                     pp_file['enstrophy(k)'] = np.mean(enstrophy_tk, axis = 0)*(4*np.pi*self.statistics['kshell']**2) / (self.statistics['dk']*self.statistics['nshell'])
                     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
+                    pp_file['renstrophy(t)'] = data_file['statistics/moments/vorticity'][ii0:ii1+1, 2, 3]/2
         for k in ['t',
                   'energy(t)',
                   'energy(k)',
@@ -303,13 +308,17 @@ class DNS(_code):
                   'enstrophy(k)',
                   'R_ij(t)',
                   'vel_max(t)',
-                  'renergy(t)']:
+                  'renergy(t)',
+                  'renstrophy(t)']:
             if k in pp_file.keys():
                 self.statistics[k] = pp_file[k][...]
         # sanity check --- Parseval theorem check
         assert(np.max(np.abs(
                 self.statistics['renergy(t)'] -
                 self.statistics['energy(t)']) / self.statistics['energy(t)']) < 1e-5)
+        assert(np.max(np.abs(
+                self.statistics['renstrophy(t)'] -
+                self.statistics['enstrophy(t)']) / self.statistics['enstrophy(t)']) < 1e-5)
         self.compute_time_averages()
         return None
     def compute_Reynolds_stress_invariants(