From 424a01ec77cabdabc6cd17ae76619dbd0e8e2906 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@ds.mpg.de>
Date: Tue, 29 May 2018 21:41:14 +0200
Subject: [PATCH] tweak spectrum computation
---
bfps/DNS.py | 17 ++++++-----------
bfps/NavierStokes.py | 17 ++++++-----------
bfps/test/test_Parseval.py | 6 ++++--
3 files changed, 16 insertions(+), 24 deletions(-)
diff --git a/bfps/DNS.py b/bfps/DNS.py
index 28291be5..538285e5 100644
--- a/bfps/DNS.py
+++ b/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']))
diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index 62a2263e..d158f5b2 100644
--- a/bfps/NavierStokes.py
+++ b/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']))
diff --git a/bfps/test/test_Parseval.py b/bfps/test/test_Parseval.py
index b2afdbee..80a16d6a 100644
--- a/bfps/test/test_Parseval.py
+++ b/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')
--
GitLab