First version where particles that collide can be removed

cpp/particles/p2p_computer.hpp

@@ -70,7 +70,9 @@ public:
     }
 
     template <int size_particle_positions, int size_particle_rhs>
-    void compute_interaction(const real_number pos_part1[], real_number rhs_part1[],
+    void compute_interaction(const partsize_t /*idx_part1*/,
+                             const real_number pos_part1[], real_number rhs_part1[],
+                             const partsize_t /*idx_part2*/,
                              const real_number pos_part2[], real_number rhs_part2[],
                              const real_number dist_pow2, const real_number cutoff,
                              const real_number /*xshift_coef*/, const real_number /*yshift_coef*/, const real_number /*zshift_coef*/) const{

cpp/particles/p2p_computer_empty.hpp

@@ -40,7 +40,9 @@ public:
     }
 
     template <int size_particle_positions, int size_particle_rhs>
-    void compute_interaction(const real_number /*pos_part1*/[], real_number /*rhs_part1*/[],
+    void compute_interaction(const partsize_t /*idx_part1*/,
+                             const real_number /*pos_part1*/[], real_number /*rhs_part1*/[],
+                             const partsize_t /*idx_part2*/,
                              const real_number /*pos_part2*/[], real_number /*rhs_part2*/[],
                              const real_number /*dist_pow2*/,  const real_number /*cutoff*/,
                              const real_number /*xshift_coef*/, const real_number /*yshift_coef*/, const real_number /*zshift_coef*/) const{

cpp/particles/p2p_cylinder_collisions.hpp

@@ -44,7 +44,9 @@ public:
     }
 
     template <int size_particle_positions, int size_particle_rhs>
-    void compute_interaction(const real_number /*pos_part1*/[], real_number /*rhs_part1*/[],
+    void compute_interaction(const partsize_t /*idx_part1*/,
+                             const real_number /*pos_part1*/[], real_number /*rhs_part1*/[],
+                             const partsize_t /*idx_part2*/,
                              const real_number /*pos_part2*/[], real_number /*rhs_part2*/[],
                              const real_number /*dist_pow2*/,  const real_number /*cutoff*/,
                              const real_number /*xshift_coef*/, const real_number /*yshift_coef*/, const real_number /*zshift_coef*/){

cpp/particles/p2p_ghost_collisions.hpp

@@ -44,7 +44,9 @@ public:
     }
 
     template <int size_particle_positions, int size_particle_rhs>
-    void compute_interaction(const real_number /*pos_part1*/[], real_number /*rhs_part1*/[],
+    void compute_interaction(const partsize_t /*idx_part1*/,
+                             const real_number /*pos_part1*/[], real_number /*rhs_part1*/[],
+                             const partsize_t /*idx_part2*/,
                              const real_number /*pos_part2*/[], real_number /*rhs_part2*/[],
                              const real_number /*dist_pow2*/,  const real_number /*cutoff*/,
                              const real_number /*xshift_coef*/, const real_number /*yshift_coef*/, const real_number /*zshift_coef*/){

cpp/particles/p2p_merge_collisions.hpp

+/******************************************************************************
+*                                                                             *
+*  Copyright 2019 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                                         *
+*                                                                             *
+******************************************************************************/
+#ifndef P2P_MERGE_COLLISIONS_HPP
+#define P2P_MERGE_COLLISIONS_HPP
+
+#include <cstring>
+
+template <class real_number, class partsize_t>
+class p2p_merge_collisions{
+    long int collision_counter;
+
+    std::vector<partsize_t> mergedParticles;
+
+public:
+    p2p_merge_collisions() {}
+
+    // Copy constructor use a counter set to zero
+    p2p_merge_collisions(const p2p_merge_collisions&){}
+
+    template <int size_particle_rhs>
+    void init_result_array(real_number /*rhs*/[], const partsize_t /*nbParticles*/) const{
+    }
+
+    template <int size_particle_rhs>
+    void reduce_particles_rhs(real_number /*rhs_dst*/[], const real_number /*rhs_src*/[], const partsize_t /*nbParticles*/) const{
+    }
+
+    template <int size_particle_positions, int size_particle_rhs>
+    void compute_interaction(const partsize_t idx_part1,
+                             const real_number /*pos_part1*/[], real_number /*rhs_part1*/[],
+                             const partsize_t idx_part2,
+                             const real_number /*pos_part2*/[], real_number /*rhs_part2*/[],
+                             const real_number /*dist_pow2*/,  const real_number /*cutoff*/,
+                             const real_number /*xshift_coef*/, const real_number /*yshift_coef*/, const real_number /*zshift_coef*/){
+        mergedParticles.emplace_back(std::max(idx_part1,idx_part2));
+    }
+
+    void merge(const p2p_merge_collisions& other){
+        collision_counter.insert(collision_counter.end(), other.collision_counter.begin(), other.collision_counter.end());
+    }
+
+    constexpr static bool isEnable() {
+        return true;
+    }
+
+    auto& get_merge_list() const{
+        return collision_counter;
+    }
+
+    void reset_merge_list(){
+        mergedParticles.clear();
+    }
+};
+
+
+#endif // P2P_GHOST_COLLISIONS_HPP

cpp/particles/particles_system.hpp

@@ -27,6 +27,7 @@
 #define PARTICLES_SYSTEM_HPP
 
 #include <array>
+#include <set>
 
 #include "abstract_particles_system.hpp"
 #include "abstract_particles_system_with_p2p.hpp"
@@ -161,6 +162,101 @@ public:
         }
     }
 
+    void delete_particles_from_indexes(const std::vector<partsize_t>& inIndexToDelete){
+        // 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));
+
+        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));
+            }
+            inIndexToDelete.resize(inIndexToDelete.size() + nbIndexesFromTop + nbIndexesFromBottom);
+            {
+                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++]));
+                }
+
+                if(nbIndexesFromTop){
+                    AssertMpi(MPI_Irecv(&inIndexToDelete[myNbIndexes], int(nbIndexesFromTop), particles_utils::GetMpiType(partsize_t()),
+                              (my_rank+1)%nb_processes, TopToBottom, mpi_com, &requests[nbRequests++]));
+                }
+                if(nbIndexesFromBottom){
+                    AssertMpi(MPI_Irecv(&inIndexToDelete[myNbIndexes+nbIndexesFromTop], 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));
+            }
+        }
+
+        if(inIndexToDelete.size()){
+            partsize_t nbDeletedParticles = 0;
+
+            std::set<partsize_t> setIndexToDelete(inIndexToDelete.begin(), inIndexToDelete.end());
+
+            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(&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>(