diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py index 122b9b3a8cb04c4b6d5da16e54a62117bb8fcab2..9c88e40111c1139026fe4849a88774e8ace7d885 100644 --- a/bfps/NavierStokes.py +++ b/bfps/NavierStokes.py @@ -388,7 +388,6 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base): def add_interpolator( self, interp_type = 'spline', - kcut = None, neighbours = 1, smoothness = 1, name = 'field_interpolator', @@ -418,12 +417,17 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base): interpolator = 'field_interpolator', frozen_particles = False, acc_name = None): - """ - :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 + """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' @@ -452,6 +456,8 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base): 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}_interpolator'.format(s0 + s)] = interpolator[s] self.parameters['tracers{0}_integration_steps'.format(s0 + s)] = integration_steps[s] self.file_datasets_grow += """ @@ -470,14 +476,18 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base): # 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' + - 'fs->ift_velocity();\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); @@ -528,10 +538,13 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base): output_vel_acc += '}\n' #### initialize, stepping and finalize code - update_fields = ('fs->compute_velocity(fs->cvorticity);\n' + - 'fs->ift_velocity();\n') - self.fluid_start += update_fields - self.fluid_loop += update_fields + 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' + @@ -556,6 +569,10 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base): '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) diff --git a/tests/base.py b/tests/base.py index dc2f5fe312a233f4237262dcf81664c774c0dfe0..79c2646a83cede058ece671aaa300a1c0526f4e3 100644 --- a/tests/base.py +++ b/tests/base.py @@ -107,15 +107,14 @@ def launch( name = 'spline', neighbours = opt.neighbours, smoothness = opt.smoothness) + c.add_particles( + 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.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') c.finalize_code() c.write_src() c.write_par()