From 0dea55d6e91a8aaf10467be932a80c9d2f44c151 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@ds.mpg.de>
Date: Thu, 2 May 2019 10:11:45 +0200
Subject: [PATCH] update spectrum computation

---
 bfps/DNS.py | 18 +++++++++++-------
 1 file changed, 11 insertions(+), 7 deletions(-)

diff --git a/bfps/DNS.py b/bfps/DNS.py
index ce0b2e03..fd22599e 100644
--- a/bfps/DNS.py
+++ b/bfps/DNS.py
@@ -276,6 +276,7 @@ class DNS(_code):
                             if k in pp_file.keys():
                                 del pp_file[k]
                 if computation_needed:
+                    #TODO figure out whether normalization is sane or not
                     pp_file['iter0'] = iter0
                     pp_file['iter1'] = iter1
                     pp_file['ii0'] = ii0
@@ -283,21 +284,24 @@ class DNS(_code):
                     pp_file['t'] = (self.parameters['dt']*
                                     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)'] = np.sum(phi_ij, axis = 1)
+                    # we have an extra division by shell_width because of the Dirac delta restricting integration to the shell
+                    phi_ij = data_file['statistics/spectra/velocity_velocity'][ii0:ii1+1] / self.statistics['dk']
+                    pp_file['R_ij(t)'] = np.sum(phi_ij*self.statistics['dk'], axis = 1)
                     energy_tk = (
                         phi_ij[:, :, 0, 0] +
                         phi_ij[:, :, 1, 1] +
                         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]
+                    pp_file['energy(t)'] = np.sum(energy_tk*self.statistics['dk'], axis = 1)
+                    # normalization factor is (4 pi * shell_width * kshell^2) / (nmodes in shell * dkx*dky*dkz)
+                    norm_factor = (4*np.pi*self.statistics['dk']*self.statistics['kshell']**2) / (self.parameters['dkx']*self.parameters['dky']*self.parameters['dkz']*self.statistics['nshell'])
+                    pp_file['energy(k)'] = np.mean(energy_tk, axis = 0)*norm_factor
+                    phi_vorticity_ij = data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1] / self.statistics['dk']
                     enstrophy_tk = (
                         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['enstrophy(t)'] = np.sum(enstrophy_tk*self.statistics['dk'], axis = 1)
+                    pp_file['enstrophy(k)'] = np.mean(enstrophy_tk, axis = 0)*norm_factor
                     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
-- 
GitLab