Commit 0dea55d6 authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

update spectrum computation

parent 35e5f37b
Pipeline #47546 failed with stage
in 3 seconds
......@@ -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
......
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