From d265e8779b27adabc85941fb9ec47c553cd90f0c Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@ds.mpg.de>
Date: Wed, 30 May 2018 13:08:31 +0200
Subject: [PATCH] use "normalized" spectrum

I'm not sure how much sense this makes, but I am now explicitly dividing
the computed spectrum by dk, such that the integral of E(k)dk is the
total energy.
---
 bfps/DNS.py          | 5 +++--
 bfps/NavierStokes.py | 5 +++--
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/bfps/DNS.py b/bfps/DNS.py
index 538285e5..90a42c79 100644
--- a/bfps/DNS.py
+++ b/bfps/DNS.py
@@ -247,6 +247,7 @@ class DNS(_code):
             ii0 = iter0 // self.parameters['niter_stat']
             ii1 = iter1 // self.parameters['niter_stat']
             self.statistics['kshell'] = data_file['kspace/kshell'].value
+            self.statistics['nshell'] = data_file['kspace/nshell'].value
             for kk in [-1, -2]:
                 if (self.statistics['kshell'][kk] == 0):
                     self.statistics['kshell'][kk] = np.nan
@@ -282,13 +283,13 @@ class DNS(_code):
                     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)
+                pp_file['energy(k)'] = np.mean(energy_tk, axis = 0)*(4*np.pi*self.statistics['kshell']**2) / (self.statistics['dk']*self.statistics['nshell'])
                 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)'] = np.sum(enstrophy_tk, axis = 1)
-                pp_file['enstrophy(k)'] = np.mean(enstrophy_tk, axis = 0)
+                pp_file['enstrophy(k)'] = np.mean(enstrophy_tk, axis = 0)*(4*np.pi*self.statistics['kshell']**2) / 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
             for k in ['t',
diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index d158f5b2..9c12ca20 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -598,6 +598,7 @@ class NavierStokes(_fluid_particle_base):
             ii0 = iter0 // self.parameters['niter_stat']
             ii1 = iter1 // self.parameters['niter_stat']
             self.statistics['kshell'] = data_file['kspace/kshell'].value
+            self.statistics['nshell'] = data_file['kspace/nshell'].value
             for kk in [-1, -2]:
                 if (self.statistics['kshell'][kk] == 0):
                     self.statistics['kshell'][kk] = np.nan
@@ -633,13 +634,13 @@ class NavierStokes(_fluid_particle_base):
                     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)
+                pp_file['energy(k)'] = np.mean(energy_tk, axis = 0)*(4*np.pi*self.statistics['kshell']**2) / (self.statistics['dk']*self.statistics['nshell'])
                 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)'] = np.sum(enstrophy_tk, axis = 1)
-                pp_file['enstrophy(k)'] = np.mean(enstrophy_tk, axis = 0)
+                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
                 if 'trS2_Q_R' in data_file['statistics/moments'].keys():
-- 
GitLab