tweak spectrum computation

bfps/DNS.py

@@ -276,23 +276,18 @@ class DNS(_code):
                                 self.parameters['niter_stat']*
                                 (np.arange(ii0, ii1+1).astype(np.float)))
                 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*(
+                pp_file['R_ij(t)'] = np.sum(phi_ij, axis = 1)
+                energy_tk = (
                     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(t)'] = np.sum(energy_tk, axis = 1)
                 pp_file['energy(k)'] = np.mean(energy_tk, axis = 0)
-                enstrophy_tk = discrete_Fourier_prefactor*(
+                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
-                pp_file['enstrophy(t)'] = (self.statistics['dk'] *
-                                           np.sum(enstrophy_tk, axis = 1))
+                pp_file['enstrophy(t)'] = 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
@@ -373,7 +368,7 @@ 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 /
+        self.statistics['Lint'] = ((np.pi /
                                     (2*self.statistics['Uint']**2)) *
                                    np.nansum(self.statistics['energy(k)'] /
                                                 self.statistics['kshell']))

bfps/NavierStokes.py

@@ -627,23 +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]
-                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*(
+                pp_file['R_ij(t)'] = np.sum(phi_ij, axis = 1)
+                energy_tk = (
                     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(t)'] = np.sum(energy_tk, axis = 1)
                 pp_file['energy(k)'] = np.mean(energy_tk, axis = 0)
-                enstrophy_tk = discrete_Fourier_prefactor*(
+                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
-                pp_file['enstrophy(t)'] = (self.statistics['dk'] *
-                                           np.sum(enstrophy_tk, axis = 1))
+                pp_file['enstrophy(t)'] = 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
@@ -727,7 +722,7 @@ class NavierStokes(_fluid_particle_base):
                                                  self.statistics['etaK' + suffix])
             if self.parameters['dealias_type'] == 1:
                 self.statistics['kMeta' + suffix] *= 0.8
-        self.statistics['Lint'] = ((self.statistics['dk']*np.pi /
+        self.statistics['Lint'] = ((np.pi /
                                     (2*self.statistics['Uint']**2)) *
                                    np.nansum(self.statistics['energy(k)'] /
                                                 self.statistics['kshell']))

bfps/test/test_Parseval.py

@@ -31,13 +31,15 @@ def main():
     energyk = c.statistics['energy(k)']
     nshell = c.get_data_file()['kspace/nshell'].value
     renergy = np.mean(c.statistics['renergy(t)'])
-    print(renergy, np.sum(energyk[:-2]) * (c.parameters['dkx']*c.parameters['dky']*c.parameters['dkz']))
+    print(renergy, np.sum(energyk[:-2]))
     print(c.parameters['dkx'], c.parameters['dky'], c.parameters['dkz'])
+    energyk_alt = (energyk / nshell)*(4*np.pi*c.statistics['kshell']**2)
+    print(renergy, np.sum(energyk_alt[:-2])*c.statistics['dk'] / (c.parameters['dkx']*c.parameters['dky']*c.parameters['dkz']))
 
     f = plt.figure()
     a = f.add_subplot(111)
     a.plot(c.statistics['kshell'], energyk, label = 'unnormalized')
-    a.plot(c.statistics['kshell'], (energyk / nshell)*(4*np.pi*c.statistics['kshell']**2), label = 'normalized')
+    a.plot(c.statistics['kshell'], energyk_alt, label = 'normalized')
     #a.set_yscale('log')
     a.set_xscale('log')
     a.legend(loc = 'best')