From dcc3b1c19107ac37765856b69c261b4167b40f7f Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@ds.mpg.de>
Date: Tue, 29 May 2018 12:32:51 +0200
Subject: [PATCH] fix energy spectrum computation for nontrivial box size
---
bfps/DNS.py | 108 ++++++++++++++++++++++++++++---------------
bfps/NavierStokes.py | 9 ++--
2 files changed, 76 insertions(+), 41 deletions(-)
diff --git a/bfps/DNS.py b/bfps/DNS.py
index 7c5fdfa7..c407a955 100644
--- a/bfps/DNS.py
+++ b/bfps/DNS.py
@@ -208,8 +208,12 @@ class DNS(_code):
return os.path.join(self.work_dir, self.simname + '_particles.h5')
def get_particle_file(self):
return h5py.File(self.get_particle_file_name(), 'r')
+ def get_cache_file_name(self):
+ return os.path.join(self.work_dir, self.simname + '_cache.h5')
+ def get_cache_file(self):
+ return h5py.File(self.get_cache_file_name(), 'r')
def get_postprocess_file_name(self):
- return os.path.join(self.work_dir, self.simname + '_postprocess.h5')
+ return self.get_cache_file_name()
def get_postprocess_file(self):
return h5py.File(self.get_postprocess_file_name(), 'r')
def compute_statistics(self, iter0 = 0, iter1 = None):
@@ -225,7 +229,7 @@ class DNS(_code):
tensors, and the enstrophy spectrum is also used to
compute the dissipation :math:`\\varepsilon(t)`.
These basic quantities are stored in a newly created HDF5 file,
- ``simname_postprocess.h5``.
+ ``simname_cache.h5``.
"""
if len(list(self.statistics.keys())) > 0:
return None
@@ -254,8 +258,15 @@ class DNS(_code):
computation_needed = not (ii0 == pp_file['ii0'].value and
ii1 == pp_file['ii1'].value)
if computation_needed:
- for k in pp_file.keys():
- del pp_file[k]
+ for k in ['t', 'vel_max(t)', 'renergy(t)',
+ 'energy(t)', 'enstrophy(t)',
+ 'energy(k)', 'enstrophy(k)',
+ 'energy(t, k)',
+ 'enstrophy(t, k)',
+ 'R_ij(t)',
+ 'ii0', 'ii1', 'iter0', 'iter1']:
+ if k in pp_file.keys():
+ del pp_file[k]
if computation_needed:
pp_file['iter0'] = iter0
pp_file['iter1'] = iter1
@@ -264,39 +275,63 @@ class DNS(_code):
pp_file['t'] = (self.parameters['dt']*
self.parameters['niter_stat']*
(np.arange(ii0, ii1+1).astype(np.float)))
- pp_file['energy(t, k)'] = (
- data_file['statistics/spectra/velocity_velocity'][ii0:ii1+1, :, 0, 0] +
- data_file['statistics/spectra/velocity_velocity'][ii0:ii1+1, :, 1, 1] +
- data_file['statistics/spectra/velocity_velocity'][ii0:ii1+1, :, 2, 2])/2
- pp_file['enstrophy(t, k)'] = (
+ 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*(
+ 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(k)'] = np.mean(energy_tk, axis = 0)
+ enstrophy_tk = discrete_Fourier_prefactor*(
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['vel_max(t)'] = data_file['statistics/moments/velocity'] [ii0:ii1+1, 9, 3]
+ pp_file['enstrophy(t)'] = (self.statistics['dk'] *
+ 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
for k in ['t',
- 'energy(t, k)',
- 'enstrophy(t, k)',
+ 'energy(t)',
+ 'energy(k)',
+ 'enstrophy(t)',
+ 'enstrophy(k)',
+ 'R_ij(t)',
'vel_max(t)',
'renergy(t)']:
if k in pp_file.keys():
self.statistics[k] = pp_file[k].value
self.compute_time_averages()
return None
+ def compute_Reynolds_stress_invariants(
+ self):
+ Rij = self.statistics['R_ij(t)']
+ Rij /= (2*self.statistics['energy(t)'][:, None, None])
+ Rij[:, 0, 0] -= 1./3
+ Rij[:, 1, 1] -= 1./3
+ Rij[:, 2, 2] -= 1./3
+ self.statistics['I2(t)'] = np.sqrt(np.einsum('...ij,...ij', Rij, Rij, optimize = True) / 6)
+ self.statistics['I3(t)'] = np.cbrt(np.einsum('...ij,...jk,...ki', Rij, Rij, Rij, optimize = True) / 6)
+ return None
def compute_time_averages(self):
"""Compute easy stats.
Further computation of statistics based on the contents of
- ``simname_postprocess.h5``.
+ ``simname_cache.h5``.
Standard quantities are as follows
(consistent with [Ishihara]_):
.. math::
U_{\\textrm{int}}(t) = \\sqrt{\\frac{2E(t)}{3}}, \\hskip .5cm
- L_{\\textrm{int}}(t) = \\frac{\pi}{2U_{int}^2(t)} \\int \\frac{dk}{k} E(t, k), \\hskip .5cm
- T_{\\textrm{int}}(t) =
- \\frac{L_{\\textrm{int}}(t)}{U_{\\textrm{int}}(t)}
+ L_{\\textrm{int}} = \\frac{\pi}{2U_{int}^2} \\int \\frac{dk}{k} E(k), \\hskip .5cm
+ T_{\\textrm{int}} =
+ \\frac{L_{\\textrm{int}}}{U_{\\textrm{int}}}
\\eta_K = \\left(\\frac{\\nu^3}{\\varepsilon}\\right)^{1/4}, \\hskip .5cm
\\tau_K = \\left(\\frac{\\nu}{\\varepsilon}\\right)^{1/2}, \\hskip .5cm
@@ -313,21 +348,14 @@ class DNS(_code):
J. Fluid Mech.,
**592**, 335-366, 2007
"""
- for key in ['energy', 'enstrophy']:
- self.statistics[key + '(t)'] = (self.statistics['dk'] *
- np.sum(self.statistics[key + '(t, k)'], axis = 1))
self.statistics['Uint(t)'] = np.sqrt(2*self.statistics['energy(t)'] / 3)
- self.statistics['Lint(t)'] = ((self.statistics['dk']*np.pi /
- (2*self.statistics['Uint(t)']**2)) *
- np.nansum(self.statistics['energy(t, k)'] /
- self.statistics['kshell'][None, :], axis = 1))
for key in ['energy',
'enstrophy',
- 'vel_max',
- 'Uint',
- 'Lint']:
+ 'mean_trS2',
+ 'Uint']:
if key + '(t)' in self.statistics.keys():
self.statistics[key] = np.average(self.statistics[key + '(t)'], axis = 0)
+ self.statistics['vel_max'] = np.max(self.statistics['vel_max(t)'])
for suffix in ['', '(t)']:
self.statistics['diss' + suffix] = (self.parameters['nu'] *
self.statistics['enstrophy' + suffix]*2)
@@ -335,9 +363,6 @@ class DNS(_code):
self.statistics['diss' + suffix])**.25
self.statistics['tauK' + suffix] = (self.parameters['nu'] /
self.statistics['diss' + suffix])**.5
- self.statistics['Re' + suffix] = (self.statistics['Uint' + suffix] *
- self.statistics['Lint' + suffix] /
- self.parameters['nu'])
self.statistics['lambda' + suffix] = (15 * self.parameters['nu'] *
self.statistics['Uint' + suffix]**2 /
self.statistics['diss' + suffix])**.5
@@ -348,6 +373,13 @@ 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 /
+ (2*self.statistics['Uint']**2)) *
+ np.nansum(self.statistics['energy(k)'] /
+ self.statistics['kshell']))
+ self.statistics['Re'] = (self.statistics['Uint'] *
+ self.statistics['Lint'] /
+ self.parameters['nu'])
self.statistics['Tint'] = self.statistics['Lint'] / self.statistics['Uint']
self.statistics['Taylor_microscale'] = self.statistics['lambda']
return None
@@ -836,21 +868,21 @@ class DNS(_code):
"""
np.random.seed(rseed)
Kdata00 = scalar_generator(
- self.parameters['nz']//2,
- self.parameters['ny']//2,
- self.parameters['nx']//2,
+ self.parameters['nz'],
+ self.parameters['ny'],
+ self.parameters['nx'],
p = spectra_slope,
amplitude = amplitude).astype(self.ctype)
Kdata01 = scalar_generator(
- self.parameters['nz']//2,
- self.parameters['ny']//2,
- self.parameters['nx']//2,
+ self.parameters['nz'],
+ self.parameters['ny'],
+ self.parameters['nx'],
p = spectra_slope,
amplitude = amplitude).astype(self.ctype)
Kdata02 = scalar_generator(
- self.parameters['nz']//2,
- self.parameters['ny']//2,
- self.parameters['nx']//2,
+ self.parameters['nz'],
+ self.parameters['ny'],
+ self.parameters['nx'],
p = spectra_slope,
amplitude = amplitude).astype(self.ctype)
Kdata0 = np.zeros(
diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index 98169ab0..62a2263e 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -627,15 +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]
- pp_file['R_ij(t)'] = self.statistics['dk']*np.sum(phi_ij, axis = 1)
- energy_tk = (
+ 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*(
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(k)'] = np.mean(energy_tk, axis = 0)
- enstrophy_tk = (
+ enstrophy_tk = discrete_Fourier_prefactor*(
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
--
GitLab