/****************************************************************************** * * * Copyright 2019 Max Planck Institute for Dynamics and Self-Organization * * * * This file is part of TurTLE. * * * * TurTLE 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. * * * * TurTLE 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 TurTLE. If not, see <http://www.gnu.org/licenses/> * * * * Contact: Cristian.Lalescu@ds.mpg.de * * * ******************************************************************************/ #ifndef PARTICLES_SYSTEM_HPP #define PARTICLES_SYSTEM_HPP #include <array> #include <set> #include <algorithm> #include <cmath> #include "particles/abstract_particles_system.hpp" #include "particles/abstract_particles_system_with_p2p.hpp" #include "particles/particles_distr_mpi.hpp" #include "particles/particles_output_hdf5.hpp" #include "particles/particles_output_mpiio.hpp" #include "particles/interpolation/particles_field_computer.hpp" #include "particles/abstract_particles_input.hpp" #include "particles/particles_adams_bashforth.hpp" #include "scope_timer.hpp" #include "particles/p2p/p2p_distr_mpi.hpp" template <class partsize_t, class real_number, class field_rnumber, class field_class, class interpolator_class, int interp_neighbours, int size_particle_positions, int size_particle_rhs, class p2p_computer_class, class particles_inner_computer_class> class particles_system : public abstract_particles_system_with_p2p<partsize_t, real_number, p2p_computer_class> { static_assert(size_particle_positions >= 3, "There should be at least the positions X,Y,Z in the state"); MPI_Comm mpi_com; const std::pair<int,int> current_partition_interval; const int partition_interval_size; interpolator_class interpolator; particles_distr_mpi<partsize_t, real_number> particles_distr; particles_adams_bashforth<partsize_t, real_number, size_particle_positions, size_particle_rhs> positions_updater; using computer_class = particles_field_computer<partsize_t, real_number, interpolator_class, interp_neighbours>; computer_class computer; const field_class& default_field; std::unique_ptr<partsize_t[]> current_my_nb_particles_per_partition; std::unique_ptr<partsize_t[]> current_offset_particles_for_partition; const std::array<real_number,3> spatial_box_width; const std::array<real_number,3> spatial_partition_width; const real_number my_spatial_low_limit; const real_number my_spatial_up_limit; std::unique_ptr<real_number[]> my_particles_positions; std::unique_ptr<partsize_t[]> my_particles_positions_indexes; partsize_t my_nb_particles; partsize_t total_nb_particles; std::vector<std::unique_ptr<real_number[]>> my_particles_rhs; std::vector<hsize_t> particle_file_layout; int step_idx; p2p_distr_mpi<partsize_t, real_number> distr_p2p; p2p_computer_class computer_p2p; particles_inner_computer_class computer_particules_inner; public: particles_system(const std::array<size_t,3>& field_grid_dim, const std::array<real_number,3>& in_spatial_box_width, const std::array<real_number,3>& in_spatial_box_offset, const std::array<real_number,3>& in_spatial_partition_width, const real_number in_my_spatial_low_limit, const real_number in_my_spatial_up_limit, const std::array<size_t,3>& in_local_field_dims, const std::array<size_t,3>& in_local_field_offset, const field_class& in_field, MPI_Comm in_mpi_com, const partsize_t in_total_nb_particles, const real_number in_cutoff, p2p_computer_class in_computer_p2p, particles_inner_computer_class in_computer_particules_inner, const int in_current_iteration = 1) : mpi_com(in_mpi_com), current_partition_interval({in_local_field_offset[IDXC_Z], in_local_field_offset[IDXC_Z] + in_local_field_dims[IDXC_Z]}), partition_interval_size(current_partition_interval.second - current_partition_interval.first), interpolator(), particles_distr(in_mpi_com, current_partition_interval, field_grid_dim), positions_updater(), computer(field_grid_dim, current_partition_interval, interpolator, in_spatial_box_width, in_spatial_box_offset, in_spatial_partition_width), default_field(in_field), spatial_box_width(in_spatial_box_width), spatial_partition_width(in_spatial_partition_width), my_spatial_low_limit(in_my_spatial_low_limit), my_spatial_up_limit(in_my_spatial_up_limit), my_nb_particles(0), total_nb_particles(in_total_nb_particles), step_idx(in_current_iteration), distr_p2p(in_mpi_com, current_partition_interval, field_grid_dim, spatial_box_width, in_spatial_box_offset, in_cutoff), computer_p2p(std::move(in_computer_p2p)), computer_particules_inner(std::move(in_computer_particules_inner)){ current_my_nb_particles_per_partition.reset(new partsize_t[partition_interval_size]); current_offset_particles_for_partition.reset(new partsize_t[partition_interval_size+1]); } ~particles_system(){ } void init(abstract_particles_input<partsize_t, real_number>& particles_input) { TIMEZONE("particles_system::init"); my_particles_positions = particles_input.getMyParticles(); my_particles_positions_indexes = particles_input.getMyParticlesIndexes(); my_particles_rhs = particles_input.getMyRhs(); my_nb_particles = particles_input.getLocalNbParticles(); for(partsize_t idx_part = 0 ; idx_part < my_nb_particles ; ++idx_part){ // TODO remove me const int partition_level = computer.pbc_field_layer(my_particles_positions[idx_part*size_particle_positions+IDXC_Z], IDXC_Z); variable_used_only_in_assert(partition_level); assert(partition_level >= current_partition_interval.first); assert(partition_level < current_partition_interval.second); } particles_utils::partition_extra_z<partsize_t, size_particle_positions>(&my_particles_positions[0], my_nb_particles, partition_interval_size, current_my_nb_particles_per_partition.get(), current_offset_particles_for_partition.get(), [&](const real_number& z_pos){ const int partition_level = computer.pbc_field_layer(z_pos, IDXC_Z); assert(current_partition_interval.first <= partition_level && partition_level < current_partition_interval.second); return partition_level - current_partition_interval.first; }, [&](const partsize_t idx1, const partsize_t idx2){ std::swap(my_particles_positions_indexes[idx1], my_particles_positions_indexes[idx2]); for(int idx_rhs = 0 ; idx_rhs < int(my_particles_rhs.size()) ; ++idx_rhs){ for(int idx_val = 0 ; idx_val < size_particle_rhs ; ++idx_val){ std::swap(my_particles_rhs[idx_rhs][idx1*size_particle_rhs + idx_val], my_particles_rhs[idx_rhs][idx2*size_particle_rhs + idx_val]); } } }); {// TODO remove for(int idxPartition = 0 ; idxPartition < partition_interval_size ; ++idxPartition){ assert(current_my_nb_particles_per_partition[idxPartition] == current_offset_particles_for_partition[idxPartition+1] - current_offset_particles_for_partition[idxPartition]); for(partsize_t idx = current_offset_particles_for_partition[idxPartition] ; idx < current_offset_particles_for_partition[idxPartition+1] ; ++idx){ assert(computer.pbc_field_layer(my_particles_positions[idx*size_particle_positions+IDXC_Z], IDXC_Z)-current_partition_interval.first == idxPartition); } } } } void delete_particles_from_indexes(const std::vector<partsize_t>& inIndexToDelete) final { // We consider that the given indexes are here or in our neighbors, // so we exchange them int my_rank; AssertMpi(MPI_Comm_rank(mpi_com, &my_rank)); int nb_processes; AssertMpi(MPI_Comm_size(mpi_com, &nb_processes)); std::set<partsize_t> setIndexToDelete(inIndexToDelete.begin(), inIndexToDelete.end()); if(nb_processes > 1){ const int TopToBottom = 0; const int BottomToTop = 1; partsize_t nbIndexesFromTop = 0; partsize_t nbIndexesFromBottom = 0; partsize_t myNbIndexes = inIndexToDelete.size(); { MPI_Request requests[4]; AssertMpi(MPI_Isend(&myNbIndexes, 1, particles_utils::GetMpiType(partsize_t()), (my_rank+1)%nb_processes, BottomToTop, mpi_com, &requests[0])); AssertMpi(MPI_Isend(&myNbIndexes, 1, particles_utils::GetMpiType(partsize_t()), (my_rank+nb_processes-1)%nb_processes, TopToBottom, mpi_com, &requests[1])); AssertMpi(MPI_Irecv(&nbIndexesFromTop, 1, particles_utils::GetMpiType(partsize_t()), (my_rank+1)%nb_processes, TopToBottom, mpi_com, &requests[2])); AssertMpi(MPI_Irecv(&nbIndexesFromBottom, 1, particles_utils::GetMpiType(partsize_t()), (my_rank+nb_processes-1)%nb_processes, BottomToTop, mpi_com, &requests[3])); AssertMpi(MPI_Waitall(4, requests, MPI_STATUSES_IGNORE)); } { MPI_Request requests[4]; int nbRequests = 0; if(myNbIndexes){ AssertMpi(MPI_Isend(&inIndexToDelete[0], myNbIndexes, particles_utils::GetMpiType(partsize_t()), (my_rank+1)%nb_processes, BottomToTop, mpi_com, &requests[nbRequests++])); AssertMpi(MPI_Isend(&inIndexToDelete[0], myNbIndexes, particles_utils::GetMpiType(partsize_t()), (my_rank+nb_processes-1)%nb_processes, TopToBottom, mpi_com, &requests[nbRequests++])); } std::vector<partsize_t> indexesFromTop(nbIndexesFromTop); std::vector<partsize_t> indexesFromBottom(nbIndexesFromTop); if(nbIndexesFromTop){ AssertMpi(MPI_Irecv(&indexesFromTop[0], int(nbIndexesFromTop), particles_utils::GetMpiType(partsize_t()), (my_rank+1)%nb_processes, TopToBottom, mpi_com, &requests[nbRequests++])); } if(nbIndexesFromBottom){ AssertMpi(MPI_Irecv(&indexesFromBottom[0], int(nbIndexesFromTop), particles_utils::GetMpiType(partsize_t()), (my_rank+nb_processes-1)%nb_processes, BottomToTop, mpi_com, &requests[nbRequests++])); } AssertMpi(MPI_Waitall(nbRequests, requests, MPI_STATUSES_IGNORE)); std::copy( indexesFromTop.begin(), indexesFromTop.end(), std::inserter( setIndexToDelete, setIndexToDelete.end() ) ); std::copy( indexesFromBottom.begin(), indexesFromBottom.end(), std::inserter( setIndexToDelete, setIndexToDelete.end() ) ); } } if(setIndexToDelete.size()){ partsize_t nbDeletedParticles = 0; for(int idxPartition = 0 ; idxPartition < partition_interval_size ; ++idxPartition){ const partsize_t nbDeletedParticlesBefore = nbDeletedParticles; for(partsize_t idx = current_offset_particles_for_partition[idxPartition] ; idx < current_offset_particles_for_partition[idxPartition+1] ; ++idx){ if(setIndexToDelete.find(my_particles_positions_indexes[idx]) != setIndexToDelete.end()){ nbDeletedParticles += 1; } else if(nbDeletedParticles){ my_particles_positions_indexes[idx-nbDeletedParticles] = my_particles_positions_indexes[idx]; for(int idx_pos = 0 ; idx_pos < size_particle_positions ; ++idx_pos){ my_particles_positions[(idx-nbDeletedParticles)*size_particle_positions+idx_pos] = my_particles_positions[idx*size_particle_positions+idx_pos]; } for(int idx_rhs = 0 ; idx_rhs < int(my_particles_rhs.size()) ; ++idx_rhs){ for(int idx_val = 0 ; idx_val < size_particle_rhs ; ++idx_val){ my_particles_rhs[idx_rhs][(idx-nbDeletedParticles)*size_particle_rhs + idx_val] = my_particles_rhs[idx_rhs][idx*size_particle_rhs + idx_val]; } } } } current_offset_particles_for_partition[idxPartition] -= nbDeletedParticlesBefore; } current_offset_particles_for_partition[partition_interval_size] -= nbDeletedParticles; my_nb_particles -= nbDeletedParticles; assert(my_nb_particles == current_offset_particles_for_partition[partition_interval_size]); } AssertMpi(MPI_Allreduce(const_cast<partsize_t*>(&my_nb_particles), &total_nb_particles, 1, particles_utils::GetMpiType(partsize_t()), MPI_SUM, mpi_com)); } void compute() final { TIMEZONE("particles_system::compute"); particles_distr.template compute_distr<computer_class, field_class, size_particle_positions, size_particle_rhs>( computer, default_field, current_my_nb_particles_per_partition.get(), my_particles_positions.get(), my_particles_rhs.front().get(), interp_neighbours); } void compute_p2p() final { // TODO P2P if(computer_p2p.isEnable() == true){ TIMEZONE("particles_system::compute_p2p"); distr_p2p.template compute_distr<p2p_computer_class, size_particle_positions, size_particle_rhs>( computer_p2p, current_my_nb_particles_per_partition.get(), my_particles_positions.get(), my_particles_rhs.data(), int(my_particles_rhs.size()), my_particles_positions_indexes.get()); } } void compute_particles_inner() final { if(computer_particules_inner.isEnable() == true){ TIMEZONE("particles_system::compute_particles_inner"); computer_particules_inner.template compute_interaction<size_particle_positions, size_particle_rhs>( my_nb_particles, my_particles_positions.get(), my_particles_rhs.front().get()); } } void add_Lagrange_multipliers() final { if(computer_particules_inner.isEnable() == true){ TIMEZONE("particles_system::add_Lagrange_multipliers"); computer_particules_inner.template add_Lagrange_multipliers<size_particle_positions, size_particle_rhs>( my_nb_particles, my_particles_positions.get(), my_particles_rhs.front().get()); } } void enforce_unit_orientation() final { if(computer_particules_inner.isEnable() == true){ TIMEZONE("particles_system::enforce_unit_orientation"); computer_particules_inner.template enforce_unit_orientation<size_particle_positions>( my_nb_particles, my_particles_positions.get()); } } void compute_sphere_particles_inner(const real_number particle_extra_field[]) final { if(computer_particules_inner.isEnable() == true){ TIMEZONE("particles_system::compute_sphere_particles_inner"); computer_particules_inner.template compute_interaction_with_extra<size_particle_positions, size_particle_rhs, 3>( my_nb_particles, my_particles_positions.get(), my_particles_rhs.front().get(), particle_extra_field); } } void compute_ellipsoid_particles_inner(const real_number particle_extra_field[]) final { if(computer_particules_inner.isEnable() == true){ TIMEZONE("particles_system::compute_ellipsoid_particles_inner"); computer_particules_inner.template compute_interaction_with_extra<size_particle_positions, size_particle_rhs, 9>( my_nb_particles, my_particles_positions.get(), my_particles_rhs.front().get(), particle_extra_field); } } template <class sample_field_class, int sample_size_particle_rhs> void sample_compute(const sample_field_class& sample_field, real_number sample_rhs[]) { TIMEZONE("particles_system::compute"); particles_distr.template compute_distr<computer_class, sample_field_class, size_particle_positions, sample_size_particle_rhs>( computer, sample_field, current_my_nb_particles_per_partition.get(), my_particles_positions.get(), sample_rhs, interp_neighbours); } //- Not generic to enable sampling begin void sample_compute_field(const field<float, FFTW, ONE>& sample_field, real_number sample_rhs[]) final { // sample_compute<decltype(sample_field), 1>(sample_field, sample_rhs); } void sample_compute_field(const field<float, FFTW, THREE>& sample_field, real_number sample_rhs[]) final { sample_compute<decltype(sample_field), 3>(sample_field, sample_rhs); } void sample_compute_field(const field<float, FFTW, THREExTHREE>& sample_field, real_number sample_rhs[]) final { sample_compute<decltype(sample_field), 9>(sample_field, sample_rhs); } void sample_compute_field(const field<double, FFTW, ONE>& sample_field, real_number sample_rhs[]) final { sample_compute<decltype(sample_field), 1>(sample_field, sample_rhs); } void sample_compute_field(const field<double, FFTW, THREE>& sample_field, real_number sample_rhs[]) final { sample_compute<decltype(sample_field), 3>(sample_field, sample_rhs); } void sample_compute_field(const field<double, FFTW, THREExTHREE>& sample_field, real_number sample_rhs[]) final { sample_compute<decltype(sample_field), 9>(sample_field, sample_rhs); } //- Not generic to enable sampling end void move(const real_number dt) final { TIMEZONE("particles_system::move"); positions_updater.move_particles(my_particles_positions.get(), my_nb_particles, my_particles_rhs.data(), std::min(step_idx, int(my_particles_rhs.size())), dt); } void redistribute() final { TIMEZONE("particles_system::redistribute"); //DEBUG_MSG("step index is %d\n", step_idx); particles_distr.template redistribute<computer_class, size_particle_positions, size_particle_rhs, 1>( computer, current_my_nb_particles_per_partition.get(), &my_nb_particles, &my_particles_positions, my_particles_rhs.data(), int(my_particles_rhs.size()), &my_particles_positions_indexes); } void inc_step_idx() final { step_idx += 1; } int get_step_idx() const final { return step_idx; } void shift_rhs_vectors() final { if(my_particles_rhs.size()){ std::unique_ptr<real_number[]> next_current(std::move(my_particles_rhs.back())); for(int idx_rhs = int(my_particles_rhs.size())-1 ; idx_rhs > 0 ; --idx_rhs){ my_particles_rhs[idx_rhs] = std::move(my_particles_rhs[idx_rhs-1]); } my_particles_rhs[0] = std::move(next_current); particles_utils::memzero(my_particles_rhs[0], size_particle_rhs*my_nb_particles); } } void completeLoop(const real_number dt) final { TIMEZONE("particles_system::completeLoop"); compute(); compute_p2p(); compute_particles_inner(); move(dt); enforce_unit_orientation(); redistribute(); inc_step_idx(); shift_rhs_vectors(); } void complete2ndOrderLoop(const real_number dt) final { assert((size_particle_positions == 6) || (size_particle_positions == 7)); assert(size_particle_rhs == size_particle_positions); std::unique_ptr<real_number[]> sampled_velocity(new real_number[getLocalNbParticles()*3]()); std::fill_n(sampled_velocity.get(), 3*getLocalNbParticles(), 0); sample_compute_field(default_field, sampled_velocity.get()); this->computer_particules_inner.template compute_interaction_with_extra<size_particle_positions, size_particle_rhs, 3>( my_nb_particles, my_particles_positions.get(), my_particles_rhs.front().get(), sampled_velocity.get()); move(dt); redistribute(); inc_step_idx(); shift_rhs_vectors(); } void completeLoopWithVorticity( const real_number dt, const real_number particle_extra_field[]) final { TIMEZONE("particles_system::completeLoopWithVorticity"); compute(); compute_p2p(); compute_sphere_particles_inner(particle_extra_field); move(dt); enforce_unit_orientation(); redistribute(); inc_step_idx(); shift_rhs_vectors(); } void completeLoopWithVelocityGradient( const real_number dt, const real_number particle_extra_field[]) final { TIMEZONE("particles_system::completeLoopWithVelocityGradient"); compute(); compute_p2p(); compute_ellipsoid_particles_inner(particle_extra_field); move(dt); enforce_unit_orientation(); redistribute(); inc_step_idx(); shift_rhs_vectors(); } const real_number* getParticlesState() const final { return my_particles_positions.get(); } std::unique_ptr<real_number[]> extractParticlesState(const int firstState, const int lastState) const final { const int nbStates = std::max(0,(std::min(lastState,size_particle_positions)-firstState)); std::unique_ptr<real_number[]> stateExtract(new real_number[my_nb_particles*nbStates]); for(partsize_t idx_part = 0 ; idx_part < my_nb_particles ; ++idx_part){ for(int idxState = 0 ; idxState < nbStates ; ++idxState){ stateExtract[idx_part*nbStates + idxState] = my_particles_positions[idx_part*size_particle_positions + idxState+firstState]; } } return stateExtract; } const std::unique_ptr<real_number[]>* getParticlesRhs() const final { return my_particles_rhs.data(); } const partsize_t* getParticlesIndexes() const final { return my_particles_positions_indexes.get(); } partsize_t getLocalNbParticles() const final { return my_nb_particles; } partsize_t getGlobalNbParticles() const final { return total_nb_particles; } int getNbRhs() const final { return int(my_particles_rhs.size()); } int setParticleFileLayout(std::vector<hsize_t> input_layout) final{ this->particle_file_layout.resize(input_layout.size()); for (unsigned int i=0; i<this->particle_file_layout.size(); i++) this->particle_file_layout[i] = input_layout[i]; return EXIT_SUCCESS; } std::vector<hsize_t> getParticleFileLayout(void) final{ return std::vector<hsize_t>(this->particle_file_layout); } void checkNan() const { // TODO remove for(partsize_t idx_part = 0 ; idx_part < my_nb_particles ; ++idx_part){ // TODO remove me assert(std::isnan(my_particles_positions[idx_part*size_particle_positions+IDXC_X]) == false); assert(std::isnan(my_particles_positions[idx_part*size_particle_positions+IDXC_Y]) == false); assert(std::isnan(my_particles_positions[idx_part*size_particle_positions+IDXC_Z]) == false); for(int idx_rhs = 0 ; idx_rhs < my_particles_rhs.size() ; ++idx_rhs){ for(int idx_rhs_val = 0 ; idx_rhs_val < size_particle_rhs ; ++idx_rhs_val){ assert(std::isnan(my_particles_rhs[idx_rhs][idx_part*size_particle_rhs+idx_rhs_val]) == false); } } } } const p2p_computer_class& getP2PComputer() const final{ return computer_p2p; } p2p_computer_class& getP2PComputer() final{ return computer_p2p; } const particles_inner_computer_class& getInnerComputer() const{ return computer_particules_inner; } particles_inner_computer_class& getInnnerComputer(){ return computer_particules_inner; } }; #endif