Commit dcc3b1c1 authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

fix energy spectrum computation for nontrivial box size

parent eaf6dbd7
Pipeline #29969 canceled with stage
in 4 minutes and 34 seconds
...@@ -208,8 +208,12 @@ class DNS(_code): ...@@ -208,8 +208,12 @@ class DNS(_code):
return os.path.join(self.work_dir, self.simname + '_particles.h5') return os.path.join(self.work_dir, self.simname + '_particles.h5')
def get_particle_file(self): def get_particle_file(self):
return h5py.File(self.get_particle_file_name(), 'r') 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): 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): def get_postprocess_file(self):
return h5py.File(self.get_postprocess_file_name(), 'r') return h5py.File(self.get_postprocess_file_name(), 'r')
def compute_statistics(self, iter0 = 0, iter1 = None): def compute_statistics(self, iter0 = 0, iter1 = None):
...@@ -225,7 +229,7 @@ class DNS(_code): ...@@ -225,7 +229,7 @@ class DNS(_code):
tensors, and the enstrophy spectrum is also used to tensors, and the enstrophy spectrum is also used to
compute the dissipation :math:`\\varepsilon(t)`. compute the dissipation :math:`\\varepsilon(t)`.
These basic quantities are stored in a newly created HDF5 file, These basic quantities are stored in a newly created HDF5 file,
``simname_postprocess.h5``. ``simname_cache.h5``.
""" """
if len(list(self.statistics.keys())) > 0: if len(list(self.statistics.keys())) > 0:
return None return None
...@@ -254,8 +258,15 @@ class DNS(_code): ...@@ -254,8 +258,15 @@ class DNS(_code):
computation_needed = not (ii0 == pp_file['ii0'].value and computation_needed = not (ii0 == pp_file['ii0'].value and
ii1 == pp_file['ii1'].value) ii1 == pp_file['ii1'].value)
if computation_needed: if computation_needed:
for k in pp_file.keys(): for k in ['t', 'vel_max(t)', 'renergy(t)',
del pp_file[k] '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: if computation_needed:
pp_file['iter0'] = iter0 pp_file['iter0'] = iter0
pp_file['iter1'] = iter1 pp_file['iter1'] = iter1
...@@ -264,39 +275,63 @@ class DNS(_code): ...@@ -264,39 +275,63 @@ class DNS(_code):
pp_file['t'] = (self.parameters['dt']* pp_file['t'] = (self.parameters['dt']*
self.parameters['niter_stat']* self.parameters['niter_stat']*
(np.arange(ii0, ii1+1).astype(np.float))) (np.arange(ii0, ii1+1).astype(np.float)))
pp_file['energy(t, k)'] = ( phi_ij = data_file['statistics/spectra/velocity_velocity'][ii0:ii1+1]
data_file['statistics/spectra/velocity_velocity'][ii0:ii1+1, :, 0, 0] + discrete_Fourier_prefactor = 1. / (self.parameters['dkx']*
data_file['statistics/spectra/velocity_velocity'][ii0:ii1+1, :, 1, 1] + self.parameters['dky']*
data_file['statistics/spectra/velocity_velocity'][ii0:ii1+1, :, 2, 2])/2 self.parameters['dkz'])
pp_file['enstrophy(t, k)'] = ( 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, :, 0, 0] +
data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 1, 1] + data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 1, 1] +
data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 2, 2])/2 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 pp_file['renergy(t)'] = data_file['statistics/moments/velocity'][ii0:ii1+1, 2, 3]/2
for k in ['t', for k in ['t',
'energy(t, k)', 'energy(t)',
'enstrophy(t, k)', 'energy(k)',
'enstrophy(t)',
'enstrophy(k)',
'R_ij(t)',
'vel_max(t)', 'vel_max(t)',
'renergy(t)']: 'renergy(t)']:
if k in pp_file.keys(): if k in pp_file.keys():
self.statistics[k] = pp_file[k].value self.statistics[k] = pp_file[k].value
self.compute_time_averages() self.compute_time_averages()
return None 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): def compute_time_averages(self):
"""Compute easy stats. """Compute easy stats.
Further computation of statistics based on the contents of Further computation of statistics based on the contents of
``simname_postprocess.h5``. ``simname_cache.h5``.
Standard quantities are as follows Standard quantities are as follows
(consistent with [Ishihara]_): (consistent with [Ishihara]_):
.. math:: .. math::
U_{\\textrm{int}}(t) = \\sqrt{\\frac{2E(t)}{3}}, \\hskip .5cm 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 L_{\\textrm{int}} = \\frac{\pi}{2U_{int}^2} \\int \\frac{dk}{k} E(k), \\hskip .5cm
T_{\\textrm{int}}(t) = T_{\\textrm{int}} =
\\frac{L_{\\textrm{int}}(t)}{U_{\\textrm{int}}(t)} \\frac{L_{\\textrm{int}}}{U_{\\textrm{int}}}
\\eta_K = \\left(\\frac{\\nu^3}{\\varepsilon}\\right)^{1/4}, \\hskip .5cm \\eta_K = \\left(\\frac{\\nu^3}{\\varepsilon}\\right)^{1/4}, \\hskip .5cm
\\tau_K = \\left(\\frac{\\nu}{\\varepsilon}\\right)^{1/2}, \\hskip .5cm \\tau_K = \\left(\\frac{\\nu}{\\varepsilon}\\right)^{1/2}, \\hskip .5cm
...@@ -313,21 +348,14 @@ class DNS(_code): ...@@ -313,21 +348,14 @@ class DNS(_code):
J. Fluid Mech., J. Fluid Mech.,
**592**, 335-366, 2007 **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['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', for key in ['energy',
'enstrophy', 'enstrophy',
'vel_max', 'mean_trS2',
'Uint', 'Uint']:
'Lint']:
if key + '(t)' in self.statistics.keys(): if key + '(t)' in self.statistics.keys():
self.statistics[key] = np.average(self.statistics[key + '(t)'], axis = 0) 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)']: for suffix in ['', '(t)']:
self.statistics['diss' + suffix] = (self.parameters['nu'] * self.statistics['diss' + suffix] = (self.parameters['nu'] *
self.statistics['enstrophy' + suffix]*2) self.statistics['enstrophy' + suffix]*2)
...@@ -335,9 +363,6 @@ class DNS(_code): ...@@ -335,9 +363,6 @@ class DNS(_code):
self.statistics['diss' + suffix])**.25 self.statistics['diss' + suffix])**.25
self.statistics['tauK' + suffix] = (self.parameters['nu'] / self.statistics['tauK' + suffix] = (self.parameters['nu'] /
self.statistics['diss' + suffix])**.5 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['lambda' + suffix] = (15 * self.parameters['nu'] *
self.statistics['Uint' + suffix]**2 / self.statistics['Uint' + suffix]**2 /
self.statistics['diss' + suffix])**.5 self.statistics['diss' + suffix])**.5
...@@ -348,6 +373,13 @@ class DNS(_code): ...@@ -348,6 +373,13 @@ class DNS(_code):
self.statistics['etaK' + suffix]) self.statistics['etaK' + suffix])
if self.parameters['dealias_type'] == 1: if self.parameters['dealias_type'] == 1:
self.statistics['kMeta' + suffix] *= 0.8 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['Tint'] = self.statistics['Lint'] / self.statistics['Uint']
self.statistics['Taylor_microscale'] = self.statistics['lambda'] self.statistics['Taylor_microscale'] = self.statistics['lambda']
return None return None
...@@ -836,21 +868,21 @@ class DNS(_code): ...@@ -836,21 +868,21 @@ class DNS(_code):
""" """
np.random.seed(rseed) np.random.seed(rseed)
Kdata00 = scalar_generator( Kdata00 = scalar_generator(
self.parameters['nz']//2, self.parameters['nz'],
self.parameters['ny']//2, self.parameters['ny'],
self.parameters['nx']//2, self.parameters['nx'],
p = spectra_slope, p = spectra_slope,
amplitude = amplitude).astype(self.ctype) amplitude = amplitude).astype(self.ctype)
Kdata01 = scalar_generator( Kdata01 = scalar_generator(
self.parameters['nz']//2, self.parameters['nz'],
self.parameters['ny']//2, self.parameters['ny'],
self.parameters['nx']//2, self.parameters['nx'],
p = spectra_slope, p = spectra_slope,
amplitude = amplitude).astype(self.ctype) amplitude = amplitude).astype(self.ctype)
Kdata02 = scalar_generator( Kdata02 = scalar_generator(
self.parameters['nz']//2, self.parameters['nz'],
self.parameters['ny']//2, self.parameters['ny'],
self.parameters['nx']//2, self.parameters['nx'],
p = spectra_slope, p = spectra_slope,
amplitude = amplitude).astype(self.ctype) amplitude = amplitude).astype(self.ctype)
Kdata0 = np.zeros( Kdata0 = np.zeros(
......
...@@ -627,15 +627,18 @@ class NavierStokes(_fluid_particle_base): ...@@ -627,15 +627,18 @@ class NavierStokes(_fluid_particle_base):
self.parameters['niter_stat']* self.parameters['niter_stat']*
(np.arange(ii0, ii1+1).astype(np.float))) (np.arange(ii0, ii1+1).astype(np.float)))
phi_ij = data_file['statistics/spectra/velocity_velocity'][ii0:ii1+1] 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) discrete_Fourier_prefactor = 1. / (self.parameters['dkx']*
energy_tk = ( 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[:, :, 0, 0] +
phi_ij[:, :, 1, 1] + phi_ij[:, :, 1, 1] +
phi_ij[:, :, 2, 2])/2 phi_ij[:, :, 2, 2])/2
pp_file['energy(t)'] = (self.statistics['dk'] * pp_file['energy(t)'] = (self.statistics['dk'] *
np.sum(energy_tk, axis = 1)) 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)
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, :, 0, 0] +
data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 1, 1] + data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 1, 1] +
data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 2, 2])/2 data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 2, 2])/2
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment