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 1/4] 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
From 922db8baa992d70c179a7b4ecc76900b34706930 Mon Sep 17 00:00:00 2001
From: Berenger Bramas <Berenger.Bramas@inria.fr>
Date: Thu, 17 Oct 2019 15:44:49 +0200
Subject: [PATCH 2/4] Update the interface and set use
---
cpp/particles/abstract_particles_system.hpp | 3 +++
cpp/particles/particles_system.hpp | 24 +++++++++++++--------
2 files changed, 18 insertions(+), 9 deletions(-)
diff --git a/cpp/particles/abstract_particles_system.hpp b/cpp/particles/abstract_particles_system.hpp
index 2f2f510f..8a1be218 100644
--- a/cpp/particles/abstract_particles_system.hpp
+++ b/cpp/particles/abstract_particles_system.hpp
@@ -27,6 +27,7 @@
#define ABSTRACT_PARTICLES_SYSTEM_HPP
#include <memory>
+#include <vector>
//- Not generic to enable sampling begin
#include "field.hpp"
@@ -39,6 +40,8 @@ class abstract_particles_system {
public:
virtual ~abstract_particles_system(){}
+ virtual void delete_particles_from_indexes(const std::vector<partsize_t>& inIndexToDelete) = 0;
+
virtual void compute() = 0;
virtual void compute_p2p() = 0;
diff --git a/cpp/particles/particles_system.hpp b/cpp/particles/particles_system.hpp
index 8756aec8..7c746965 100644
--- a/cpp/particles/particles_system.hpp
+++ b/cpp/particles/particles_system.hpp
@@ -28,6 +28,7 @@
#include <array>
#include <set>
+#include <algorithm>
#include "abstract_particles_system.hpp"
#include "abstract_particles_system_with_p2p.hpp"
@@ -73,7 +74,7 @@ class particles_system : public abstract_particles_system_with_p2p<partsize_t, r
std::unique_ptr<real_number[]> my_particles_positions;
std::unique_ptr<partsize_t[]> my_particles_positions_indexes;
partsize_t my_nb_particles;
- const partsize_t total_nb_particles;
+ partsize_t total_nb_particles;
std::vector<std::unique_ptr<real_number[]>> my_particles_rhs;
std::vector<hsize_t> particle_file_layout;
@@ -162,7 +163,7 @@ public:
}
}
- void delete_particles_from_indexes(const std::vector<partsize_t>& inIndexToDelete){
+ 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;
@@ -170,6 +171,8 @@ public:
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;
@@ -191,7 +194,6 @@ public:
AssertMpi(MPI_Waitall(4, requests, MPI_STATUSES_IGNORE));
}
- inIndexToDelete.resize(inIndexToDelete.size() + nbIndexesFromTop + nbIndexesFromBottom);
{
MPI_Request requests[4];
int nbRequests = 0;
@@ -203,24 +205,28 @@ public:
(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(&inIndexToDelete[myNbIndexes], int(nbIndexesFromTop), particles_utils::GetMpiType(partsize_t()),
+ 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(&inIndexToDelete[myNbIndexes+nbIndexesFromTop], int(nbIndexesFromTop), particles_utils::GetMpiType(partsize_t()),
+ 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(inIndexToDelete.size()){
+ if(setIndexToDelete.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){
@@ -253,7 +259,7 @@ public:
assert(my_nb_particles == current_offset_particles_for_partition[partition_interval_size]);
}
- AssertMpi(MPI_Allreduce(&my_nb_particles, &total_nb_particles, 1,
+ AssertMpi(MPI_Allreduce(const_cast<partsize_t*>(&my_nb_particles), &total_nb_particles, 1,
particles_utils::GetMpiType(partsize_t()), MPI_SUM, mpi_com));
}
--
GitLab
From 630e7e8aed890ef03f3c1c4ae723b9acfbda2a2a Mon Sep 17 00:00:00 2001
From: Berenger Bramas <Berenger.Bramas@inria.fr>
Date: Mon, 4 Nov 2019 16:46:12 +0100
Subject: [PATCH 3/4] Debug (forgot to allocate the right array for indexes)
---
cpp/particles/p2p_distr_mpi.hpp | 5 ++++-
1 file changed, 4 insertions(+), 1 deletion(-)
diff --git a/cpp/particles/p2p_distr_mpi.hpp b/cpp/particles/p2p_distr_mpi.hpp
index e567e86f..274dc8a6 100644
--- a/cpp/particles/p2p_distr_mpi.hpp
+++ b/cpp/particles/p2p_distr_mpi.hpp
@@ -580,6 +580,7 @@ public:
const partsize_t NbParticlesToReceive = descriptor.nbParticlesToExchange;
assert(NbParticlesToReceive != -1);
assert(descriptor.toCompute == nullptr);
+ assert(descriptor.indexes == nullptr);
if(NbParticlesToReceive){
descriptor.toCompute.reset(new real_number[NbParticlesToReceive*size_particle_positions]);
@@ -590,7 +591,7 @@ public:
particles_utils::GetMpiType(real_number()), destProc, TAG_POSITION_PARTICLES,
current_com, &mpiRequests.back()));
- descriptor.toCompute.reset(new real_number[NbParticlesToReceive*size_particle_positions]);
+ descriptor.indexes.reset(new partsize_t[NbParticlesToReceive]);
whatNext.emplace_back(std::pair<Action,int>{COMPUTE_PARTICLES, releasedAction.second});
mpiRequests.emplace_back();
assert(NbParticlesToReceive*size_particle_positions < std::numeric_limits<int>::max());
@@ -615,6 +616,7 @@ public:
const partsize_t NbParticlesToReceive = descriptor.nbParticlesToExchange;
assert(descriptor.toCompute != nullptr);
+ assert(descriptor.indexes != nullptr);
descriptor.results.reset(new real_number[NbParticlesToReceive*size_particle_rhs]);
computer_thread.template init_result_array<size_particle_rhs>(descriptor.results.get(), NbParticlesToReceive);
@@ -685,6 +687,7 @@ public:
particles_utils::GetMpiType(real_number()), destProc, TAG_RESULT_PARTICLES,
current_com, &mpiRequests.back()));
delete[] descriptor.toCompute.release();
+ delete[] descriptor.indexes.release();
}
}
//////////////////////////////////////////////////////////////////////
--
GitLab
From 395528bd2b244b5812ce0adf8bd1f5148a89cccf Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@ds.mpg.de>
Date: Tue, 5 Nov 2019 16:05:31 +0100
Subject: [PATCH 4/4] adds hdf5 single value write method
---
cpp/hdf5_tools.cpp | 42 ++++++++++++++++++++++++++++++++++++++++++
cpp/hdf5_tools.hpp | 6 ++++++
2 files changed, 48 insertions(+)
diff --git a/cpp/hdf5_tools.cpp b/cpp/hdf5_tools.cpp
index 5a3aef39..9300f7ef 100644
--- a/cpp/hdf5_tools.cpp
+++ b/cpp/hdf5_tools.cpp
@@ -169,6 +169,36 @@ number hdf5_tools::read_value(
return result;
}
+template <typename number>
+int hdf5_tools::write_value_with_single_rank(
+ const hid_t group,
+ const std::string dset_name,
+ const number value)
+{
+ hid_t dset, mem_dtype;
+ if (typeid(number) == typeid(int))
+ mem_dtype = H5Tcopy(H5T_NATIVE_INT);
+ else if (typeid(number) == typeid(double))
+ mem_dtype = H5Tcopy(H5T_NATIVE_DOUBLE);
+ if (H5Lexists(group, dset_name.c_str(), H5P_DEFAULT))
+ {
+ dset = H5Dopen(group, dset_name.c_str(), H5P_DEFAULT);
+ }
+ else
+ {
+ hid_t fspace;
+ hsize_t count[1];
+ count[0] = 1;
+ fspace = H5Screate_simple(1, count, NULL);
+ dset = H5Dcreate(group, dset_name.c_str(), mem_dtype, fspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ H5Sclose(fspace);
+ }
+ H5Dwrite(dset, mem_dtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &value);
+ H5Dclose(dset);
+ H5Tclose(mem_dtype);
+ return EXIT_SUCCESS;
+}
+
template <typename dtype>
std::vector<dtype> hdf5_tools::read_vector_with_single_rank(
const int myrank,
@@ -269,3 +299,15 @@ double hdf5_tools::read_value<double>(
const hid_t,
const std::string);
+template
+int hdf5_tools::write_value_with_single_rank<int>(
+ const hid_t group,
+ const std::string dset_name,
+ const int value);
+
+template
+int hdf5_tools::write_value_with_single_rank<double>(
+ const hid_t group,
+ const std::string dset_name,
+ const double value);
+
diff --git a/cpp/hdf5_tools.hpp b/cpp/hdf5_tools.hpp
index 99ba45a1..82c00486 100644
--- a/cpp/hdf5_tools.hpp
+++ b/cpp/hdf5_tools.hpp
@@ -84,6 +84,12 @@ namespace hdf5_tools
number read_value(
const hid_t group,
const std::string dset_name);
+
+ template <typename number>
+ int write_value_with_single_rank(
+ const hid_t group,
+ const std::string dset_name,
+ const number value);
}
#endif//HDF5_TOOLS_HPP
--
GitLab