diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index 920b8413a95b5929255a5f06f4e1afb6a4c6cbf0..445dda463ae4fc3bd48ce37dd09fe2f15a53b766 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -119,80 +119,76 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
self.fluid_includes += '#include "fftw_tools.hpp"\n'
self.stat_src += """
//begincpp
+ double *velocity_moments = new double[10*4];
+ double *vorticity_moments = new double[10*4];
+ ptrdiff_t *hist_velocity = new ptrdiff_t[histogram_bins*4];
+ ptrdiff_t *hist_vorticity = new ptrdiff_t[histogram_bins*4];
+ double max_estimates[4];
+ fs->compute_velocity(fs->cvorticity);
+ double *spec_velocity = new double[fs->nshells*9];
+ double *spec_vorticity = new double[fs->nshells*9];
+ fs->cospectrum(fs->cvelocity, fs->cvelocity, spec_velocity);
+ fs->cospectrum(fs->cvorticity, fs->cvorticity, spec_vorticity);
//endcpp
"""
- self.stat_src += """
- //begincpp
- double *velocity_moments = new double[10*4];
- double *vorticity_moments = new double[10*4];
- ptrdiff_t *hist_velocity = new ptrdiff_t[histogram_bins*4];
- ptrdiff_t *hist_vorticity = new ptrdiff_t[histogram_bins*4];
- double max_estimates[4];
- fs->compute_velocity(fs->cvorticity);
- double *spec_velocity = new double[fs->nshells*9];
- double *spec_vorticity = new double[fs->nshells*9];
- fs->cospectrum(fs->cvelocity, fs->cvelocity, spec_velocity);
- fs->cospectrum(fs->cvorticity, fs->cvorticity, spec_vorticity);
- //endcpp
- """
if self.QR_stats_on:
self.stat_src += """
//begincpp
- double *trS2_Q_R_moments = new double[10*3];
- double *gradu_moments = new double[10*9];
- ptrdiff_t *hist_trS2_Q_R = new ptrdiff_t[histogram_bins*3];
- ptrdiff_t *hist_gradu = new ptrdiff_t[histogram_bins*9];
- ptrdiff_t *hist_QR2D = new ptrdiff_t[QR2D_histogram_bins*QR2D_histogram_bins];
- double trS2QR_max_estimates[3];
- double gradu_max_estimates[9];
- trS2QR_max_estimates[0] = max_trS2_estimate;
- trS2QR_max_estimates[1] = max_Q_estimate;
- trS2QR_max_estimates[2] = max_R_estimate;
- std::fill_n(gradu_max_estimates, 9, sqrt(3*max_trS2_estimate));
- fs->compute_gradient_statistics(
- fs->cvelocity,
- gradu_moments,
- trS2_Q_R_moments,
- hist_gradu,
- hist_trS2_Q_R,
- hist_QR2D,
- trS2QR_max_estimates,
- gradu_max_estimates,
- histogram_bins,
- QR2D_histogram_bins);
- //endcpp
- """
+ double *trS2_Q_R_moments = new double[10*3];
+ double *gradu_moments = new double[10*9];
+ ptrdiff_t *hist_trS2_Q_R = new ptrdiff_t[histogram_bins*3];
+ ptrdiff_t *hist_gradu = new ptrdiff_t[histogram_bins*9];
+ ptrdiff_t *hist_QR2D = new ptrdiff_t[QR2D_histogram_bins*QR2D_histogram_bins];
+ double trS2QR_max_estimates[3];
+ double gradu_max_estimates[9];
+ trS2QR_max_estimates[0] = max_trS2_estimate;
+ trS2QR_max_estimates[1] = max_Q_estimate;
+ trS2QR_max_estimates[2] = max_R_estimate;
+ std::fill_n(gradu_max_estimates, 9, sqrt(3*max_trS2_estimate));
+ fs->compute_gradient_statistics(
+ fs->cvelocity,
+ gradu_moments,
+ trS2_Q_R_moments,
+ hist_gradu,
+ hist_trS2_Q_R,
+ hist_QR2D,
+ trS2QR_max_estimates,
+ gradu_max_estimates,
+ histogram_bins,
+ QR2D_histogram_bins);
+ //endcpp
+ """
self.stat_src += """
//begincpp
- fs->ift_velocity();
- max_estimates[0] = max_velocity_estimate/sqrt(3);
- max_estimates[1] = max_estimates[0];
- max_estimates[2] = max_estimates[0];
- max_estimates[3] = max_velocity_estimate;
- fs->compute_rspace_stats4(fs->rvelocity,
- velocity_moments,
- hist_velocity,
- max_estimates,
- histogram_bins);
- fs->ift_vorticity();
- max_estimates[0] = max_vorticity_estimate/sqrt(3);
- max_estimates[1] = max_estimates[0];
- max_estimates[2] = max_estimates[0];
- max_estimates[3] = max_vorticity_estimate;
- fs->compute_rspace_stats4(fs->rvorticity,
- vorticity_moments,
- hist_vorticity,
- max_estimates,
- histogram_bins);
- if (fs->cd->myrank == 0)
- {{
- hid_t Cdset, wspace, mspace;
- int ndims;
- hsize_t count[4], offset[4], dims[4];
- offset[0] = fs->iteration/niter_stat;
- offset[1] = 0;
- offset[2] = 0;
- offset[3] = 0;
+ fs->ift_velocity();
+ max_estimates[0] = max_velocity_estimate/sqrt(3);
+ max_estimates[1] = max_estimates[0];
+ max_estimates[2] = max_estimates[0];
+ max_estimates[3] = max_velocity_estimate;
+ fs->compute_rspace_stats4(fs->rvelocity,
+ velocity_moments,
+ hist_velocity,
+ max_estimates,
+ histogram_bins);
+ fs->ift_vorticity();
+ max_estimates[0] = max_vorticity_estimate/sqrt(3);
+ max_estimates[1] = max_estimates[0];
+ max_estimates[2] = max_estimates[0];
+ max_estimates[3] = max_vorticity_estimate;
+ fs->compute_rspace_stats4(fs->rvorticity,
+ vorticity_moments,
+ hist_vorticity,
+ max_estimates,
+ histogram_bins);
+ if (fs->cd->myrank == 0)
+ {{
+ hid_t Cdset, wspace, mspace;
+ int ndims;
+ hsize_t count[4], offset[4], dims[4];
+ offset[0] = fs->iteration/niter_stat;
+ offset[1] = 0;
+ offset[2] = 0;
+ offset[3] = 0;
//endcpp
""".format(self.C_dtype)
if self.dtype == np.float32:
@@ -300,23 +296,23 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
""")
self.stat_src += """
//begincpp
- }
- delete[] spec_velocity;
- delete[] spec_vorticity;
- delete[] velocity_moments;
- delete[] vorticity_moments;
- delete[] hist_velocity;
- delete[] hist_vorticity;
+ }
+ delete[] spec_velocity;
+ delete[] spec_vorticity;
+ delete[] velocity_moments;
+ delete[] vorticity_moments;
+ delete[] hist_velocity;
+ delete[] hist_vorticity;
//endcpp
"""
if self.QR_stats_on:
self.stat_src += """
//begincpp
- delete[] trS2_Q_R_moments;
- delete[] gradu_moments;
- delete[] hist_trS2_Q_R;
- delete[] hist_gradu;
- delete[] hist_QR2D;
+ delete[] trS2_Q_R_moments;
+ delete[] gradu_moments;
+ delete[] hist_trS2_Q_R;
+ delete[] hist_gradu;
+ delete[] hist_QR2D;
//endcpp
"""
return None
@@ -399,7 +395,8 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
field_name = 'fs->rvelocity'):
name = quantity + '_' + name
self.fluid_includes += '#include "rFFTW_interpolator.hpp"\n'
- self.fluid_variables += 'rFFTW_interpolator <{0}, {1}> *{2};\n'.format(self.C_dtype, neighbours, name)
+ self.fluid_variables += 'rFFTW_interpolator <{0}, {1}> *{2};\n'.format(
+ self.C_dtype, neighbours, name)
self.parameters[name + '_type'] = interp_type
self.parameters[name + '_neighbours'] = neighbours
if interp_type == 'spline':
@@ -434,7 +431,8 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
field_class = 'rFFTW_interpolator',
acc_field_name = 'rFFTW_acc'):
self.fluid_includes += '#include "{0}.hpp"\n'.format(field_class)
- self.fluid_variables += field_class + '<{0}, {1}> *vel_{2}, *acc_{2};\n'.format(self.C_dtype, neighbours, name)
+ self.fluid_variables += field_class + '<{0}, {1}> *vel_{2}, *acc_{2};\n'.format(
+ self.C_dtype, neighbours, name)
self.parameters[name + '_type'] = interp_type
self.parameters[name + '_neighbours'] = neighbours
if interp_type == 'spline':
@@ -442,13 +440,22 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
beta_name = 'beta_n{0}_m{1}'.format(neighbours, smoothness)
elif interp_type == 'Lagrange':
beta_name = 'beta_Lagrange_n{0}'.format(neighbours)
- self.fluid_start += ('vel_{0} = new {1}<{2}, {3}>(fs, {4}, fs->rvelocity);\n' +
- 'acc_{0} = new {1}<{2}, {3}>(fs, {4}, {5});\n').format(name,
- field_class,
- self.C_dtype,
- neighbours,
- beta_name,
- acc_field_name)
+ if field_class == 'rFFTW_interpolator':
+ self.fluid_start += ('vel_{0} = new {1}<{2}, {3}>(fs, {4}, fs->rvelocity);\n' +
+ 'acc_{0} = new {1}<{2}, {3}>(fs, {4}, {5});\n').format(name,
+ field_class,
+ self.C_dtype,
+ neighbours,
+ beta_name,
+ acc_field_name)
+ elif field_class == 'interpolator':
+ self.fluid_start += ('vel_{0} = new {1}<{2}, {3}>(fs, {4});\n' +
+ 'acc_{0} = new {1}<{2}, {3}>(fs, {4});\n').format(name,
+ field_class,
+ self.C_dtype,
+ neighbours,
+ beta_name,
+ acc_field_name)
self.fluid_end += ('delete vel_{0};\n' +
'delete acc_{0};\n').format(name)
update_fields = 'fs->compute_velocity(fs->cvorticity);\n'
@@ -474,9 +481,11 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
integration_steps = 2
neighbours = self.parameters[fields_name + '_neighbours']
self.parameters['tracers{0}_field'.format(self.particle_species)] = fields_name
- self.parameters['tracers{0}_integration_method'.format(self.particle_species)] = integration_method
+ self.parameters['tracers{0}_integration_method'.format(
+ self.particle_species)] = integration_method
self.parameters['tracers{0}_kcut'.format(self.particle_species)] = kcut
- self.parameters['tracers{0}_integration_steps'.format(self.particle_species)] = integration_steps
+ self.parameters['tracers{0}_integration_steps'.format(
+ self.particle_species)] = integration_steps
self.file_datasets_grow += """
//begincpp
temp_string = (std::string("/particles/") +
@@ -491,48 +500,48 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
elif self.dtype == np.float64:
FFTW = 'fftw'
output_vel_acc = """
- //begincpp
- {{
- double *acceleration = new double[ps{0}->array_size];
- double *velocity = new double[ps{0}->array_size];
- ps{0}->sample_vec_field(vel_{1}, velocity);
- ps{0}->sample_vec_field(acc_{1}, acceleration);
- if (ps{0}->myrank == 0)
- {{
- //VELOCITY begin
- std::string temp_string = (std::string("/particles/") +
- std::string(ps{0}->name) +
- std::string("/velocity"));
- hid_t Cdset = H5Dopen(stat_file, temp_string.c_str(), H5P_DEFAULT);
- hid_t mspace, wspace;
- int ndims;
- hsize_t count[3], offset[3];
- wspace = H5Dget_space(Cdset);
- ndims = H5Sget_simple_extent_dims(wspace, count, NULL);
- count[0] = 1;
- offset[0] = ps{0}->iteration / ps{0}->traj_skip;
- offset[1] = 0;
- offset[2] = 0;
- mspace = H5Screate_simple(ndims, count, NULL);
- H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL);
- H5Dwrite(Cdset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, velocity);
- H5Dclose(Cdset);
- //VELOCITY end
- //ACCELERATION begin
- temp_string = (std::string("/particles/") +
- std::string(ps{0}->name) +
- std::string("/acceleration"));
- Cdset = H5Dopen(stat_file, temp_string.c_str(), H5P_DEFAULT);
- H5Dwrite(Cdset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, acceleration);
- H5Sclose(mspace);
- H5Sclose(wspace);
- H5Dclose(Cdset);
- //ACCELERATION end
- }}
- delete[] acceleration;
- delete[] velocity;
- }}
- //endcpp
+ //begincpp
+ {{
+ double *acceleration = new double[ps{0}->array_size];
+ double *velocity = new double[ps{0}->array_size];
+ ps{0}->sample_vec_field(vel_{1}, velocity);
+ ps{0}->sample_vec_field(acc_{1}, acceleration);
+ if (ps{0}->myrank == 0)
+ {{
+ //VELOCITY begin
+ std::string temp_string = (std::string("/particles/") +
+ std::string(ps{0}->name) +
+ std::string("/velocity"));
+ hid_t Cdset = H5Dopen(stat_file, temp_string.c_str(), H5P_DEFAULT);
+ hid_t mspace, wspace;
+ int ndims;
+ hsize_t count[3], offset[3];
+ wspace = H5Dget_space(Cdset);
+ ndims = H5Sget_simple_extent_dims(wspace, count, NULL);
+ count[0] = 1;
+ offset[0] = ps{0}->iteration / ps{0}->traj_skip;
+ offset[1] = 0;
+ offset[2] = 0;
+ mspace = H5Screate_simple(ndims, count, NULL);
+ H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL);
+ H5Dwrite(Cdset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, velocity);
+ H5Dclose(Cdset);
+ //VELOCITY end
+ //ACCELERATION begin
+ temp_string = (std::string("/particles/") +
+ std::string(ps{0}->name) +
+ std::string("/acceleration"));
+ Cdset = H5Dopen(stat_file, temp_string.c_str(), H5P_DEFAULT);
+ H5Dwrite(Cdset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, acceleration);
+ H5Sclose(mspace);
+ H5Sclose(wspace);
+ H5Dclose(Cdset);
+ //ACCELERATION end
+ }}
+ delete[] acceleration;
+ delete[] velocity;
+ }}
+ //endcpp
""".format(self.particle_species, fields_name)
self.particle_start += 'sprintf(fname, "tracers{0}");\n'.format(self.particle_species)
self.particle_end += ('ps{0}->write(stat_file);\n' +
@@ -552,7 +561,11 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
'fname, fs, vel_{4},\n' +
'nparticles,\n' +
'niter_part, tracers{0}_integration_steps);\n').format(
- self.particle_species, self.C_dtype, multistep, neighbours, fields_name)
+ self.particle_species,
+ self.C_dtype,
+ multistep,
+ neighbours,
+ fields_name)
else:
self.particle_variables += 'rFFTW_particles<VELOCITY_TRACER, {0}, {1}> *ps{2};\n'.format(
self.C_dtype,
@@ -562,7 +575,10 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
'fname, vel_{3},\n' +
'nparticles,\n' +
'niter_part, tracers{0}_integration_steps);\n').format(
- self.particle_species, self.C_dtype, neighbours, fields_name)
+ self.particle_species,
+ self.C_dtype,
+ neighbours,
+ fields_name)
self.particle_start += ('ps{0}->dt = dt;\n' +
'ps{0}->iteration = iteration;\n' +
'ps{0}->read(stat_file);\n').format(self.particle_species)
@@ -601,7 +617,8 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
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 = min((data_file['statistics/moments/velocity'].shape[0] *
+ self.parameters['niter_stat']-1),
iter0)
if type(iter1) == type(None):
iter1 = data_file['iteration'].value
@@ -613,9 +630,9 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
self.statistics['kM'] = data_file['kspace/kM'].value
self.statistics['dk'] = data_file['kspace/dk'].value
if self.particle_species > 0:
- self.trajectories = [
- data_file['particles/' + key + '/state'][
- iter0//self.parameters['niter_part']:iter1//self.parameters['niter_part']+1]
+ self.trajectories = [data_file['particles/' + key + '/state'][
+ iter0//self.parameters['niter_part'] :
+ iter1//self.parameters['niter_part']+1]
for key in data_file['particles'].keys()]
computation_needed = True
pp_file = h5py.File(self.get_postprocess_file_name(), 'a')
@@ -660,7 +677,8 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
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)) *
+ 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',
@@ -687,7 +705,8 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
self.statistics['Rlambda' + suffix] = (self.statistics['Uint' + suffix] *
self.statistics['lambda' + suffix] /
self.parameters['nu'])
- self.statistics['kMeta' + suffix] = self.statistics['kM']*self.statistics['etaK' + suffix]
+ self.statistics['kMeta' + suffix] = (self.statistics['kM'] *
+ self.statistics['etaK' + suffix])
if self.parameters['dealias_type'] == 1:
self.statistics['kMeta' + suffix] *= 0.8
self.statistics['Tint'] = self.statistics['Lint'] / self.statistics['Uint']