diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index f3fa39d04c76ac5b400b481ecf423095b4135bcd..70a4426d6f81c0610aff559b7cb8620194f74a0a 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
@@ -378,16 +374,27 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
'}\n' +
'delete fs;\n')
return None
- def add_particle_fields(
+ def add_3D_rFFTW_field(
+ self,
+ name = 'rFFTW_acc'):
+ if self.dtype == np.float32:
+ FFTW = 'fftwf'
+ elif self.dtype == np.float64:
+ FFTW = 'fftw'
+ self.fluid_variables += '{0} *{1};\n'.format(self.C_dtype, name)
+ self.fluid_start += '{0} = {1}_alloc_real(2*fs->cd->local_size);\n'.format(name, FFTW)
+ self.fluid_end += '{0}_free({1});\n'.format(FFTW, name)
+ return None
+ def add_interpolator(
self,
interp_type = 'spline',
- kcut = None,
neighbours = 1,
smoothness = 1,
- name = 'particle_field',
- field_class = 'rFFTW_interpolator'):
- 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)
+ name = 'field_interpolator',
+ field_name = 'fs->rvelocity'):
+ self.fluid_includes += '#include "rFFTW_interpolator.hpp"\n'
+ 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':
@@ -395,41 +402,66 @@ 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});\n' +
- 'acc_{0} = new {1}<{2}, {3}>(fs, {4});\n').format(name,
- field_class,
- self.C_dtype,
- neighbours,
- beta_name)
- self.fluid_end += ('delete vel_{0};\n' +
- 'delete acc_{0};\n').format(name)
- update_fields = 'fs->compute_velocity(fs->cvorticity);\n'
- if not type(kcut) == type(None):
- update_fields += 'fs->low_pass_Fourier(fs->cvelocity, 3, {0});\n'.format(kcut)
- update_fields += ('fs->ift_velocity();\n' +
- 'vel_{0}->read_rFFTW(fs->rvelocity);\n' +
- 'fs->compute_Lagrangian_acceleration(acc_{0}->field);\n').format(name)
- self.fluid_start += update_fields
- self.fluid_loop += update_fields
+ self.fluid_start += '{0} = new rFFTW_interpolator<{1}, {2}>(fs, {3}, {4});\n'.format(
+ name,
+ self.C_dtype,
+ neighbours,
+ beta_name,
+ field_name)
+ self.fluid_end += 'delete {0};\n'.format(name)
return None
def add_particles(
self,
- integration_method = 'AdamsBashforth',
integration_steps = 2,
- kcut = 'fs->kM',
+ kcut = None,
+ interpolator = 'field_interpolator',
frozen_particles = False,
- fields_name = None,
- particle_class = 'rFFTW_particles'):
- if integration_method == 'cRK4':
- integration_steps = 4
- elif integration_method == 'Heun':
- 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}_kcut'.format(self.particle_species)] = kcut
- self.parameters['tracers{0}_integration_steps'.format(self.particle_species)] = integration_steps
- self.file_datasets_grow += """
+ acc_name = None):
+ """Adds code for tracking a series of particle species, each
+ consisting of `nparticles` particles.
+
+ :type integration_steps: int, list of int
+ :type kcut: None (default), str, list of str
+ :type interpolator: str, list of str
+ :type frozen_particles: bool
+ :type acc_name: str
+
+ .. warning :: if not None, kcut must be a list of decreasing
+ wavenumbers, since filtering is done on the same field...
+ """
+ if self.dtype == np.float32:
+ FFTW = 'fftwf'
+ elif self.dtype == np.float64:
+ FFTW = 'fftw'
+ s0 = self.particle_species
+ if type(integration_steps) == int:
+ integration_steps = [integration_steps]
+ if type(kcut) == str:
+ kcut = [kcut]
+ if type(interpolator) == str:
+ interpolator = [interpolator]
+ nspecies = max(len(integration_steps), len(interpolator))
+ if type(kcut) == list:
+ nspecies = max(nspecies, len(kcut))
+ if len(integration_steps) == 1:
+ integration_steps = [integration_steps[0] for s in range(nspecies)]
+ if len(interpolator) == 1:
+ interpolator = [interpolator[0] for s in range(nspecies)]
+ if type(kcut) == list:
+ if len(kcut) == 1:
+ kcut = [kcut[0] for s in range(nspecies)]
+ assert(len(integration_steps) == nspecies)
+ assert(len(interpolator) == nspecies)
+ if type(kcut) == list:
+ assert(len(kcut) == nspecies)
+ for s in range(nspecies):
+ neighbours = self.parameters[interpolator[s] + '_neighbours']
+ if type(kcut) == list:
+ self.parameters['tracers{0}_kcut'.format(s0 + s)] = kcut[s]
+ self.parameters['tracers{0}_acc_on'.format(s0 + s)] = not type(acc_name) == type(None)
+ self.parameters['tracers{0}_interpolator'.format(s0 + s)] = interpolator[s]
+ self.parameters['tracers{0}_integration_steps'.format(s0 + s)] = integration_steps[s]
+ self.file_datasets_grow += """
//begincpp
temp_string = (std::string("/particles/") +
std::string(ps{0}->name));
@@ -437,110 +469,118 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
grow_particle_datasets(group, temp_string.c_str(), NULL, NULL);
H5Gclose(group);
//endcpp
- """.format(self.particle_species)
- if self.dtype == np.float32:
- FFTW = 'fftwf'
- 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
- """.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' +
- 'delete ps{0};\n').format(self.particle_species)
- self.particle_includes += '#include "{0}.hpp"\n'.format(particle_class)
- if particle_class == 'particles':
- if integration_method == 'AdamsBashforth':
- multistep = 'true'
- else:
- multistep = 'false'
- self.particle_variables += 'particles<VELOCITY_TRACER, {0}, {1}, {2}> *ps{3};\n'.format(
- self.C_dtype,
- multistep,
- neighbours,
- self.particle_species)
- self.particle_start += ('ps{0} = new particles<VELOCITY_TRACER, {1}, {2},{3}>(\n' +
- '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)
+ """.format(s0 + s)
+
+ #### code that outputs statistics
+ output_vel_acc = '{\n'
+ # array for putting sampled velocity in
+ # must compute velocity, just in case it was messed up by some
+ # other particle species before the stats
+ output_vel_acc += ('double *velocity = new double[3*nparticles];\n' +
+ 'fs->compute_velocity(fs->cvorticity);\n')
+ if not type(kcut) == list:
+ output_vel_acc += 'fs->ift_velocity();\n'
+ if not type(acc_name) == type(None):
+ # array for putting sampled acceleration in
+ # must compute acceleration
+ output_vel_acc += 'double *acceleration = new double[3*nparticles];\n'
+ output_vel_acc += 'fs->compute_Lagrangian_acceleration({0});\n'.format(acc_name)
+ for s in range(nspecies):
+ if type(kcut) == list:
+ output_vel_acc += 'fs->low_pass_Fourier(fs->cvelocity, 3, {0});\n'.format(kcut[s])
+ output_vel_acc += 'fs->ift_velocity();\n'
+ output_vel_acc += """
+ {0}->field = fs->rvelocity;
+ ps{1}->sample_vec_field({0}, velocity);
+ """.format(interpolator[s], s0 + s)
+ if not type(acc_name) == type(None):
+ output_vel_acc += """
+ {0}->field = {1};
+ ps{2}->sample_vec_field({0}, acceleration);
+ """.format(interpolator[s], acc_name, s0 + s)
+ output_vel_acc += """
+ if (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\n""".format(s0 + s)
+ if not type(acc_name) == type(None):
+ output_vel_acc += """
+ //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\n""".format(s0 + s)
+ output_vel_acc += '}\n'
+ output_vel_acc += 'delete[] velocity;\n'
+ if not type(acc_name) == type(None):
+ output_vel_acc += 'delete[] acceleration;\n'
+ output_vel_acc += '}\n'
+
+ #### initialize, stepping and finalize code
+ if not type(kcut) == list:
+ update_fields = ('fs->compute_velocity(fs->cvorticity);\n' +
+ 'fs->ift_velocity();\n')
+ self.particle_start += update_fields
+ self.particle_loop += update_fields
else:
+ self.particle_loop += 'fs->compute_velocity(fs->cvorticity);\n'
+ self.particle_includes += '#include "rFFTW_particles.hpp"\n'
+ self.particle_stat_src += (
+ 'if (ps0->iteration % niter_part == 0)\n' +
+ '{\n')
+ for s in range(nspecies):
+ self.particle_start += 'sprintf(fname, "tracers{0}");\n'.format(s0 + s)
+ self.particle_end += ('ps{0}->write(stat_file);\n' +
+ 'delete ps{0};\n').format(s0 + s)
self.particle_variables += 'rFFTW_particles<VELOCITY_TRACER, {0}, {1}> *ps{2};\n'.format(
self.C_dtype,
neighbours,
- self.particle_species)
+ s0 + s)
self.particle_start += ('ps{0} = new rFFTW_particles<VELOCITY_TRACER, {1}, {2}>(\n' +
- 'fname, vel_{3},\n' +
+ 'fname, {3},\n' +
'nparticles,\n' +
'niter_part, tracers{0}_integration_steps);\n').format(
- 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)
+ s0 + s,
+ self.C_dtype,
+ neighbours,
+ interpolator[s])
+ self.particle_start += ('ps{0}->dt = dt;\n' +
+ 'ps{0}->iteration = iteration;\n' +
+ 'ps{0}->read(stat_file);\n').format(s0 + s)
+ if not frozen_particles:
+ if type(kcut) == list:
+ update_field = ('fs->low_pass_Fourier(fs->cvelocity, 3, {0});\n'.format(kcut[s]) +
+ 'fs->ift_velocity();\n')
+ self.particle_loop += update_field
+ self.particle_loop += '{0}->field = fs->rvelocity;\n'.format(interpolator[s])
+ self.particle_loop += 'ps{0}->step();\n'.format(s0 + s)
+ self.particle_stat_src += 'ps{0}->write(stat_file, false);\n'.format(s0 + s)
self.particle_start += output_vel_acc
- if not frozen_particles:
- if particle_class == 'particles':
- if integration_method == 'AdamsBashforth':
- self.particle_loop += 'ps{0}->AdamsBashforth((ps{0}->iteration < ps{0}->integration_steps) ? ps{0}->iteration+1 : ps{0}->integration_steps);\n'.format(self.particle_species)
- elif integration_method == 'Euler':
- self.particle_loop += 'ps{0}->Euler();\n'.format(self.particle_species)
- elif integration_method == 'Heun':
- assert(integration_steps == 2)
- self.particle_loop += 'ps{0}->Heun();\n'.format(self.particle_species)
- elif integration_method == 'cRK4':
- assert(integration_steps == 4)
- self.particle_loop += 'ps{0}->cRK4();\n'.format(self.particle_species)
- self.particle_loop += 'ps{0}->iteration++;\n'.format(self.particle_species)
- self.particle_loop += 'ps{0}->synchronize();\n'.format(self.particle_species)
- elif particle_class == 'rFFTW_particles':
- self.particle_loop += 'ps{0}->step();\n'.format(self.particle_species)
- self.particle_stat_src += (
- ('if (ps{0}->iteration % niter_part == 0)\n' +
- '{{\n' +
- 'ps{0}->write(stat_file, false);\n').format(self.particle_species) +
- output_vel_acc + '}\n')
- self.particle_species += 1
+ self.particle_stat_src += output_vel_acc
+ self.particle_stat_src += '}\n'
+ self.particle_species += nspecies
return None
def get_data_file(self):
return h5py.File(os.path.join(self.work_dir, self.simname + '.h5'), 'r')
@@ -553,7 +593,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
@@ -565,9 +606,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')
@@ -612,7 +653,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',
@@ -639,11 +681,12 @@ 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])
+ if self.parameters['dealias_type'] == 1:
+ self.statistics['kMeta' + suffix] *= 0.8
self.statistics['Tint'] = self.statistics['Lint'] / self.statistics['Uint']
self.statistics['Taylor_microscale'] = self.statistics['lambda']
- self.statistics['kMeta'] = self.statistics['kM']*self.statistics['etaK']
- if self.parameters['dealias_type'] == 1:
- self.statistics['kMeta'] *= 0.8
return None
def set_plt_style(
self,
@@ -731,6 +774,51 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
dtype = np.int64,
compression = 'gzip')
return None
+ def add_particle_fields(
+ self,
+ interp_type = 'spline',
+ kcut = None,
+ neighbours = 1,
+ smoothness = 1,
+ name = 'particle_field',
+ 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.parameters[name + '_type'] = interp_type
+ self.parameters[name + '_neighbours'] = neighbours
+ if interp_type == 'spline':
+ self.parameters[name + '_smoothness'] = smoothness
+ beta_name = 'beta_n{0}_m{1}'.format(neighbours, smoothness)
+ elif interp_type == 'Lagrange':
+ beta_name = 'beta_Lagrange_n{0}'.format(neighbours)
+ 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'
+ if not type(kcut) == type(None):
+ update_fields += 'fs->low_pass_Fourier(fs->cvelocity, 3, {0});\n'.format(kcut)
+ update_fields += ('fs->ift_velocity();\n' +
+ 'fs->compute_Lagrangian_acceleration(acc_{0}->field);\n').format(name)
+ self.fluid_start += update_fields
+ self.fluid_loop += update_fields
+ return None
def launch(
opt,
diff --git a/bfps/cpp/interpolator.hpp b/bfps/cpp/interpolator.hpp
index 299fef4910693b878d3f2734636b43dbf8f98a6d..e083bd56c6c27c1cc048d0ecfdfb7ee416082651 100644
--- a/bfps/cpp/interpolator.hpp
+++ b/bfps/cpp/interpolator.hpp
@@ -27,23 +27,12 @@
#include "field_descriptor.hpp"
#include "fftw_tools.hpp"
#include "fluid_solver_base.hpp"
-#include "spline_n1.hpp"
-#include "spline_n2.hpp"
-#include "spline_n3.hpp"
-#include "spline_n4.hpp"
-#include "spline_n5.hpp"
-#include "spline_n6.hpp"
-#include "Lagrange_polys.hpp"
+#include "interpolator_base.hpp"
#ifndef INTERPOLATOR
#define INTERPOLATOR
-typedef void (*base_polynomial_values)(
- const int derivative,
- const double fraction,
- double *__restrict__ destination);
-
template <class rnumber, int interp_neighbours>
class interpolator
{
diff --git a/bfps/cpp/interpolator_base.hpp b/bfps/cpp/interpolator_base.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..69c2d4f27bb3b029fc5fd4c4e85d6160b23fe32c
--- /dev/null
+++ b/bfps/cpp/interpolator_base.hpp
@@ -0,0 +1,45 @@
+/**********************************************************************
+* *
+* Copyright 2015 Max Planck Institute *
+* for Dynamics and Self-Organization *
+* *
+* This file is part of bfps. *
+* *
+* bfps is free software: you can redistribute it and/or modify *
+* it under the terms of the GNU General Public License as published *
+* by the Free Software Foundation, either version 3 of the License, *
+* or (at your option) any later version. *
+* *
+* bfps is distributed in the hope that it will be useful, *
+* but WITHOUT ANY WARRANTY; without even the implied warranty of *
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+* GNU General Public License for more details. *
+* *
+* You should have received a copy of the GNU General Public License *
+* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* *
+* Contact: Cristian.Lalescu@ds.mpg.de *
+* *
+**********************************************************************/
+
+
+
+#include "spline_n1.hpp"
+#include "spline_n2.hpp"
+#include "spline_n3.hpp"
+#include "spline_n4.hpp"
+#include "spline_n5.hpp"
+#include "spline_n6.hpp"
+#include "Lagrange_polys.hpp"
+
+#ifndef INTERPOLATOR_BASE
+
+#define INTERPOLATOR_BASE
+
+typedef void (*base_polynomial_values)(
+ const int derivative,
+ const double fraction,
+ double *__restrict__ destination);
+
+#endif//INTERPOLATOR_BASE
+
diff --git a/bfps/cpp/particles.hpp b/bfps/cpp/particles.hpp
index 4b743a006c15c1c83b65d66b9fd8056caee6e2e9..6fa379c458f38d364765efa4d7e58bb25e9daff5 100644
--- a/bfps/cpp/particles.hpp
+++ b/bfps/cpp/particles.hpp
@@ -29,6 +29,7 @@
#include <iostream>
#include <hdf5.h>
#include "base.hpp"
+#include "particles_base.hpp"
#include "fluid_solver_base.hpp"
#include "interpolator.hpp"
@@ -36,9 +37,6 @@
#define PARTICLES
-/* particle types */
-enum particle_types {VELOCITY_TRACER};
-
template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
class particles
{
diff --git a/bfps/cpp/particles_base.hpp b/bfps/cpp/particles_base.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..cf61f38efadf0030a67a10ff94c21a92e79bdd64
--- /dev/null
+++ b/bfps/cpp/particles_base.hpp
@@ -0,0 +1,37 @@
+/**********************************************************************
+* *
+* Copyright 2015 Max Planck Institute *
+* for Dynamics and Self-Organization *
+* *
+* This file is part of bfps. *
+* *
+* bfps is free software: you can redistribute it and/or modify *
+* it under the terms of the GNU General Public License as published *
+* by the Free Software Foundation, either version 3 of the License, *
+* or (at your option) any later version. *
+* *
+* bfps is distributed in the hope that it will be useful, *
+* but WITHOUT ANY WARRANTY; without even the implied warranty of *
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+* GNU General Public License for more details. *
+* *
+* You should have received a copy of the GNU General Public License *
+* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* *
+* Contact: Cristian.Lalescu@ds.mpg.de *
+* *
+**********************************************************************/
+
+
+
+#include "interpolator_base.hpp"
+
+#ifndef PARTICLES_BASE
+
+#define PARTICLES_BASE
+
+/* particle types */
+enum particle_types {VELOCITY_TRACER};
+
+#endif//PARTICLES_BASE
+
diff --git a/bfps/cpp/rFFTW_interpolator.cpp b/bfps/cpp/rFFTW_interpolator.cpp
index a858eb104abf19837f7ccd30ae673caaf841a493..c701fbd64bba8e1c68ab2de953eae0795baead2f 100644
--- a/bfps/cpp/rFFTW_interpolator.cpp
+++ b/bfps/cpp/rFFTW_interpolator.cpp
@@ -32,15 +32,13 @@
template <class rnumber, int interp_neighbours>
rFFTW_interpolator<rnumber, interp_neighbours>::rFFTW_interpolator(
fluid_solver_base<rnumber> *fs,
- base_polynomial_values BETA_POLYS)
+ base_polynomial_values BETA_POLYS,
+ rnumber *FIELD)
{
this->descriptor = fs->rd;
this->field_size = 2*fs->cd->local_size;
this->compute_beta = BETA_POLYS;
- if (sizeof(rnumber) == 4)
- this->field = (rnumber*)((void*)fftwf_alloc_real(this->field_size));
- else if (sizeof(rnumber) == 8)
- this->field = (rnumber*)((void*)fftw_alloc_real(this->field_size));
+ this->field = FIELD;
// compute dx, dy, dz;
this->dx = 4*acos(0) / (fs->dkx*this->descriptor->sizes[2]);
@@ -59,24 +57,9 @@ rFFTW_interpolator<rnumber, interp_neighbours>::rFFTW_interpolator(
template <class rnumber, int interp_neighbours>
rFFTW_interpolator<rnumber, interp_neighbours>::~rFFTW_interpolator()
{
- if (sizeof(rnumber) == 4)
- fftwf_free((float*)((void*)this->field));
- else if (sizeof(rnumber) == 8)
- fftw_free((double*)((void*)this->field));
delete[] this->compute;
}
-template <class rnumber, int interp_neighbours>
-int rFFTW_interpolator<rnumber, interp_neighbours>::read_rFFTW(void *void_src)
-{
- rnumber *src = (rnumber*)void_src;
- /* do big copy of middle stuff */
- std::copy(src,
- src + this->field_size,
- this->field);
- return EXIT_SUCCESS;
-}
-
template <class rnumber, int interp_neighbours>
void rFFTW_interpolator<rnumber, interp_neighbours>::get_grid_coordinates(
const int nparticles,
diff --git a/bfps/cpp/rFFTW_interpolator.hpp b/bfps/cpp/rFFTW_interpolator.hpp
index 0864bd964b77b5b3527a58898b8c22da63c33c2e..3dccd7150d17bea612128bf6d70b37033d281db8 100644
--- a/bfps/cpp/rFFTW_interpolator.hpp
+++ b/bfps/cpp/rFFTW_interpolator.hpp
@@ -27,14 +27,7 @@
#include "field_descriptor.hpp"
#include "fftw_tools.hpp"
#include "fluid_solver_base.hpp"
-#include "spline_n1.hpp"
-#include "spline_n2.hpp"
-#include "spline_n3.hpp"
-#include "spline_n4.hpp"
-#include "spline_n5.hpp"
-#include "spline_n6.hpp"
-#include "Lagrange_polys.hpp"
-#include "interpolator.hpp"
+#include "interpolator_base.hpp"
#ifndef RFFTW_INTERPOLATOR
@@ -67,7 +60,8 @@ class rFFTW_interpolator
rFFTW_interpolator(
fluid_solver_base<rnumber> *FSOLVER,
- base_polynomial_values BETA_POLYS);
+ base_polynomial_values BETA_POLYS,
+ rnumber *FIELD_DATA);
~rFFTW_interpolator();
/* map real locations to grid coordinates */
@@ -90,7 +84,6 @@ class rFFTW_interpolator
const double *__restrict__ xx,
double *__restrict__ dest,
const int *deriv = NULL);
- int read_rFFTW(void *src);
};
#endif//RFFTW_INTERPOLATOR
diff --git a/bfps/cpp/rFFTW_particles.hpp b/bfps/cpp/rFFTW_particles.hpp
index 68f8badea0359436d3aaf3f06b02543ee9ff0d89..b577da055b77f58703cbf175c8215c6178cd7600 100644
--- a/bfps/cpp/rFFTW_particles.hpp
+++ b/bfps/cpp/rFFTW_particles.hpp
@@ -29,9 +29,9 @@
#include <iostream>
#include <hdf5.h>
#include "base.hpp"
+#include "particles_base.hpp"
#include "fluid_solver_base.hpp"
#include "rFFTW_interpolator.hpp"
-#include "particles.hpp"
#ifndef RFFTW_PARTICLES
diff --git a/bfps/fluid_base.py b/bfps/fluid_base.py
index ba44f306bce344d73ef39c51618e51b48bf98a6c..ba0be9ed8e3660dc66c9fa0a89a2eb2b8979a99f 100644
--- a/bfps/fluid_base.py
+++ b/bfps/fluid_base.py
@@ -130,9 +130,10 @@ class fluid_particle_base(bfps.code):
'hsize_t dset;\n')
for key in ['state', 'velocity', 'acceleration']:
self.definitions += ('full_name = (std::string(name) + std::string("/{0}"));\n'.format(key) +
+ 'if (H5Lexists(g_id, full_name.c_str(), H5P_DEFAULT))\n{\n' +
'dset = H5Dopen(g_id, full_name.c_str(), H5P_DEFAULT);\n' +
'grow_single_dataset(dset, niter_todo/niter_part);\n' +
- 'H5Dclose(dset);\n')
+ 'H5Dclose(dset);\n}\n')
self.definitions += ('full_name = (std::string(name) + std::string("/rhs"));\n' +
'if (H5Lexists(g_id, full_name.c_str(), H5P_DEFAULT))\n{\n' +
'dset = H5Dopen(g_id, full_name.c_str(), H5P_DEFAULT);\n' +
@@ -445,25 +446,24 @@ class fluid_particle_base(bfps.code):
dtype = np.int64,
compression = 'gzip')
for s in range(self.particle_species):
- if self.parameters['tracers{0}_integration_method'.format(s)] == 'AdamsBashforth':
- time_chunk = 2**20 // (8*3*
- self.parameters['nparticles']*
- self.parameters['tracers{0}_integration_steps'.format(s)])
- time_chunk = max(time_chunk, 1)
- ofile.create_dataset('particles/tracers{0}/rhs'.format(s),
- (1,
- self.parameters['tracers{0}_integration_steps'.format(s)],
- self.parameters['nparticles'],
- 3),
- maxshape = (None,
- self.parameters['tracers{0}_integration_steps'.format(s)],
- self.parameters['nparticles'],
- 3),
- chunks = (time_chunk,
- self.parameters['tracers{0}_integration_steps'.format(s)],
- self.parameters['nparticles'],
- 3),
- dtype = np.float64)
+ time_chunk = 2**20 // (8*3*
+ self.parameters['nparticles']*
+ self.parameters['tracers{0}_integration_steps'.format(s)])
+ time_chunk = max(time_chunk, 1)
+ ofile.create_dataset('particles/tracers{0}/rhs'.format(s),
+ (1,
+ self.parameters['tracers{0}_integration_steps'.format(s)],
+ self.parameters['nparticles'],
+ 3),
+ maxshape = (None,
+ self.parameters['tracers{0}_integration_steps'.format(s)],
+ self.parameters['nparticles'],
+ 3),
+ chunks = (time_chunk,
+ self.parameters['tracers{0}_integration_steps'.format(s)],
+ self.parameters['nparticles'],
+ 3),
+ dtype = np.float64)
time_chunk = 2**20 // (8*3*self.parameters['nparticles'])
time_chunk = max(time_chunk, 1)
ofile.create_dataset(
@@ -474,14 +474,15 @@ class fluid_particle_base(bfps.code):
chunks = (time_chunk, self.parameters['nparticles'], 3),
maxshape = (None, self.parameters['nparticles'], 3),
dtype = np.float64)
- ofile.create_dataset(
- '/particles/tracers{0}/acceleration'.format(s),
- (1,
- self.parameters['nparticles'],
- 3),
- chunks = (time_chunk, self.parameters['nparticles'], 3),
- maxshape = (None, self.parameters['nparticles'], 3),
- dtype = np.float64)
+ if self.parameters['tracers{0}_acc_on'.format(s)]:
+ ofile.create_dataset(
+ '/particles/tracers{0}/acceleration'.format(s),
+ (1,
+ self.parameters['nparticles'],
+ 3),
+ chunks = (time_chunk, self.parameters['nparticles'], 3),
+ maxshape = (None, self.parameters['nparticles'], 3),
+ dtype = np.float64)
ofile.close()
return None
diff --git a/done.txt b/done.txt
index 344832ae21046d40f47702000ab8b31d0128114b..79d3a56f93a3ca8227643655bae6e46bfa4b2ee0 100644
--- a/done.txt
+++ b/done.txt
@@ -3,3 +3,10 @@ x 2015-12-23 decide on versioning system
x 2015-12-24 move get_grid coords to interpolator @optimization +v1.0
x 2015-12-25 get rid of temporal interpolation @optimization +v1.0
x 2015-12-26 call interpolation only when needed @optimization +v1.0
+x 2015-12-26 clean up tox files, make sure all tests run @tests +v1.0
+x 2016-01-03 check divfree function
+x 2016-01-03 compute kMeta(t) as well
+x 2016-01-03 split library into core and extra @optimization +v1.0
+x 2016-01-07 FFTW interpolator doesn't need its own field @optimization +v1.0 +particle_api
+x 2016-01-08 simplify tracer/field addition mechanism @design +v1.0 +particle_api
+x 2016-01-08 add stat choice parameter to add_particles @design +v1.0 +particle_api
diff --git a/setup.py b/setup.py
index 3fdf80d57f41f6620de62d62ce25c00f6603e3b4..c60f568f5a4f0baa5b14b6548f13bc3d89c2321e 100644
--- a/setup.py
+++ b/setup.py
@@ -68,9 +68,7 @@ print(VERSION)
src_file_list = ['field_descriptor',
'fluid_solver_base',
'fluid_solver',
- 'interpolator',
'rFFTW_interpolator',
- 'particles',
'rFFTW_particles',
'fftw_tools',
'spline_n1',
@@ -81,7 +79,11 @@ src_file_list = ['field_descriptor',
'spline_n6',
'Lagrange_polys']
-header_list = ['cpp/base.hpp'] + ['cpp/' + fname + '.hpp' for fname in src_file_list]
+header_list = (['cpp/base.hpp',
+ 'cpp/particles_base.hpp',
+ 'cpp/interpolator_base.hpp'] +
+ ['cpp/' + fname + '.hpp'
+ for fname in src_file_list])
with open('MANIFEST.in', 'w') as manifest_in_file:
for fname in ['bfps/cpp/' + fname + '.cpp' for fname in src_file_list] + header_list:
diff --git a/tests/base.py b/tests/base.py
old mode 100755
new mode 100644
index c050245b85910bd35e2c62bf937945ba8d897e52..abe5cf861ecc35b9211b76bd32a7f43e84a82148
--- a/tests/base.py
+++ b/tests/base.py
@@ -102,29 +102,19 @@ def launch(
c.parameters['famplitude'] = 0.2
c.fill_up_fluid_code()
if c.parameters['nparticles'] > 0:
- c.add_particle_fields(
- name = 'regular',
+ c.add_3D_rFFTW_field(name = 'rFFTW_acc')
+ c.add_interpolator(
+ name = 'spline',
neighbours = opt.neighbours,
smoothness = opt.smoothness)
- c.add_particle_fields(kcut = 'fs->kM/2', name = 'filtered', neighbours = opt.neighbours)
c.add_particles(
- kcut = 'fs->kM/2',
- integration_steps = 1,
- fields_name = 'filtered')
- #for integr_steps in range(1, 7):
- # c.add_particles(
- # integration_steps = integr_steps,
- # neighbours = opt.neighbours,
- # smoothness = opt.smoothness,
- # fields_name = 'regular')
- for info in [(2, 'AdamsBashforth'),
- (3, 'AdamsBashforth'),
- (4, 'AdamsBashforth'),
- (6, 'AdamsBashforth')]:
- c.add_particles(
- integration_steps = info[0],
- integration_method = info[1],
- fields_name = 'regular')
+ kcut = ['fs->kM/2', 'fs->kM/3'],
+ integration_steps = 3,
+ interpolator = 'spline')
+ c.add_particles(
+ integration_steps = [2, 3, 4, 6],
+ interpolator = 'spline',
+ acc_name = 'rFFTW_acc')
c.finalize_code()
c.write_src()
c.write_par()
@@ -158,6 +148,8 @@ def launch(
def acceleration_test(c, m = 3, species = 0):
+ if not c.parameters['tracers{0}_acc_on'.format(species)]:
+ return None
import numpy as np
import matplotlib.pyplot as plt
from bfps.tools import get_fornberg_coeffs
@@ -182,13 +174,12 @@ def acceleration_test(c, m = 3, species = 0):
pid = np.argmin(SNR(num_acc1, acc[n+1:-n-1]))
pars = d['parameters']
to_print = (
- 'integration={0}, steps={1}, interp={2}, neighbours={3}, '.format(
- pars['tracers{0}_integration_method'.format(species)].value,
+ 'steps={0}, interp={1}, neighbours={2}, '.format(
pars['tracers{0}_integration_steps'.format(species)].value,
- pars[str(pars['tracers{0}_field'.format(species)].value) + '_type'].value,
- pars[str(pars['tracers{0}_field'.format(species)].value) + '_neighbours'].value))
- if 'spline' in pars['tracers{0}_field'.format(species)].value:
- to_print += 'smoothness = {0}, '.format(pars[str(pars['tracers{0}_field'.format(species)].value) + '_smoothness'].value)
+ pars[str(pars['tracers{0}_interpolator'.format(species)].value) + '_type'].value,
+ pars[str(pars['tracers{0}_interpolator'.format(species)].value) + '_neighbours'].value))
+ if 'spline' in pars['tracers{0}_interpolator'.format(species)].value:
+ to_print += 'smoothness = {0}, '.format(pars[str(pars['tracers{0}_interpolator'.format(species)].value) + '_smoothness'].value)
to_print += (
'SNR d1p-vel={0:.3f}, d1v-acc={1:.3f}, d2p-acc={2:.3f}'.format(
np.mean(SNR(num_vel1, vel[n+1:-n-1])),
diff --git a/tests/test_convergence.py b/tests/test_convergence.py
old mode 100755
new mode 100644
index 636dea9bc9487241b62166d35ac7367272e513b3..f00e8df7e46388c461b59c8f5b34c62936344baf
--- a/tests/test_convergence.py
+++ b/tests/test_convergence.py
@@ -1,4 +1,3 @@
-#! /usr/bin/env python2
#######################################################################
# #
# Copyright 2015 Max Planck Institute #
@@ -30,7 +29,7 @@ import h5py
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
-from test_base import *
+from base import *
def convergence_test(
opt,
@@ -224,6 +223,15 @@ def convergence_test(
return c0, c1, c2
if __name__ == '__main__':
- opt = parser.parse_args()
+ opt = parser.parse_args(
+ ['-n', '16',
+ '--run',
+ '--initialize',
+ '--ncpu', '2',
+ '--nparticles', '1000',
+ '--niter_todo', '32',
+ '--precision', 'single',
+ '--wd', 'data/single'] +
+ sys.argv[1:])
convergence_test(opt, launch)
diff --git a/tests/test_frozen_field.py b/tests/test_frozen_field.py
old mode 100755
new mode 100644
index ca42527c977f16477faeca651202f23cd3f4589c..4bf10325a348f59533ac309d8407a1eecbbeb989
--- a/tests/test_frozen_field.py
+++ b/tests/test_frozen_field.py
@@ -1,4 +1,3 @@
-#! /usr/bin/env python2
#######################################################################
# #
# Copyright 2015 Max Planck Institute #
@@ -25,7 +24,7 @@
-from test_base import *
+from base import *
class FrozenFieldParticles(bfps.NavierStokes):
def __init__(
@@ -33,12 +32,16 @@ class FrozenFieldParticles(bfps.NavierStokes):
name = 'FrozenFieldParticles',
work_dir = './',
simname = 'test',
- fluid_precision = 'single'):
+ frozen_fields = True,
+ fluid_precision = 'single',
+ use_fftw_wisdom = False):
super(FrozenFieldParticles, self).__init__(
name = name,
work_dir = work_dir,
simname = simname,
- fluid_precision = fluid_precision)
+ fluid_precision = fluid_precision,
+ frozen_fields = frozen_fields,
+ use_fftw_wisdom = use_fftw_wisdom)
def fill_up_fluid_code(self):
self.fluid_includes += '#include <cstring>\n'
self.fluid_variables += 'fluid_solver<{0}> *fs;\n'.format(self.C_dtype)
diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py
old mode 100755
new mode 100644
index 928ae115370019c6b32a8720efddb77376d86a09..29755bd6a5207c322c7e442df89be28b7a996543
--- a/tests/test_interpolation.py
+++ b/tests/test_interpolation.py
@@ -1,4 +1,3 @@
-#! /usr/bin/env python2
#######################################################################
# #
# Copyright 2015 Max Planck Institute #
@@ -26,35 +25,44 @@
import numpy as np
-from test_base import *
+from base import *
import matplotlib.pyplot as plt
-class test_interpolation(bfps.NavierStokes):
- def __init__(
- self,
- name = 'test_interpolation',
- work_dir = './',
- simname = 'test'):
- super(test_interpolation, self).__init__(
- work_dir = work_dir,
- simname = simname,
- name = name,
- frozen_fields = True)
- return None
-
if __name__ == '__main__':
- opt = parser.parse_args()
- c = test_interpolation(work_dir = opt.work_dir + '/io')
+ opt = parser.parse_args(
+ ['-n', '32',
+ '--run',
+ '--ncpu', '2',
+ '--initialize',
+ '--neighbours', '4',
+ '--smoothness', '1',
+ '--nparticles', '4225',
+ '--niter_todo', '1',
+ '--niter_stat', '1',
+ '--niter_part', '1',
+ '--niter_out', '1',
+ '--wd', 'interp'] +
+ sys.argv[1:])
+ c = bfps.NavierStokes(
+ name = 'test_interpolation',
+ use_fftw_wisdom = False,
+ frozen_fields = True)
c.pars_from_namespace(opt)
+ c.fill_up_fluid_code()
+ c.add_particle_fields(
+ name = 'n{0}m{1}'.format(opt.neighbours, opt.smoothness),
+ neighbours = opt.neighbours,
+ smoothness = opt.smoothness)
+ c.add_particle_fields(
+ name = 'n{0}'.format(opt.neighbours),
+ neighbours = opt.neighbours,
+ interp_type = 'Lagrange')
c.add_particles(
integration_steps = 1,
- neighbours = 4,
- smoothness = 1)
+ fields_name = 'n{0}m{1}'.format(opt.neighbours, opt.smoothness))
c.add_particles(
integration_steps = 1,
- neighbours = 4,
- interp_type = 'Lagrange')
- c.fill_up_fluid_code()
+ fields_name = 'n{0}'.format(opt.neighbours))
c.finalize_code()
c.write_src()
c.write_par()
diff --git a/tests/test_io.py b/tests/test_io.py
old mode 100755
new mode 100644
index 0d8b8c055eae1deea8fbaef0644f050d8a900e3c..447cb405406a61573cc7d7b2b5942d012e6f55f6
--- a/tests/test_io.py
+++ b/tests/test_io.py
@@ -1,4 +1,3 @@
-#! /usr/bin/env python2
#######################################################################
# #
# Copyright 2015 Max Planck Institute #
@@ -25,7 +24,7 @@
-from test_base import *
+from base import *
class test_io(bfps.code):
def __init__(
@@ -43,7 +42,12 @@ class test_io(bfps.code):
return None
if __name__ == '__main__':
- opt = parser.parse_args()
+ opt = parser.parse_args(
+ ['-n', '32',
+ '--run',
+ '--initialize',
+ '--ncpu', '2'] +
+ sys.argv[1:])
c = test_io(work_dir = opt.work_dir + '/io')
c.write_src()
c.write_par()
diff --git a/tests/test_particles.py b/tests/test_particles.py
old mode 100755
new mode 100644
index 93801df1c7263964c8bf1ae2a189a826b4cc2763..92f82b901624ad0b48cf9ba5e3661de7f5071f03
--- a/tests/test_particles.py
+++ b/tests/test_particles.py
@@ -28,7 +28,7 @@
import numpy as np
import matplotlib.pyplot as plt
-from test_base import *
+from base import *
from test_frozen_field import FrozenFieldParticles
from test_convergence import convergence_test
@@ -147,7 +147,17 @@ class err_finder:
return errlist
if __name__ == '__main__':
- opt = parser.parse_args()
+ opt = parser.parse_args(
+ ['-n', '16',
+ '--run',
+ '--initialize',
+ '--frozen',
+ '--ncpu', '2',
+ '--nparticles', '128',
+ '--niter_todo', '16',
+ '--precision', 'single',
+ '--wd', 'data/single'] +
+ sys.argv[1:])
if opt.precision == 'single':
dtype = np.complex64
elif opt.precision == 'double':
diff --git a/tests/test_plain.py b/tests/test_plain.py
old mode 100755
new mode 100644
index 8ada988ad741d9a2fe177c8e2ea09369a9e69d39..5d526cfd3412b097068f9ab1b93b2ca86523b34d
--- a/tests/test_plain.py
+++ b/tests/test_plain.py
@@ -65,6 +65,7 @@ if __name__ == '__main__':
opt = parser.parse_args(
['-n', '32',
'--run',
+ '--initialize',
'--ncpu', '2',
'--nparticles', '1000',
'--niter_todo', '48',
diff --git a/tests/test_scaling.py b/tests/test_scaling.py
old mode 100755
new mode 100644
diff --git a/tests/test_time_step.py b/tests/test_time_step.py
old mode 100755
new mode 100644
index fd1662c8d812f59dd9faa9e4bcf9f60ef289350b..9a24619b94bc06d77b712018544aa663e32665c0
--- a/tests/test_time_step.py
+++ b/tests/test_time_step.py
@@ -1,4 +1,3 @@
-#! /usr/bin/env python2
#######################################################################
# #
# Copyright 2015 Max Planck Institute #
@@ -30,7 +29,7 @@ import h5py
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
-from test_base import *
+from base import *
def convergence_test(
opt,
@@ -93,10 +92,19 @@ def convergence_test(
a.set_xscale('log')
a.set_yscale('log')
a.set_xlabel('$\\|u\\|_\\infty \\frac{\\Delta t}{\\Delta x}$')
- fig.savefig('spec_err_vs_dt_{0}.pdf'.format(opt.precision))
+ fig.savefig('vel_err_vs_dt_{0}.pdf'.format(opt.precision))
return None
if __name__ == '__main__':
- opt = parser.parse_args()
+ opt = parser.parse_args(
+ ['-n', '32',
+ '--run',
+ '--initialize',
+ '--ncpu', '2',
+ '--nparticles', '1000',
+ '--niter_todo', '16',
+ '--precision', 'single',
+ '--wd', 'data/single'] +
+ sys.argv[1:])
convergence_test(opt, launch)
diff --git a/todo.txt b/todo.txt
index 2356510f3170e2841465fa2231af6b11411fcd5b..41ab0f009a43dd95179363bdb499aab43624983c 100644
--- a/todo.txt
+++ b/todo.txt
@@ -1,16 +1,18 @@
-(B) FFTW interpolator doesn't need its own field @optimization +v1.0
-(B) clean up tox files, make sure all tests run @tests +v1.0
-(B) compute z polynomials only when needed @optimization +v1.0
-(B) use less memory @optimization
-(B) split library into core and extra @optimization +v1.0
-(C) clean up machine_settings mess @design @documentation +v2.0
-(C) code overview @documentation
-(C) move stat I/O to cpp lib @design @HDF5
-(C) set up mechanism for adding in new PDEs @design +v2.0 +alternate_algorithms
-(C) test involving hydrodynamic similarity @tests
-(C) use HDF5 io for fields @design @HDF5 +I/O
-(D) test anisotropic grids @tests
-(D) test non-cubic domains @tests
-(E) add u-equation algorithm for testing purposes @tests +alternate_algorithms
-(E) pure python DNS addon: pros and cons @tests +alternate_algorithms
+(B) compute z polynomials only when needed @optimization +v1.0
+(B) read https://www.xsede.org/documents/271087/369161/ExtScale-Koziol.pdf @optimization @HDF5 +I/O
+(B) set up mechanism for adding in new PDEs @design +v2.0 +alternate_algorithms
+(B) tweak HDF5 settings @optimization @HDF5 +I/O
+(B) use less memory @optimization
+(C) clean up machine_settings mess @design @documentation +v2.0
+(C) code overview @documentation
+(C) move stat I/O to cpp lib @design @HDF5
+(C) test involving hydrodynamic similarity @tests
+(C) use HDF5 io for fields @design @HDF5 +I/O
+(D) generalize interpolation comparison test @tests
+(D) generate separate lib(s) with extra classes @tests +alternate_algorithms
+(D) test anisotropic grids @tests
+(D) test non-cubic domains @tests
+(D) tests should not overwrite other tests (tox_full) @tests
+(E) add u-equation algorithm for testing purposes @tests +alternate_algorithms
+(E) pure python DNS addon: pros and cons @tests +alternate_algorithms
(F) add switch to turn off simulation
diff --git a/tox_convergence.ini b/tox_convergence.ini
deleted file mode 100644
index 39a368ff0d1954cb841cbaab41e5a8d98bb68031..0000000000000000000000000000000000000000
--- a/tox_convergence.ini
+++ /dev/null
@@ -1,31 +0,0 @@
-[tox]
-envlist = py34
-[testenv]
-deps = matplotlib
-whitelist_externals =
- echo
- cp
-passenv = LD_LIBRARY_PATH
-changedir =
- {envtmpdir}
-commands =
- cp -r {toxinidir}/tests {envtmpdir}
- python tests/test_convergence.py \
- -n 32 \
- --run \
- --initialize \
- --ncpu 2 \
- --nparticles 1 \
- --niter_todo 16 \
- --precision single \
- --wd "data/single"
- python tests/test_convergence.py \
- -n 32 \
- --run \
- --initialize \
- --ncpu 2 \
- --nparticles 1 \
- --niter_todo 16 \
- --precision double \
- --wd "data/double"
-
diff --git a/tox_full.ini b/tox_full.ini
new file mode 100644
index 0000000000000000000000000000000000000000..d102b91ef0cb7b305d93d08d3dd899f80178a2f4
--- /dev/null
+++ b/tox_full.ini
@@ -0,0 +1,25 @@
+[tox]
+envlist = py34
+[testenv]
+sitepackages = True
+whitelist_externals =
+ echo
+ cp
+passenv = LD_LIBRARY_PATH
+changedir =
+ {envtmpdir}
+commands =
+ cp -r {toxinidir}/tests {envtmpdir}
+ python tests/test_io.py
+ python tests/test_interpolation.py
+ python tests/test_interpolation.py --precision double --wd "interp_double"
+ python tests/test_plain.py
+ python tests/test_plain.py --precision double --wd "data/double"
+ python tests/test_scaling.py
+ python tests/test_scaling.py --precision double --wd "data/double"
+ python tests/test_particles.py
+ python tests/test_particles.py --precision double --wd "data/double"
+ python tests/test_time_step.py
+ python tests/test_time_step.py --precision double --wd "data/double"
+ python tests/test_convergence.py
+ python tests/test_convergence.py --precision double --wd "data/double"
diff --git a/tox_interpolation.ini b/tox_interpolation.ini
deleted file mode 100644
index 71b94b2220f43496050b09855cf7e41b9e2da51e..0000000000000000000000000000000000000000
--- a/tox_interpolation.ini
+++ /dev/null
@@ -1,21 +0,0 @@
-[tox]
-envlist = py27
-[testenv]
-whitelist_externals =
- echo
- cp
-passenv = LD_LIBRARY_PATH
-changedir =
- {envtmpdir}
-commands =
- cp -r {toxinidir}/tests {envtmpdir}
- python tests/test_interpolation.py \
- -n 32 --ncpu 2 \
- --initialize \
- --nparticles 4225 \
- --niter_todo 1 \
- --niter_stat 1 \
- --niter_part 1 \
- --niter_out 1 \
- --run
-
diff --git a/tox_io.ini b/tox_io.ini
deleted file mode 100644
index de8c2fc5d0bb5fa7498527ecb4b96b32ca000c69..0000000000000000000000000000000000000000
--- a/tox_io.ini
+++ /dev/null
@@ -1,14 +0,0 @@
-[tox]
-envlist = py27, py34
-[testenv]
-deps = matplotlib
-whitelist_externals =
- echo
- cp
-passenv = LD_LIBRARY_PATH
-changedir =
- {envtmpdir}
-commands =
- cp -r {toxinidir}/tests {envtmpdir}
- python tests/test_io.py -n 32 --run --initialize --ncpu 2
-
diff --git a/tox_particles.ini b/tox_particles.ini
deleted file mode 100644
index e76f1b387753e12f9e6b6270f2959c1e1bca4ba5..0000000000000000000000000000000000000000
--- a/tox_particles.ini
+++ /dev/null
@@ -1,30 +0,0 @@
-[tox]
-envlist = py27
-[testenv]
-whitelist_externals =
- echo
- cp
-passenv = LD_LIBRARY_PATH
-changedir =
- {envtmpdir}
-commands =
- cp -r {toxinidir}/tests {envtmpdir}
- python tests/test_particles.py \
- -n 16 \
- --run \
- --initialize \
- --ncpu 2 \
- --nparticles 32 \
- --niter_todo 16 \
- --precision single \
- --wd "data/single"
- python tests/test_particles.py \
- -n 16 \
- --run \
- --initialize \
- --ncpu 2 \
- --nparticles 32 \
- --niter_todo 16 \
- --precision double \
- --wd "data/double"
-
diff --git a/tox_plain.ini b/tox_plain.ini
deleted file mode 100644
index 38f34b15fee4d479f34cbf805adebb61be8ec218..0000000000000000000000000000000000000000
--- a/tox_plain.ini
+++ /dev/null
@@ -1,30 +0,0 @@
-[tox]
-envlist = py34
-[testenv]
-sitepackages = True
-whitelist_externals =
- echo
- cp
- g++
-passenv = LD_LIBRARY_PATH
-changedir =
- {envtmpdir}
-commands =
- cp -r {toxinidir}/tests {envtmpdir}
- #python tests/test_plain.py -n 256 --run --initialize --ncpu 8 --niter_todo 8 --precision single --wd "data/single"
- #python tests/test_plain.py \
- # -n 128 --run --initialize --ncpu 2 \
- # --nparticles 100 --niter_todo 24 \
- # --precision single --wd "data/single" #\
- # #--multiplejob
- python tests/test_plain.py -n 32 \
- --run --initialize --multiplejob \
- --ncpu 2 --nparticles 16 \
- --niter_todo 48 \
- --precision single --wd "data/single"
- python tests/test_plain.py -n 32 \
- --run --initialize --multiplejob \
- --ncpu 2 --nparticles 16 \
- --niter_todo 48 \
- --precision double --wd "data/double"
-
diff --git a/tox_time_step.ini b/tox_time_step.ini
deleted file mode 100644
index 87ec18d1af8b8290d18a3cd1e8c3e9b5d671bb6c..0000000000000000000000000000000000000000
--- a/tox_time_step.ini
+++ /dev/null
@@ -1,30 +0,0 @@
-[tox]
-envlist = py27
-[testenv]
-whitelist_externals =
- echo
- cp
-passenv = LD_LIBRARY_PATH
-changedir =
- {envtmpdir}
-commands =
- cp -r {toxinidir}/tests {envtmpdir}
- python tests/test_time_step.py \
- -n 32 \
- --run \
- --initialize \
- --ncpu 2 \
- --nparticles 32 \
- --niter_todo 16 \
- --precision single \
- --wd "data/single"
- python tests/test_time_step.py \
- -n 32 \
- --run \
- --initialize \
- --ncpu 2 \
- --nparticles 32 \
- --niter_todo 16 \
- --precision double \
- --wd "data/double"
-