diff --git a/bfps/DNS.py b/bfps/DNS.py index 90988019734e777e8b355eff4cd33763e1f71857..7eb15af64dbc5eb1b3f76b6f8309030cf32d2d69 100644 --- a/bfps/DNS.py +++ b/bfps/DNS.py @@ -127,7 +127,7 @@ class DNS(_code): template_class = '{0}<{1}>::'.format(self.dns_type, rnumber), template_prefix = 'template '.format(rnumber), just_declaration = True) + '\n\n') - if self.dns_type in ['NSVEparticles', 'NSVE_no_output', 'NSVEparticles_no_output', 'NSVEparticlesP2P']: + if self.dns_type in ['NSVEparticles', 'NSVE_no_output', 'NSVEparticles_no_output', 'NSVEcomplex_particles']: outfile.write('template <typename rnumber> int NSVE<rnumber>::read_parameters(){return EXIT_SUCCESS;}\n') outfile.write('template int NSVE<float>::read_parameters();\n') outfile.write('template int NSVE<double>::read_parameters();\n\n') @@ -167,9 +167,9 @@ class DNS(_code): self.NSVEp_extra_parameters['tracers0_neighbours'] = int(1) self.NSVEp_extra_parameters['tracers0_smoothness'] = int(1) #self.extra_parameters = {} - #for key in ['NSVE', 'NSVE_no_output', 'NSVEparticles', 'NSVEparticles_no_output', 'NSVEparticlesP2P']: + #for key in ['NSVE', 'NSVE_no_output', 'NSVEparticles', 'NSVEparticles_no_output', 'NSVEcomplex_particles']: # self.extra_parameters[key] = {} - #for key in ['NSVEparticles', 'NSVEparticles_no_output', 'NSVEparticlesP2P']: + #for key in ['NSVEparticles', 'NSVEparticles_no_output', 'NSVEcomplex_particles']: # self.extra_parameters[key].update(self.NSVEp_extra_parameters) return None def get_kspace(self): @@ -387,7 +387,7 @@ class DNS(_code): assert (self.parameters['niter_todo'] % self.parameters['niter_stat'] == 0) assert (self.parameters['niter_todo'] % self.parameters['niter_out'] == 0) assert (self.parameters['niter_out'] % self.parameters['niter_stat'] == 0) - if self.dns_type in ['NSVEparticles_no_output', 'NSVEparticlesP2P', 'NSVEparticles']: + if self.dns_type in ['NSVEparticles_no_output', 'NSVEcomplex_particles', 'NSVEparticles']: assert (self.parameters['niter_todo'] % self.parameters['niter_part'] == 0) assert (self.parameters['niter_out'] % self.parameters['niter_part'] == 0) _code.write_par(self, iter0 = iter0) @@ -443,7 +443,7 @@ class DNS(_code): for val in pbase_shape[1:]: number_of_particles *= val ncomponents = 3 - if self.dns_type in ['NSVEparticlesP2P']: + if self.dns_type in ['NSVEcomplex_particles']: ncomponents = 6 with h5py.File(self.get_checkpoint_0_fname(), 'a') as ofile: s = 0 @@ -622,8 +622,8 @@ class DNS(_code): help = 'plain Navier-Stokes vorticity formulation, with basic fluid tracers') parser_NSVEp2p = subparsers.add_parser( - 'NSVEparticlesP2P', - help = 'plain Navier-Stokes vorticity formulation, with basic fluid tracers') + 'NSVEcomplex_particles', + help = 'plain Navier-Stokes vorticity formulation, with oriented active particles') for parser in ['NSVEparticles_no_output', 'NSVEp2', 'NSVEp2p']: eval('self.simulation_parser_arguments({0})'.format('parser_' + parser)) @@ -668,7 +668,7 @@ class DNS(_code): self.dns_type = opt.DNS_class self.name = self.dns_type + '-' + self.fluid_precision + '-v' + bfps.__version__ # merge parameters if needed - if self.dns_type in ['NSVEparticles', 'NSVEparticlesP2P', 'NSVEparticles_no_output']: + if self.dns_type in ['NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output']: for k in self.NSVEp_extra_parameters.keys(): self.parameters[k] = self.NSVEp_extra_parameters[k] if type(extra_parameters) != type(None): @@ -728,7 +728,7 @@ class DNS(_code): # hardcoded FFTW complex representation size field_size = 3*(opt.nx+2)*opt.ny*opt.nz*self.fluid_dtype.itemsize checkpoint_size = field_size - if self.dns_type in ['NSVEparticles', 'NSVEparticlesP2P', 'NSVEparticles_no_output']: + if self.dns_type in ['NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output']: rhs_size = self.parameters['tracers0_integration_steps'] if type(opt.tracers0_integration_steps) != type(None): rhs_size = opt.tracers0_integration_steps @@ -932,9 +932,9 @@ class DNS(_code): particle_file.create_group('tracers0/position') particle_file.create_group('tracers0/velocity') particle_file.create_group('tracers0/acceleration') - if self.dns_type in ['NSVEparticlesP2P']: + if self.dns_type in ['NSVEcomplex_particles']: particle_file.create_group('tracers0/orientation') - particle_file.create_group('tracers0/vorticity') + particle_file.create_group('tracers0/velocity_gradient') return None def launch_jobs( self, @@ -994,7 +994,7 @@ class DNS(_code): # particle_initial_condition[..., 2] += onedarray[None, :, None, None] self.write_par( particle_ic = None) - if self.dns_type in ['NSVEparticles', 'NSVEparticlesP2P', 'NSVEparticles_no_output']: + if self.dns_type in ['NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output']: self.generate_particle_data(opt = opt) self.run( nb_processes = opt.nb_processes, diff --git a/bfps/cpp/full_code/NSVEcomplex_particles.cpp b/bfps/cpp/full_code/NSVEcomplex_particles.cpp index a6b7082e6b2ffa16e456f32d439e8b4cfecf451d..93d0edc17e24925b25513dd3a4832f687255329a 100644 --- a/bfps/cpp/full_code/NSVEcomplex_particles.cpp +++ b/bfps/cpp/full_code/NSVEcomplex_particles.cpp @@ -1,13 +1,13 @@ #include <string> #include <cmath> -#include "NSVEparticlesP2P.hpp" +#include "NSVEcomplex_particles.hpp" #include "scope_timer.hpp" #include "particles/particles_sampling.hpp" #include "particles/p2p_computer.hpp" #include "particles/particles_inner_computer.hpp" template <typename rnumber> -int NSVEparticlesP2P<rnumber>::initialize(void) +int NSVEcomplex_particles<rnumber>::initialize(void) { this->NSVE<rnumber>::initialize(); @@ -50,14 +50,18 @@ int NSVEparticlesP2P<rnumber>::initialize(void) (this->simname + "_particles.h5"), "tracers0", "position/0"); - // TODO: remove the following testing initial condition, and use a proper - // way to initialize with 0 (i.e. generate a 0 field as the initial condition). - //*this->fs->cvorticity = 0.0; + + + /// allocate grad vel field + this->nabla_u = new field<rnumber, FFTW, THREExTHREE>( + this->nx, this->ny, this->nz, + this->comm, + this->fs->cvorticity->fftw_plan_rigor); return EXIT_SUCCESS; } template <typename rnumber> -int NSVEparticlesP2P<rnumber>::step(void) +int NSVEcomplex_particles<rnumber>::step(void) { this->fs->compute_velocity(this->fs->cvorticity); this->fs->cvelocity->ift(); @@ -74,7 +78,7 @@ int NSVEparticlesP2P<rnumber>::step(void) } template <typename rnumber> -int NSVEparticlesP2P<rnumber>::write_checkpoint(void) +int NSVEcomplex_particles<rnumber>::write_checkpoint(void) { this->NSVE<rnumber>::write_checkpoint(); this->particles_output_writer_mpi->open_file(this->fs->get_current_fname()); @@ -90,8 +94,9 @@ int NSVEparticlesP2P<rnumber>::write_checkpoint(void) } template <typename rnumber> -int NSVEparticlesP2P<rnumber>::finalize(void) +int NSVEcomplex_particles<rnumber>::finalize(void) { + delete this->nabla_u; delete this->particles_output_writer_mpi; delete this->particles_sample_writer_mpi; this->NSVE<rnumber>::finalize(); @@ -102,7 +107,7 @@ int NSVEparticlesP2P<rnumber>::finalize(void) */ template <typename rnumber> -int NSVEparticlesP2P<rnumber>::do_stats() +int NSVEcomplex_particles<rnumber>::do_stats() { /// perform fluid stats this->NSVE<rnumber>::do_stats(); @@ -115,6 +120,7 @@ int NSVEparticlesP2P<rnumber>::do_stats() /// allocate temporary data array std::unique_ptr<double[]> pdata0 = this->ps->extractParticlesState(0, 3); std::unique_ptr<double[]> pdata1 = this->ps->extractParticlesState(3, 6); + std::unique_ptr<double[]> pdata2(new double[9*this->ps->getLocalNbParticles()]); /// sample position this->particles_sample_writer_mpi->template save_dataset<3>( @@ -147,15 +153,20 @@ int NSVEparticlesP2P<rnumber>::do_stats() this->ps->getLocalNbParticles(), this->ps->get_step_idx()-1); - /// sample vorticity - *this->tmp_vec_field = this->fs->cvorticity->get_cdata(); - this->tmp_vec_field->ift(); - this->ps->sample_compute_field(*this->tmp_vec_field, pdata1.get()); - this->particles_sample_writer_mpi->template save_dataset<3>( + /// sample velocity gradient + /// fs->cvelocity should contain the velocity in Fourier space + this->fs->compute_velocity(this->fs->cvorticity); + compute_gradient( + this->fs->kk, + this->fs->cvelocity, + this->nabla_u); + this->nabla_u->ift(); + this->ps->sample_compute_field(*this->nabla_u, pdata2.get()); + this->particles_sample_writer_mpi->template save_dataset<9>( "tracers0", - "vorticity", + "velocity_gradient", pdata0.get(), - &pdata1, + &pdata2, this->ps->getParticlesIndexes(), this->ps->getLocalNbParticles(), this->ps->get_step_idx()-1); @@ -176,6 +187,6 @@ int NSVEparticlesP2P<rnumber>::do_stats() return EXIT_SUCCESS; } -template class NSVEparticlesP2P<float>; -template class NSVEparticlesP2P<double>; +template class NSVEcomplex_particles<float>; +template class NSVEcomplex_particles<double>; diff --git a/bfps/cpp/full_code/NSVEcomplex_particles.hpp b/bfps/cpp/full_code/NSVEcomplex_particles.hpp index b74169f5b35dabf66433cbab953e290fb6ceb943..2015ec5bbec4c9e978c1ab658e4c91b97c85ba9e 100644 --- a/bfps/cpp/full_code/NSVEcomplex_particles.hpp +++ b/bfps/cpp/full_code/NSVEcomplex_particles.hpp @@ -24,8 +24,8 @@ -#ifndef NSVEPARTICLESP2P_HPP -#define NSVEPARTICLESP2P_HPP +#ifndef NSVECOMPLEX_PARTICLES_HPP +#define NSVECOMPLEX_PARTICLES_HPP @@ -37,15 +37,18 @@ #include "particles/particles_output_hdf5.hpp" #include "particles/particles_sampling.hpp" -/** \brief Navier-Stokes solver that includes simple Lagrangian tracers. +/** \brief Navier-Stokes solver that includes complex particles. * * Child of Navier Stokes vorticity equation solver, this class calls all the - * methods from `NSVE`, and in addition integrates simple Lagrangian tracers + * methods from `NSVE`, and in addition integrates `complex particles` * in the resulting velocity field. + * By `complex particles` we mean neutrally buoyant, very small particles, + * which have an orientation and actively swim in that direction, and they may + * also interact with each other, trying to reorient to a common orientation. */ template <typename rnumber> -class NSVEparticlesP2P: public NSVE<rnumber> +class NSVEcomplex_particles: public NSVE<rnumber> { public: @@ -67,16 +70,18 @@ class NSVEparticlesP2P: public NSVE<rnumber> // TODO P2P use a reader with particle data particles_output_hdf5<long long int, double,6> *particles_output_writer_mpi; particles_output_sampling_hdf5<long long int, double, 3> *particles_sample_writer_mpi; + // field for sampling velocity gradient + field<rnumber, FFTW, THREExTHREE> *nabla_u; - NSVEparticlesP2P( + NSVEcomplex_particles( const MPI_Comm COMMUNICATOR, const std::string &simulation_name): NSVE<rnumber>( COMMUNICATOR, simulation_name), cutoff(10), inner_v0(1), enable_p2p(true), enable_inner(true), enable_vorticity_omega(true){} - ~NSVEparticlesP2P(){} + ~NSVEcomplex_particles(){} int initialize(void); int step(void); @@ -87,5 +92,5 @@ class NSVEparticlesP2P: public NSVE<rnumber> int do_stats(void); }; -#endif//NSVEPARTICLESP2P_HPP +#endif//NSVECOMPLEX_PARTICLES_HPP diff --git a/bfps/test/test_particles.py b/bfps/test/test_particles.py index 40c428bd25ca5d315dde53c4d68f9ac5a79df738..6c12fac17383967619a9054c55d85dd8db7f3d4b 100644 --- a/bfps/test/test_particles.py +++ b/bfps/test/test_particles.py @@ -20,7 +20,7 @@ def main(): if sys.argv[2] == 'on': c = DNS() c.launch( - ['NSVEparticlesP2P', + ['NSVEcomplex_particles', '-n', '32', '--src-simname', 'B32p1e4', '--src-wd', bfps.lib_dir + '/test', diff --git a/setup.py b/setup.py index a160ce4ba86bf85f509ccdabbacb13c855c027de..ff73b945fac67f4db1f67e224a56c0a1fed27a79 100644 --- a/setup.py +++ b/setup.py @@ -88,7 +88,7 @@ print('This is bfps version ' + VERSION) ### lists of files and MANIFEST.in -src_file_list = ['full_code/NSVEparticlesP2P', +src_file_list = ['full_code/NSVEcomplex_particles', 'full_code/joint_acc_vel_stats', 'full_code/test', 'full_code/filter_test',