From 690b53a939088b448e60cd853c1a96156ee2e35c Mon Sep 17 00:00:00 2001
From: Berenger Bramas <Berenger.Bramas@inria.fr>
Date: Thu, 17 Oct 2019 15:30:37 +0200
Subject: [PATCH] First version where particles that collide can be removed

---
 cpp/particles/p2p_computer.hpp            |   4 +-
 cpp/particles/p2p_computer_empty.hpp      |   4 +-
 cpp/particles/p2p_cylinder_collisions.hpp |   4 +-
 cpp/particles/p2p_distr_mpi.hpp           | 163 +++++++++++++---------
 cpp/particles/p2p_ghost_collisions.hpp    |   4 +-
 cpp/particles/p2p_merge_collisions.hpp    |  76 ++++++++++
 cpp/particles/particles_system.hpp        |  96 +++++++++++++
 7 files changed, 283 insertions(+), 68 deletions(-)
 create mode 100644 cpp/particles/p2p_merge_collisions.hpp

diff --git a/cpp/particles/p2p_computer.hpp b/cpp/particles/p2p_computer.hpp
index 5fd2aef8..e5ddda47 100644
--- a/cpp/particles/p2p_computer.hpp
+++ b/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{
diff --git a/cpp/particles/p2p_computer_empty.hpp b/cpp/particles/p2p_computer_empty.hpp
index c7c0057f..cb9710cc 100644
--- a/cpp/particles/p2p_computer_empty.hpp
+++ b/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{
diff --git a/cpp/particles/p2p_cylinder_collisions.hpp b/cpp/particles/p2p_cylinder_collisions.hpp
index 46763ddd..c1fc4b97 100644
--- a/cpp/particles/p2p_cylinder_collisions.hpp
+++ b/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*/){
diff --git a/cpp/particles/p2p_distr_mpi.hpp b/cpp/particles/p2p_distr_mpi.hpp
index 64b7b605..e567e86f 100644
--- a/cpp/particles/p2p_distr_mpi.hpp
+++ b/cpp/particles/p2p_distr_mpi.hpp
@@ -49,6 +49,7 @@ protected:
     enum MpiTag{
         TAG_NB_PARTICLES,
         TAG_POSITION_PARTICLES,
+        TAG_INDEX_PARTICLES,
         TAG_RESULT_PARTICLES,
     };
 
@@ -57,10 +58,12 @@ protected:
         int destProc;
         int nbLevelsToExchange;
         bool isRecv;
+        int nbReceived;
 
         std::unique_ptr<real_number[]> toRecvAndMerge;
         std::unique_ptr<real_number[]> toCompute;
         std::unique_ptr<real_number[]> results;
+        std::unique_ptr<partsize_t[]> indexes;
     };
 
     enum Action{
@@ -416,6 +419,7 @@ public:
                 descriptor.nbLevelsToExchange = nb_levels_to_send;
                 descriptor.nbParticlesToExchange = particles_offset_layers[my_nb_cell_levels] - particles_offset_layers[my_nb_cell_levels-nb_levels_to_send];
                 descriptor.isRecv = false;
+                descriptor.nbReceived = 0;
 
                 neigDescriptors.emplace_back(std::move(descriptor));
 
@@ -438,6 +442,7 @@ public:
                 descriptor.nbLevelsToExchange = nb_levels_to_recv;
                 descriptor.nbParticlesToExchange = -1;
                 descriptor.isRecv = true;
+                descriptor.nbReceived = 0;
 
                 neigDescriptors.emplace_back(std::move(descriptor));
 
@@ -479,6 +484,14 @@ public:
                               descriptor.destProc, TAG_POSITION_PARTICLES,
                               current_com, &mpiRequests.back()));
 
+                    whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
+                    mpiRequests.emplace_back();
+                    assert(descriptor.nbParticlesToExchange*size_particle_positions < std::numeric_limits<int>::max());
+                    AssertMpi(MPI_Isend(const_cast<partsize_t*>(&inout_index_particles[particles_offset_layers[my_nb_cell_levels-descriptor.nbLevelsToExchange]]),
+                              int(descriptor.nbParticlesToExchange), particles_utils::GetMpiType(partsize_t()),
+                              descriptor.destProc, TAG_INDEX_PARTICLES,
+                              current_com, &mpiRequests.back()));
+
                     assert(descriptor.toRecvAndMerge == nullptr);
                     descriptor.toRecvAndMerge.reset(new real_number[descriptor.nbParticlesToExchange*size_particle_rhs]);
                     whatNext.emplace_back(std::pair<Action,int>{MERGE_PARTICLES, idxDescr});
@@ -576,6 +589,14 @@ public:
                             AssertMpi(MPI_Irecv(descriptor.toCompute.get(), int(NbParticlesToReceive*size_particle_positions),
                                                 particles_utils::GetMpiType(real_number()), destProc, TAG_POSITION_PARTICLES,
                                                 current_com, &mpiRequests.back()));
+
+                            descriptor.toCompute.reset(new real_number[NbParticlesToReceive*size_particle_positions]);
+                            whatNext.emplace_back(std::pair<Action,int>{COMPUTE_PARTICLES, releasedAction.second});
+                            mpiRequests.emplace_back();
+                            assert(NbParticlesToReceive*size_particle_positions < std::numeric_limits<int>::max());
+                            AssertMpi(MPI_Irecv(descriptor.indexes.get(), int(NbParticlesToReceive),
+                                                particles_utils::GetMpiType(partsize_t()), destProc, TAG_INDEX_PARTICLES,
+                                                current_com, &mpiRequests.back()));
                         }
                     }
 
@@ -583,80 +604,88 @@ public:
                     /// Computation
                     //////////////////////////////////////////////////////////////////////
                     if(releasedAction.first == COMPUTE_PARTICLES){
-                        TIMEZONE("compute-particles");
                         NeighborDescriptor& descriptor = neigDescriptors[releasedAction.second];
-                        assert(descriptor.isRecv);
-                        const partsize_t NbParticlesToReceive = descriptor.nbParticlesToExchange;
-
-                        assert(descriptor.toCompute != nullptr);
-                        descriptor.results.reset(new real_number[NbParticlesToReceive*size_particle_rhs]);
-                        computer_thread.template init_result_array<size_particle_rhs>(descriptor.results.get(), NbParticlesToReceive);
-
-                        // Compute
-                        partsize_t idxPart = 0;
-                        while(idxPart != NbParticlesToReceive){
-                            const long int current_cell_idx = get_cell_idx(descriptor.toCompute[idxPart*size_particle_positions + IDXC_X],
-                                                                           descriptor.toCompute[idxPart*size_particle_positions + IDXC_Y],
-                                                                           descriptor.toCompute[idxPart*size_particle_positions + IDXC_Z]);
-                            partsize_t nb_parts_in_cell = 1;
-                            while(idxPart+nb_parts_in_cell != NbParticlesToReceive
-                                  && current_cell_idx == get_cell_idx(descriptor.toCompute[(idxPart+nb_parts_in_cell)*size_particle_positions + IDXC_X],
-                                                                     descriptor.toCompute[(idxPart+nb_parts_in_cell)*size_particle_positions + IDXC_Y],
-                                                                     descriptor.toCompute[(idxPart+nb_parts_in_cell)*size_particle_positions + IDXC_Z])){
-                                nb_parts_in_cell += 1;
-                            }
+                        descriptor.nbReceived += 1;
+                        assert(descriptor.nbReceived <= 2);
+
+                        if(descriptor.nbReceived == 2){
+                            TIMEZONE("compute-particles");
+                            NeighborDescriptor& descriptor = neigDescriptors[releasedAction.second];
+                            assert(descriptor.isRecv);
+                            const partsize_t NbParticlesToReceive = descriptor.nbParticlesToExchange;
+
+                            assert(descriptor.toCompute != nullptr);
+                            descriptor.results.reset(new real_number[NbParticlesToReceive*size_particle_rhs]);
+                            computer_thread.template init_result_array<size_particle_rhs>(descriptor.results.get(), NbParticlesToReceive);
+
+                            // Compute
+                            partsize_t idxPart = 0;
+                            while(idxPart != NbParticlesToReceive){
+                                const long int current_cell_idx = get_cell_idx(descriptor.toCompute[idxPart*size_particle_positions + IDXC_X],
+                                                                               descriptor.toCompute[idxPart*size_particle_positions + IDXC_Y],
+                                                                               descriptor.toCompute[idxPart*size_particle_positions + IDXC_Z]);
+                                partsize_t nb_parts_in_cell = 1;
+                                while(idxPart+nb_parts_in_cell != NbParticlesToReceive
+                                      && current_cell_idx == get_cell_idx(descriptor.toCompute[(idxPart+nb_parts_in_cell)*size_particle_positions + IDXC_X],
+                                                                         descriptor.toCompute[(idxPart+nb_parts_in_cell)*size_particle_positions + IDXC_Y],
+                                                                         descriptor.toCompute[(idxPart+nb_parts_in_cell)*size_particle_positions + IDXC_Z])){
+                                    nb_parts_in_cell += 1;
+                                }
 
-                            #pragma omp task default(shared) firstprivate(idxPart, nb_parts_in_cell, current_cell_idx)
-                            {
-                                const std::vector<std::pair<partsize_t,partsize_t>>* neighbors[27];
-                                long int neighbors_indexes[27];
-                                std::array<real_number,3> shift[27];
-                                const int nbNeighbors = my_tree.getNeighbors(current_cell_idx, neighbors, neighbors_indexes, shift, true);
-
-                                // with other interval
-                                for(int idx_neighbor = 0 ; idx_neighbor < nbNeighbors ; ++idx_neighbor){
-                                    cells_locker.lock(neighbors_indexes[idx_neighbor]);
-
-                                    for(size_t idx_2 = 0 ; idx_2 < (*neighbors[idx_neighbor]).size() ; ++idx_2){
-                                        for(partsize_t idx_p1 = 0 ; idx_p1 < nb_parts_in_cell ; ++idx_p1){
-                                            for(partsize_t idx_p2 = 0 ; idx_p2 < (*neighbors[idx_neighbor])[idx_2].second ; ++idx_p2){
-                                                const real_number dist_r2 = compute_distance_r2(descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions + IDXC_X],
-                                                                                                descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions + IDXC_Y],
-                                                                                                descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions + IDXC_Z],
-                                                                                                particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDXC_X],
-                                                                                                particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDXC_Y],
-                                                                                                particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDXC_Z],
-                                                                                                shift[idx_neighbor][IDXC_X], shift[idx_neighbor][IDXC_Y], shift[idx_neighbor][IDXC_Z]);
-                                                if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
-                                                    computer_thread.template compute_interaction<size_particle_positions, size_particle_rhs>(
-                                                                        &descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions],
-                                                                        &descriptor.results[(idxPart+idx_p1)*size_particle_rhs],
-                                                                        &particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions],
-                                                                        &particles_current_rhs[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_rhs],
-                                                                        dist_r2, cutoff_radius_compute, shift[idx_neighbor][IDXC_X], shift[idx_neighbor][IDXC_Y], shift[idx_neighbor][IDXC_Z]);
+                                #pragma omp task default(shared) firstprivate(idxPart, nb_parts_in_cell, current_cell_idx)
+                                {
+                                    const std::vector<std::pair<partsize_t,partsize_t>>* neighbors[27];
+                                    long int neighbors_indexes[27];
+                                    std::array<real_number,3> shift[27];
+                                    const int nbNeighbors = my_tree.getNeighbors(current_cell_idx, neighbors, neighbors_indexes, shift, true);
+
+                                    // with other interval
+                                    for(int idx_neighbor = 0 ; idx_neighbor < nbNeighbors ; ++idx_neighbor){
+                                        cells_locker.lock(neighbors_indexes[idx_neighbor]);
+
+                                        for(size_t idx_2 = 0 ; idx_2 < (*neighbors[idx_neighbor]).size() ; ++idx_2){
+                                            for(partsize_t idx_p1 = 0 ; idx_p1 < nb_parts_in_cell ; ++idx_p1){
+                                                for(partsize_t idx_p2 = 0 ; idx_p2 < (*neighbors[idx_neighbor])[idx_2].second ; ++idx_p2){
+                                                    const real_number dist_r2 = compute_distance_r2(descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions + IDXC_X],
+                                                                                                    descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions + IDXC_Y],
+                                                                                                    descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions + IDXC_Z],
+                                                                                                    particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDXC_X],
+                                                                                                    particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDXC_Y],
+                                                                                                    particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDXC_Z],
+                                                                                                    shift[idx_neighbor][IDXC_X], shift[idx_neighbor][IDXC_Y], shift[idx_neighbor][IDXC_Z]);
+                                                    if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
+                                                        computer_thread.template compute_interaction<size_particle_positions, size_particle_rhs>(
+                                                                            descriptor.indexes[(idxPart+idx_p1)],
+                                                                            &descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions],
+                                                                            &descriptor.results[(idxPart+idx_p1)*size_particle_rhs],
+                                                                            inout_index_particles[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)],
+                                                                            &particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions],
+                                                                            &particles_current_rhs[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_rhs],
+                                                                            dist_r2, cutoff_radius_compute, shift[idx_neighbor][IDXC_X], shift[idx_neighbor][IDXC_Y], shift[idx_neighbor][IDXC_Z]);
+                                                    }
                                                 }
                                             }
                                         }
-                                    }
 
-                                    cells_locker.unlock(neighbors_indexes[idx_neighbor]);
+                                        cells_locker.unlock(neighbors_indexes[idx_neighbor]);
+                                    }
                                 }
-                            }
 
-                            idxPart += nb_parts_in_cell;
-                        }
+                                idxPart += nb_parts_in_cell;
+                            }
 
-                        #pragma omp taskwait
+                            #pragma omp taskwait
 
-                        // Send back
-                        const int destProc = descriptor.destProc;
-                        whatNext.emplace_back(std::pair<Action,int>{RELEASE_BUFFER_PARTICLES, releasedAction.second});
-                        mpiRequests.emplace_back();
-                        assert(NbParticlesToReceive*size_particle_rhs < std::numeric_limits<int>::max());
-                        AssertMpi(MPI_Isend(descriptor.results.get(), int(NbParticlesToReceive*size_particle_rhs),
-                                            particles_utils::GetMpiType(real_number()), destProc, TAG_RESULT_PARTICLES,
-                                            current_com, &mpiRequests.back()));
-                        delete[] descriptor.toCompute.release();
+                            // Send back
+                            const int destProc = descriptor.destProc;
+                            whatNext.emplace_back(std::pair<Action,int>{RELEASE_BUFFER_PARTICLES, releasedAction.second});
+                            mpiRequests.emplace_back();
+                            assert(NbParticlesToReceive*size_particle_rhs < std::numeric_limits<int>::max());
+                            AssertMpi(MPI_Isend(descriptor.results.get(), int(NbParticlesToReceive*size_particle_rhs),
+                                                particles_utils::GetMpiType(real_number()), destProc, TAG_RESULT_PARTICLES,
+                                                current_com, &mpiRequests.back()));
+                            delete[] descriptor.toCompute.release();
+                        }
                     }
                     //////////////////////////////////////////////////////////////////////
                     /// Release memory that was sent back
@@ -713,8 +742,10 @@ public:
                                                                                 0, 0, 0);
                                 if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
                                     computer_thread.template compute_interaction<size_particle_positions,size_particle_rhs>(
+                                                        inout_index_particles[(intervals[idx_1].first+idx_p1)],
                                                         &particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions],
                                                         &particles_current_rhs[(intervals[idx_1].first+idx_p1)*size_particle_rhs],
+                                                        inout_index_particles[(intervals[idx_1].first+idx_p2)],
                                                         &particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions],
                                                         &particles_current_rhs[(intervals[idx_1].first+idx_p2)*size_particle_rhs],
                                                         dist_r2, cutoff_radius_compute, 0, 0, 0);
@@ -735,8 +766,10 @@ public:
                                                                                     0, 0, 0);
                                     if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
                                         computer_thread.template compute_interaction<size_particle_positions,size_particle_rhs>(
+                                                            inout_index_particles[(intervals[idx_1].first+idx_p1)],
                                                             &particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions],
                                                             &particles_current_rhs[(intervals[idx_1].first+idx_p1)*size_particle_rhs],
+                                                            inout_index_particles[(intervals[idx_2].first+idx_p2)],
                                                             &particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions],
                                                             &particles_current_rhs[(intervals[idx_2].first+idx_p2)*size_particle_rhs],
                                                             dist_r2, cutoff_radius_compute, 0, 0, 0);
@@ -769,8 +802,10 @@ public:
                                                                                             shift[idx_neighbor][IDXC_X], shift[idx_neighbor][IDXC_Y], shift[idx_neighbor][IDXC_Z]);
                                             if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
                                                 computer_thread.template compute_interaction<size_particle_positions,size_particle_rhs>(
+                                                                    inout_index_particles[(intervals[idx_1].first+idx_p1)],
                                                                     &particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions],
                                                                     &particles_current_rhs[(intervals[idx_1].first+idx_p1)*size_particle_rhs],
+                                                                    inout_index_particles[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)],
                                                                     &particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions],
                                                                     &particles_current_rhs[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_rhs],
                                                                     dist_r2, cutoff_radius_compute, shift[idx_neighbor][IDXC_X], shift[idx_neighbor][IDXC_Y], shift[idx_neighbor][IDXC_Z]);
diff --git a/cpp/particles/p2p_ghost_collisions.hpp b/cpp/particles/p2p_ghost_collisions.hpp
index 604ac120..63c46981 100644
--- a/cpp/particles/p2p_ghost_collisions.hpp
+++ b/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*/){
diff --git a/cpp/particles/p2p_merge_collisions.hpp b/cpp/particles/p2p_merge_collisions.hpp
new file mode 100644
index 00000000..96fff9c2
--- /dev/null
+++ b/cpp/particles/p2p_merge_collisions.hpp
@@ -0,0 +1,76 @@
+/******************************************************************************
+*                                                                             *
+*  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
diff --git a/cpp/particles/particles_system.hpp b/cpp/particles/particles_system.hpp
index cdb523bb..8756aec8 100644
--- a/cpp/particles/particles_system.hpp
+++ b/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>(
-- 
GitLab