diff --git a/bfps/DNS.py b/bfps/DNS.py
index 90a42c795e2b3ccd6661f3c7ba9c77dad0901664..495278bfde4fd6cce49d4530b26403d0d36f9d9d 100644
--- a/bfps/DNS.py
+++ b/bfps/DNS.py
@@ -233,76 +233,102 @@ class DNS(_code):
"""
if len(list(self.statistics.keys())) > 0:
return None
- self.read_parameters()
- with self.get_data_file() as data_file:
- if 'moments' not in data_file['statistics'].keys():
- return None
- iter0 = min((data_file['statistics/moments/velocity'].shape[0] *
- self.parameters['niter_stat']-1),
- iter0)
- if type(iter1) == type(None):
- iter1 = data_file['iteration'].value
- else:
- iter1 = min(data_file['iteration'].value, iter1)
- 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
- self.statistics['kM'] = data_file['kspace/kM'].value
- self.statistics['dk'] = data_file['kspace/dk'].value
- computation_needed = True
- pp_file = h5py.File(self.get_postprocess_file_name(), 'a')
- if 'ii0' in pp_file.keys():
- computation_needed = not (ii0 == pp_file['ii0'].value and
- ii1 == pp_file['ii1'].value)
- if computation_needed:
- for k in ['t', 'vel_max(t)', 'renergy(t)',
- 'energy(t)', 'enstrophy(t)',
- 'energy(k)', 'enstrophy(k)',
- 'energy(t, k)',
- 'enstrophy(t, k)',
+ if not os.path.exists(self.get_data_file_name()):
+ if os.path.exists(self.get_cache_file_name()):
+ self.read_parameters(fname = self.get_cache_file_name())
+ with self.get_cache_file() as pp_file:
+ for k in ['t',
+ 'energy(t)',
+ 'energy(k)',
+ 'enstrophy(t)',
+ 'enstrophy(k)',
'R_ij(t)',
- 'ii0', 'ii1', 'iter0', 'iter1']:
+ 'vel_max(t)',
+ 'renergy(t)']:
if k in pp_file.keys():
- del pp_file[k]
- if computation_needed:
- pp_file['iter0'] = iter0
- pp_file['iter1'] = iter1
- pp_file['ii0'] = ii0
- pp_file['ii1'] = ii1
- 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)
- 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'])
- 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)*(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',
- '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()
+ self.statistics[k] = pp_file[k].value
+ self.statistics['kM'] = pp_file['kspace/kM'].value
+ self.statistics['dk'] = pp_file['kspace/dk'].value
+ self.statistics['kshell'] = pp_file['kspace/kshell'].value
+ self.statistics['nshell'] = pp_file['kspace/nshell'].value
+ else:
+ self.read_parameters()
+ with self.get_data_file() as data_file:
+ if 'moments' not in data_file['statistics'].keys():
+ return None
+ iter0 = min((data_file['statistics/moments/velocity'].shape[0] *
+ self.parameters['niter_stat']-1),
+ iter0)
+ if type(iter1) == type(None):
+ iter1 = data_file['iteration'].value
+ else:
+ iter1 = min(data_file['iteration'].value, iter1)
+ 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
+ self.statistics['kM'] = data_file['kspace/kM'].value
+ self.statistics['dk'] = data_file['kspace/dk'].value
+ computation_needed = True
+ pp_file = h5py.File(self.get_postprocess_file_name(), 'a')
+ if not ('parameters' in pp_file.keys()):
+ data_file.copy('parameters', pp_file)
+ data_file.copy('kspace', pp_file)
+ if 'ii0' in pp_file.keys():
+ computation_needed = not (ii0 == pp_file['ii0'].value and
+ ii1 == pp_file['ii1'].value)
+ if computation_needed:
+ 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
+ pp_file['ii0'] = ii0
+ pp_file['ii1'] = ii1
+ 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)
+ 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'])
+ 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)*(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
+ for k in ['t',
+ '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
+ # sanity check --- Parseval theorem check
+ assert(np.max(np.abs(
+ self.statistics['renergy(t)'] -
+ self.statistics['energy(t)']) / self.statistics['energy(t)']) < 1e-5)
+ self.compute_time_averages()
return None
def compute_Reynolds_stress_invariants(
self):
diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index 9c12ca201b928566545405c6ecaea12fb9da553d..c30adbe2ec41dac86993399a0235a18d20820269 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -584,79 +584,105 @@ class NavierStokes(_fluid_particle_base):
"""
if len(list(self.statistics.keys())) > 0:
return None
- self.read_parameters()
- with self.get_data_file() as data_file:
- if 'moments' not in data_file['statistics'].keys():
- return None
- iter0 = min((data_file['statistics/moments/velocity'].shape[0] *
- self.parameters['niter_stat']-1),
- iter0)
- if type(iter1) == type(None):
- iter1 = data_file['iteration'].value
- else:
- iter1 = min(data_file['iteration'].value, iter1)
- 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
- self.statistics['kM'] = data_file['kspace/kM'].value
- self.statistics['dk'] = data_file['kspace/dk'].value
- computation_needed = True
- pp_file = h5py.File(self.get_postprocess_file_name(), 'a')
- if 'ii0' in pp_file.keys():
- computation_needed = not (ii0 == pp_file['ii0'].value and
- ii1 == pp_file['ii1'].value)
- if computation_needed:
- for k in ['t', 'vel_max(t)', 'renergy(t)',
- 'energy(t)', 'enstrophy(t)',
- 'energy(k)', 'enstrophy(k)',
- 'energy(t, k)',
- 'enstrophy(t, k)',
+ if not os.path.exists(self.get_data_file_name()):
+ if os.path.exists(self.get_cache_file_name()):
+ self.read_parameters(fname = self.get_cache_file_name())
+ with self.get_cache_file() as pp_file:
+ for k in ['t',
+ 'energy(t)',
+ 'energy(k)',
+ 'enstrophy(t)',
+ 'enstrophy(k)',
'R_ij(t)',
- 'ii0', 'ii1', 'iter0', 'iter1']:
+ 'vel_max(t)',
+ 'renergy(t)']:
if k in pp_file.keys():
- del pp_file[k]
- if computation_needed:
- pp_file['iter0'] = iter0
- pp_file['iter1'] = iter1
- pp_file['ii0'] = ii0
- pp_file['ii1'] = ii1
- 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)
- 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'])
- 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)*(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():
- pp_file['mean_trS2(t)'] = data_file['statistics/moments/trS2_Q_R'][:, 1, 0]
- for k in ['t',
- 'energy(t)',
- 'energy(k)',
- 'enstrophy(t)',
- 'enstrophy(k)',
- 'R_ij(t)',
- 'vel_max(t)',
- 'renergy(t)',
- 'mean_trS2(t)']:
- if k in pp_file.keys():
- self.statistics[k] = pp_file[k].value
- self.compute_time_averages()
+ self.statistics[k] = pp_file[k].value
+ self.statistics['kM'] = pp_file['kspace/kM'].value
+ self.statistics['dk'] = pp_file['kspace/dk'].value
+ self.statistics['kshell'] = pp_file['kspace/kshell'].value
+ self.statistics['nshell'] = pp_file['kspace/nshell'].value
+ else:
+ self.read_parameters()
+ with self.get_data_file() as data_file:
+ if 'moments' not in data_file['statistics'].keys():
+ return None
+ iter0 = min((data_file['statistics/moments/velocity'].shape[0] *
+ self.parameters['niter_stat']-1),
+ iter0)
+ if type(iter1) == type(None):
+ iter1 = data_file['iteration'].value
+ else:
+ iter1 = min(data_file['iteration'].value, iter1)
+ 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
+ self.statistics['kM'] = data_file['kspace/kM'].value
+ self.statistics['dk'] = data_file['kspace/dk'].value
+ computation_needed = True
+ pp_file = h5py.File(self.get_postprocess_file_name(), 'a')
+ if not ('parameters' in pp_file.keys()):
+ data_file.copy('parameters', pp_file)
+ data_file.copy('kspace', pp_file)
+ if 'ii0' in pp_file.keys():
+ computation_needed = not (ii0 == pp_file['ii0'].value and
+ ii1 == pp_file['ii1'].value)
+ if computation_needed:
+ 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
+ pp_file['ii0'] = ii0
+ pp_file['ii1'] = ii1
+ 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)
+ 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'])
+ 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)*(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():
+ pp_file['mean_trS2(t)'] = data_file['statistics/moments/trS2_Q_R'][:, 1, 0]
+ for k in ['t',
+ 'energy(t)',
+ 'energy(k)',
+ 'enstrophy(t)',
+ 'enstrophy(k)',
+ 'R_ij(t)',
+ 'vel_max(t)',
+ 'renergy(t)',
+ 'mean_trS2(t)']:
+ if k in pp_file.keys():
+ self.statistics[k] = pp_file[k].value
+ # sanity check --- Parseval theorem check
+ assert(np.max(np.abs(
+ self.statistics['renergy(t)'] -
+ self.statistics['energy(t)']) / self.statistics['energy(t)']) < 1e-5)
+ self.compute_time_averages()
return None
def compute_Reynolds_stress_invariants(
self):
diff --git a/bfps/_base.py b/bfps/_base.py
index 2ea8fbbbf139cfdef6d5678b2def4aeef76180fc..575067272f5e1ef171eeb6a374fbf52a732ebc61 100644
--- a/bfps/_base.py
+++ b/bfps/_base.py
@@ -224,8 +224,10 @@ class _base(object):
ofile[group + '/' + k][...] = parameters[k]
ofile.close()
return None
- def read_parameters(self):
- with h5py.File(os.path.join(self.work_dir, self.simname + '.h5'), 'r') as data_file:
+ def read_parameters(self, fname = None):
+ if type(fname) == type(None):
+ fname = os.path.join(self.work_dir, self.simname + '.h5')
+ with h5py.File(fname, 'r') as data_file:
for k in data_file['parameters'].keys():
if k in self.parameters.keys():
if type(self.parameters[k]) in [int, str, float]: