From 99d5fe727a8825da6f17693196e2e4f3e9866d77 Mon Sep 17 00:00:00 2001
From: Berenger Bramas <berenger.bramas@mpcdf.mpg.de>
Date: Fri, 22 Sep 2017 15:53:10 +0200
Subject: [PATCH] Update BFPS interface to include the new P2P classes
---
bfps/cpp/full_code/NSVEparticlesP2P.cpp | 113 ++++++++++++++++++
bfps/cpp/full_code/NSVEparticlesP2P.hpp | 82 +++++++++++++
.../particles/abstract_particles_system.hpp | 2 +
bfps/cpp/particles/lock_free_bool_array.hpp | 4 +-
bfps/cpp/particles/p2p_computer.hpp | 6 +
bfps/cpp/particles/p2p_computer_empty.hpp | 31 +++++
bfps/cpp/particles/p2p_distr_mpi.hpp | 6 +-
bfps/cpp/particles/particles_distr_mpi.hpp | 112 ++++++++++++++---
bfps/cpp/particles/particles_system.hpp | 37 +++++-
.../particles/particles_system_builder.hpp | 30 ++++-
setup.py | 3 +-
11 files changed, 393 insertions(+), 33 deletions(-)
create mode 100644 bfps/cpp/full_code/NSVEparticlesP2P.cpp
create mode 100644 bfps/cpp/full_code/NSVEparticlesP2P.hpp
create mode 100644 bfps/cpp/particles/p2p_computer_empty.hpp
diff --git a/bfps/cpp/full_code/NSVEparticlesP2P.cpp b/bfps/cpp/full_code/NSVEparticlesP2P.cpp
new file mode 100644
index 00000000..08326119
--- /dev/null
+++ b/bfps/cpp/full_code/NSVEparticlesP2P.cpp
@@ -0,0 +1,113 @@
+#include <string>
+#include <cmath>
+#include "NSVEparticlesP2P.hpp"
+#include "scope_timer.hpp"
+#include "particles/particles_sampling.hpp"
+#include "particles/p2p_computer.hpp"
+
+template <typename rnumber>
+int NSVEparticlesP2P<rnumber>::initialize(void)
+{
+ this->NSVE<rnumber>::initialize();
+
+ this->ps = particles_system_builder(
+ this->fs->cvelocity, // (field object)
+ this->fs->kk, // (kspace object, contains dkx, dky, dkz)
+ tracers0_integration_steps, // to check coherency between parameters and hdf input file (nb rhs)
+ (long long int)nparticles, // to check coherency between parameters and hdf input file
+ this->fs->get_current_fname(), // particles input filename
+ std::string("/tracers0/state/") + std::to_string(this->fs->iteration), // dataset name for initial input
+ std::string("/tracers0/rhs/") + std::to_string(this->fs->iteration), // dataset name for initial input
+ tracers0_neighbours, // parameter (interpolation no neighbours)
+ tracers0_smoothness, // parameter
+ this->comm,
+ this->fs->iteration+1);
+ // TODO P2P write particle data too
+ this->particles_output_writer_mpi = new particles_output_hdf5<
+ long long int, double, 3, 3>(
+ MPI_COMM_WORLD,
+ "tracers0",
+ nparticles,
+ tracers0_integration_steps);
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVEparticlesP2P<rnumber>::step(void)
+{
+ this->fs->compute_velocity(this->fs->cvorticity);
+ this->fs->cvelocity->ift();
+ this->ps->completeLoop(this->dt);
+ this->NSVE<rnumber>::step();
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVEparticlesP2P<rnumber>::write_checkpoint(void)
+{
+ this->NSVE<rnumber>::write_checkpoint();
+ this->particles_output_writer_mpi->open_file(this->fs->get_current_fname());
+ // TODO P2P write particle data too
+ this->particles_output_writer_mpi->save(
+ this->ps->getParticlesPositions(),
+ this->ps->getParticlesRhs(),
+ this->ps->getParticlesIndexes(),
+ this->ps->getLocalNbParticles(),
+ this->fs->iteration);
+ this->particles_output_writer_mpi->close_file();
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVEparticlesP2P<rnumber>::finalize(void)
+{
+ this->ps.release();
+ delete this->particles_output_writer_mpi;
+ this->NSVE<rnumber>::finalize();
+ return EXIT_SUCCESS;
+}
+
+/** \brief Compute fluid stats and sample fields at particle locations.
+ */
+
+template <typename rnumber>
+int NSVEparticlesP2P<rnumber>::do_stats()
+{
+ /// fluid stats go here
+ this->NSVE<rnumber>::do_stats();
+
+
+ if (!(this->iteration % this->niter_part == 0))
+ return EXIT_SUCCESS;
+
+ /// sample position
+ sample_particles_system_position(
+ this->ps,
+ (this->simname + "_particles.h5"), // filename
+ "tracers0", // hdf5 parent group
+ "position" // dataset basename TODO
+ );
+
+ /// sample velocity
+ sample_from_particles_system(*this->tmp_vec_field, // field to save
+ this->ps,
+ (this->simname + "_particles.h5"), // filename
+ "tracers0", // hdf5 parent group
+ "velocity" // dataset basename TODO
+ );
+
+ /// compute acceleration and sample it
+ this->fs->compute_Lagrangian_acceleration(this->tmp_vec_field);
+ this->tmp_vec_field->ift();
+ sample_from_particles_system(*this->tmp_vec_field,
+ this->ps,
+ (this->simname + "_particles.h5"),
+ "tracers0",
+ "acceleration");
+
+ return EXIT_SUCCESS;
+}
+
+template class NSVEparticlesP2P<float>;
+template class NSVEparticlesP2P<double>;
+
diff --git a/bfps/cpp/full_code/NSVEparticlesP2P.hpp b/bfps/cpp/full_code/NSVEparticlesP2P.hpp
new file mode 100644
index 00000000..9b015659
--- /dev/null
+++ b/bfps/cpp/full_code/NSVEparticlesP2P.hpp
@@ -0,0 +1,82 @@
+/**********************************************************************
+* *
+* Copyright 2017 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 NSVEPARTICLESP2P_HPP
+#define NSVEPARTICLESP2P_HPP
+
+
+
+#include <cstdlib>
+#include "base.hpp"
+#include "vorticity_equation.hpp"
+#include "full_code/NSVE.hpp"
+#include "particles/particles_system_builder.hpp"
+#include "particles/particles_output_hdf5.hpp"
+
+/** \brief Navier-Stokes solver that includes simple Lagrangian tracers.
+ *
+ * Child of Navier Stokes vorticity equation solver, this class calls all the
+ * methods from `NSVE`, and in addition integrates simple Lagrangian tracers
+ * in the resulting velocity field.
+ */
+
+template <typename rnumber>
+class NSVEparticlesP2P: public NSVE<rnumber>
+{
+ public:
+
+ /* parameters that are read in read_parameters */
+ int niter_part;
+ int nparticles;
+ int tracers0_integration_steps;
+ int tracers0_neighbours;
+ int tracers0_smoothness;
+
+ /* other stuff */
+ std::unique_ptr<abstract_particles_system<long long int, double>> ps;
+ // TODO P2P use a reader with particle data
+ particles_output_hdf5<long long int, double,3,3> *particles_output_writer_mpi;
+
+
+ NSVEparticlesP2P(
+ const MPI_Comm COMMUNICATOR,
+ const std::string &simulation_name):
+ NSVE<rnumber>(
+ COMMUNICATOR,
+ simulation_name){}
+ ~NSVEparticlesP2P(){}
+
+ int initialize(void);
+ int step(void);
+ int finalize(void);
+
+ int read_parameters(void);
+ int write_checkpoint(void);
+ int do_stats(void);
+};
+
+#endif//NSVEPARTICLESP2P_HPP
+
diff --git a/bfps/cpp/particles/abstract_particles_system.hpp b/bfps/cpp/particles/abstract_particles_system.hpp
index 1c8592f3..26ce7198 100644
--- a/bfps/cpp/particles/abstract_particles_system.hpp
+++ b/bfps/cpp/particles/abstract_particles_system.hpp
@@ -14,6 +14,8 @@ class abstract_particles_system {
public:
virtual void compute() = 0;
+ virtual void compute_p2p() = 0;
+
virtual void move(const real_number dt) = 0;
virtual void redistribute() = 0;
diff --git a/bfps/cpp/particles/lock_free_bool_array.hpp b/bfps/cpp/particles/lock_free_bool_array.hpp
index d0b2aa2b..928c1753 100644
--- a/bfps/cpp/particles/lock_free_bool_array.hpp
+++ b/bfps/cpp/particles/lock_free_bool_array.hpp
@@ -15,7 +15,7 @@ public:
}
}
- void lock(const int inKey){
+ void lock(const long int inKey){
volatile long int* k = keys[inKey%keys.size()].get();
long int res = 1;
while(res == 1){
@@ -23,7 +23,7 @@ public:
}
}
- void unlock(const int inKey){
+ void unlock(const long int inKey){
volatile long int* k = keys[inKey%keys.size()].get();
assert(k && *k);
(*k) = 0;
diff --git a/bfps/cpp/particles/p2p_computer.hpp b/bfps/cpp/particles/p2p_computer.hpp
index cb384950..de024ff7 100644
--- a/bfps/cpp/particles/p2p_computer.hpp
+++ b/bfps/cpp/particles/p2p_computer.hpp
@@ -6,6 +6,8 @@
template <class real_number, class partsize_t>
class p2p_computer{
public:
+ constexpr static int size_data = 3;
+
template <int size_particle_rhs>
void init_result_array(real_number rhs[], const partsize_t nbParticles) const{
memset(rhs, 0, sizeof(real_number)*nbParticles*size_particle_rhs);
@@ -28,6 +30,10 @@ public:
// TODO put the kernel here
static_assert(size_particle_positions == 3, "This kernel works only with 3 values for one position");
}
+
+ constexpr static bool isEmpty() {
+ return false;
+ }
};
#endif
diff --git a/bfps/cpp/particles/p2p_computer_empty.hpp b/bfps/cpp/particles/p2p_computer_empty.hpp
new file mode 100644
index 00000000..7076061e
--- /dev/null
+++ b/bfps/cpp/particles/p2p_computer_empty.hpp
@@ -0,0 +1,31 @@
+#ifndef P2P_COMPUTER_EMPTY_HPP
+#define P2P_COMPUTER_EMPTY_HPP
+
+#include <cstring>
+
+template <class real_number, class partsize_t>
+class p2p_computer_empty{
+public:
+ constexpr int static size_data = 0;
+
+ 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_data, int size_particle_rhs>
+ void compute_interaction(const real_number /*pos_part1*/[], const real_number /*data_part1*/[], real_number /*rhs_part1*/[],
+ const real_number /*pos_part2*/[], const real_number /*data_part2*/[], real_number /*rhs_part2*/[],
+ const real_number /*dist_pow2*/,
+ const real_number /*xshift_coef*/, const real_number /*yshift_coef*/, const real_number /*zshift_coef*/) const{
+ }
+
+ constexpr static bool isEmpty() {
+ return true;
+ }
+};
+
+#endif
diff --git a/bfps/cpp/particles/p2p_distr_mpi.hpp b/bfps/cpp/particles/p2p_distr_mpi.hpp
index 12b349dc..ccc236d6 100644
--- a/bfps/cpp/particles/p2p_distr_mpi.hpp
+++ b/bfps/cpp/particles/p2p_distr_mpi.hpp
@@ -300,7 +300,7 @@ public:
part_to_sort.back().second = idxPart;
}
- assert(part_to_sort.size() == (current_my_nb_particles_per_partition[idxPartition]));
+ assert(partsize_t(part_to_sort.size()) == (current_my_nb_particles_per_partition[idxPartition]));
std::sort(part_to_sort.begin(), part_to_sort.end(),
[](const std::pair<long int,partsize_t>& p1,
@@ -575,7 +575,7 @@ public:
const int nbNeighbors = my_tree.getNeighbors(current_cell_idx, neighbors, neighbors_indexes, shift, true);
// with other interval
- for(size_t idx_neighbor = 0 ; idx_neighbor < nbNeighbors ; ++idx_neighbor){
+ 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){
@@ -726,7 +726,7 @@ public:
for(size_t idx_1 = 0 ; idx_1 < intervals.size() ; ++idx_1){
// with other interval
- for(size_t idx_neighbor = 0 ; idx_neighbor < nbNeighbors ; ++idx_neighbor){
+ for(int idx_neighbor = 0 ; idx_neighbor < nbNeighbors ; ++idx_neighbor){
if(currenct_cell_idx < neighbors_indexes[idx_neighbor]){
cells_locker.lock(neighbors_indexes[idx_neighbor]);
diff --git a/bfps/cpp/particles/particles_distr_mpi.hpp b/bfps/cpp/particles/particles_distr_mpi.hpp
index 48559518..f0c09fd9 100644
--- a/bfps/cpp/particles/particles_distr_mpi.hpp
+++ b/bfps/cpp/particles/particles_distr_mpi.hpp
@@ -35,6 +35,9 @@ protected:
TAG_LOW_UP_MOVED_PARTICLES_INDEXES,
TAG_UP_LOW_MOVED_PARTICLES_INDEXES,
+ TAG_LOW_UP_MOVED_PARTICLES_DATA,
+ TAG_UP_LOW_MOVED_PARTICLES_DATA,
+
TAG_LOW_UP_MOVED_PARTICLES_RHS,
TAG_LOW_UP_MOVED_PARTICLES_RHS_MAX = TAG_LOW_UP_MOVED_PARTICLES_RHS+MaxNbRhs,
@@ -503,13 +506,14 @@ public:
////////////////////////////////////////////////////////////////////////////
- template <class computer_class, int size_particle_positions, int size_particle_rhs, int size_particle_index>
+ template <class computer_class, int size_particle_positions, int size_particle_data, int size_particle_rhs, int size_particle_index>
void redistribute(computer_class& in_computer,
partsize_t current_my_nb_particles_per_partition[],
partsize_t* nb_particles,
std::unique_ptr<real_number[]>* inout_positions_particles,
std::unique_ptr<real_number[]> inout_rhs_particles[], const int in_nb_rhs,
- std::unique_ptr<partsize_t[]>* inout_index_particles){
+ std::unique_ptr<partsize_t[]>* inout_index_particles,
+ std::unique_ptr<real_number[]>* inout_data_particles){
TIMEZONE("redistribute");
// Some latest processes might not be involved
@@ -537,7 +541,13 @@ public:
},
[&](const partsize_t idx1, const partsize_t idx2){
for(int idx_val = 0 ; idx_val < size_particle_index ; ++idx_val){
- std::swap((*inout_index_particles)[idx1], (*inout_index_particles)[idx2]);
+ std::swap((*inout_index_particles)[size_particle_index*idx1+idx_val],
+ (*inout_index_particles)[size_particle_index*idx2+idx_val]);
+ }
+
+ for(int idx_val = 0 ; idx_val < size_particle_data ; ++idx_val){
+ std::swap((*inout_data_particles)[size_particle_data*idx1+idx_val],
+ (*inout_data_particles)[size_particle_data*idx2+idx_val]);
}
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
@@ -562,7 +572,13 @@ public:
},
[&](const partsize_t idx1, const partsize_t idx2){
for(int idx_val = 0 ; idx_val < size_particle_index ; ++idx_val){
- std::swap((*inout_index_particles)[idx1], (*inout_index_particles)[idx2]);
+ std::swap((*inout_index_particles)[size_particle_index*idx1+idx_val],
+ (*inout_index_particles)[size_particle_index*idx2+idx_val]);
+ }
+
+ for(int idx_val = 0 ; idx_val < size_particle_data ; ++idx_val){
+ std::swap((*inout_data_particles)[size_particle_data*idx1+idx_val],
+ (*inout_data_particles)[size_particle_data*idx2+idx_val]);
}
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
@@ -581,6 +597,8 @@ public:
std::unique_ptr<real_number[]> newParticlesUp;
std::unique_ptr<partsize_t[]> newParticlesLowIndexes;
std::unique_ptr<partsize_t[]> newParticlesUpIndexes;
+ std::unique_ptr<real_number[]> newParticlesLowData;
+ std::unique_ptr<real_number[]> newParticlesUpData;
std::vector<std::unique_ptr<real_number[]>> newParticlesLowRhs(in_nb_rhs);
std::vector<std::unique_ptr<real_number[]>> newParticlesUpRhs(in_nb_rhs);
@@ -607,13 +625,24 @@ public:
assert(nbOutLower*size_particle_positions < std::numeric_limits<int>::max());
AssertMpi(MPI_Isend(&(*inout_positions_particles)[0], int(nbOutLower*size_particle_positions), particles_utils::GetMpiType(real_number()), (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_LOW_UP_MOVED_PARTICLES,
MPI_COMM_WORLD, &mpiRequests.back()));
+
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
- assert(nbOutLower < std::numeric_limits<int>::max());
- AssertMpi(MPI_Isend(&(*inout_index_particles)[0], int(nbOutLower), particles_utils::GetMpiType(partsize_t()),
+ assert(nbOutLower*size_particle_index < std::numeric_limits<int>::max());
+ AssertMpi(MPI_Isend(&(*inout_index_particles)[0], int(nbOutLower*size_particle_index), particles_utils::GetMpiType(partsize_t()),
(my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_LOW_UP_MOVED_PARTICLES_INDEXES,
MPI_COMM_WORLD, &mpiRequests.back()));
+ if(size_particle_data){
+ whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
+ mpiRequests.emplace_back();
+ assert(nbOutLower*size_particle_data < std::numeric_limits<int>::max());
+ AssertMpi(MPI_Isend(&(*inout_data_particles)[0], int(nbOutLower*size_particle_data),
+ particles_utils::GetMpiType(partsize_t()),
+ (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_LOW_UP_MOVED_PARTICLES_DATA,
+ MPI_COMM_WORLD, &mpiRequests.back()));
+ }
+
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
@@ -643,13 +672,23 @@ public:
AssertMpi(MPI_Isend(&(*inout_positions_particles)[(myTotalNbParticles-nbOutUpper)*size_particle_positions],
int(nbOutUpper*size_particle_positions), particles_utils::GetMpiType(real_number()), (my_rank+1)%nb_processes_involved, TAG_UP_LOW_MOVED_PARTICLES,
MPI_COMM_WORLD, &mpiRequests.back()));
+
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
- assert(nbOutUpper < std::numeric_limits<int>::max());
- AssertMpi(MPI_Isend(&(*inout_index_particles)[(myTotalNbParticles-nbOutUpper)], int(nbOutUpper),
+ assert(nbOutUpper*size_particle_index < std::numeric_limits<int>::max());
+ AssertMpi(MPI_Isend(&(*inout_index_particles)[(myTotalNbParticles-nbOutUpper)*size_particle_index], int(nbOutUpper*size_particle_index),
particles_utils::GetMpiType(partsize_t()), (my_rank+1)%nb_processes_involved, TAG_UP_LOW_MOVED_PARTICLES_INDEXES,
MPI_COMM_WORLD, &mpiRequests.back()));
+ if(size_particle_data){
+ whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
+ mpiRequests.emplace_back();
+ assert(nbOutUpper*size_particle_data < std::numeric_limits<int>::max());
+ AssertMpi(MPI_Isend(&(*inout_data_particles)[(myTotalNbParticles-nbOutUpper)*size_particle_data],
+ int(nbOutUpper*size_particle_data),
+ particles_utils::GetMpiType(partsize_t()), (my_rank+1)%nb_processes_involved, TAG_UP_LOW_MOVED_PARTICLES_DATA,
+ MPI_COMM_WORLD, &mpiRequests.back()));
+ }
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
@@ -684,14 +723,26 @@ public:
(my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_UP_LOW_MOVED_PARTICLES,
MPI_COMM_WORLD, &mpiRequests.back()));
- newParticlesLowIndexes.reset(new partsize_t[nbNewFromLow]);
+ newParticlesLowIndexes.reset(new partsize_t[nbNewFromLow*size_particle_index]);
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
- assert(nbNewFromLow < std::numeric_limits<int>::max());
- AssertMpi(MPI_Irecv(&newParticlesLowIndexes[0], int(nbNewFromLow), particles_utils::GetMpiType(partsize_t()),
+ assert(nbNewFromLow*size_particle_index < std::numeric_limits<int>::max());
+ AssertMpi(MPI_Irecv(&newParticlesLowIndexes[0], int(nbNewFromLow*size_particle_index),
+ particles_utils::GetMpiType(partsize_t()),
(my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_UP_LOW_MOVED_PARTICLES_INDEXES,
MPI_COMM_WORLD, &mpiRequests.back()));
+ if(size_particle_data){
+ newParticlesLowData.reset(new real_number[nbNewFromLow*size_particle_data]);
+ whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
+ mpiRequests.emplace_back();
+ assert(nbNewFromLow*size_particle_data < std::numeric_limits<int>::max());
+ AssertMpi(MPI_Irecv(&newParticlesLowData[0], int(nbNewFromLow*size_particle_data),
+ particles_utils::GetMpiType(real_number()),
+ (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_UP_LOW_MOVED_PARTICLES_DATA,
+ MPI_COMM_WORLD, &mpiRequests.back()));
+ }
+
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
newParticlesLowRhs[idx_rhs].reset(new real_number[nbNewFromLow*size_particle_rhs]);
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
@@ -713,14 +764,26 @@ public:
AssertMpi(MPI_Irecv(&newParticlesUp[0], int(nbNewFromUp*size_particle_positions), particles_utils::GetMpiType(real_number()), (my_rank+1)%nb_processes_involved, TAG_LOW_UP_MOVED_PARTICLES,
MPI_COMM_WORLD, &mpiRequests.back()));
- newParticlesUpIndexes.reset(new partsize_t[nbNewFromUp]);
+ newParticlesUpIndexes.reset(new partsize_t[nbNewFromUp*size_particle_index]);
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
- assert(nbNewFromUp < std::numeric_limits<int>::max());
- AssertMpi(MPI_Irecv(&newParticlesUpIndexes[0], int(nbNewFromUp), particles_utils::GetMpiType(partsize_t()),
+ assert(nbNewFromUp*size_particle_index < std::numeric_limits<int>::max());
+ AssertMpi(MPI_Irecv(&newParticlesUpIndexes[0], int(nbNewFromUp*size_particle_index),
+ particles_utils::GetMpiType(partsize_t()),
(my_rank+1)%nb_processes_involved, TAG_LOW_UP_MOVED_PARTICLES_INDEXES,
MPI_COMM_WORLD, &mpiRequests.back()));
+ if(size_particle_data){
+ newParticlesUpData.reset(new real_number[nbNewFromUp*size_particle_data]);
+ whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
+ mpiRequests.emplace_back();
+ assert(nbNewFromUp*size_particle_data < std::numeric_limits<int>::max());
+ AssertMpi(MPI_Irecv(&newParticlesUpData[0], int(nbNewFromUp*size_particle_data),
+ particles_utils::GetMpiType(real_number()),
+ (my_rank+1)%nb_processes_involved, TAG_LOW_UP_MOVED_PARTICLES_DATA,
+ MPI_COMM_WORLD, &mpiRequests.back()));
+ }
+
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
newParticlesUpRhs[idx_rhs].reset(new real_number[nbNewFromUp*size_particle_rhs]);
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
@@ -750,7 +813,8 @@ public:
const partsize_t myTotalNewNbParticles = nbOldParticlesInside + nbNewFromLow + nbNewFromUp;
std::unique_ptr<real_number[]> newArray(new real_number[myTotalNewNbParticles*size_particle_positions]);
- std::unique_ptr<partsize_t[]> newArrayIndexes(new partsize_t[myTotalNewNbParticles]);
+ std::unique_ptr<partsize_t[]> newArrayIndexes(new partsize_t[myTotalNewNbParticles*size_particle_index]);
+ std::unique_ptr<real_number[]> newArrayData(new real_number[myTotalNewNbParticles*size_particle_data]);
std::vector<std::unique_ptr<real_number[]>> newArrayRhs(in_nb_rhs);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
newArrayRhs[idx_rhs].reset(new real_number[myTotalNewNbParticles*size_particle_rhs]);
@@ -760,7 +824,8 @@ public:
if(nbNewFromLow){
const particles_utils::fixed_copy fcp(0, 0, nbNewFromLow);
fcp.copy(newArray, newParticlesLow, size_particle_positions);
- fcp.copy(newArrayIndexes, newParticlesLowIndexes);
+ fcp.copy(newArrayIndexes, newParticlesLowIndexes, size_particle_index);
+ fcp.copy(newArrayData, newParticlesLowData, size_particle_data);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
fcp.copy(newArrayRhs[idx_rhs], newParticlesLowRhs[idx_rhs], size_particle_rhs);
}
@@ -770,7 +835,8 @@ public:
{
const particles_utils::fixed_copy fcp(nbNewFromLow, nbOutLower, nbOldParticlesInside);
fcp.copy(newArray, (*inout_positions_particles), size_particle_positions);
- fcp.copy(newArrayIndexes, (*inout_index_particles));
+ fcp.copy(newArrayIndexes, (*inout_index_particles), size_particle_index);
+ fcp.copy(newArrayData, (*inout_data_particles), size_particle_data);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
fcp.copy(newArrayRhs[idx_rhs], inout_rhs_particles[idx_rhs], size_particle_rhs);
}
@@ -780,7 +846,8 @@ public:
if(nbNewFromUp){
const particles_utils::fixed_copy fcp(nbNewFromLow+nbOldParticlesInside, 0, nbNewFromUp);
fcp.copy(newArray, newParticlesUp, size_particle_positions);
- fcp.copy(newArrayIndexes, newParticlesUpIndexes);
+ fcp.copy(newArrayIndexes, newParticlesUpIndexes, size_particle_index);
+ fcp.copy(newArrayData, newParticlesUpData, size_particle_data);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
fcp.copy(newArrayRhs[idx_rhs], newParticlesUpRhs[idx_rhs], size_particle_rhs);
}
@@ -788,6 +855,7 @@ public:
(*inout_positions_particles) = std::move(newArray);
(*inout_index_particles) = std::move(newArrayIndexes);
+ (*inout_data_particles) = std::move(newArrayData);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
inout_rhs_particles[idx_rhs] = std::move(newArrayRhs[idx_rhs]);
}
@@ -808,7 +876,13 @@ public:
},
[&](const partsize_t idx1, const partsize_t idx2){
for(int idx_val = 0 ; idx_val < size_particle_index ; ++idx_val){
- std::swap((*inout_index_particles)[idx1], (*inout_index_particles)[idx2]);
+ std::swap((*inout_index_particles)[size_particle_index*idx1 + idx_val],
+ (*inout_index_particles)[size_particle_index*idx2 + idx_val]);
+ }
+
+ for(int idx_val = 0 ; idx_val < size_particle_data ; ++idx_val){
+ std::swap((*inout_data_particles)[size_particle_data*idx1 + idx_val],
+ (*inout_data_particles)[size_particle_data*idx2 + idx_val]);
}
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
diff --git a/bfps/cpp/particles/particles_system.hpp b/bfps/cpp/particles/particles_system.hpp
index 02767a8b..12bf6c29 100644
--- a/bfps/cpp/particles/particles_system.hpp
+++ b/bfps/cpp/particles/particles_system.hpp
@@ -12,8 +12,10 @@
#include "particles_adams_bashforth.hpp"
#include "scope_timer.hpp"
+#include "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_rhs>
+ int size_particle_rhs, class p2p_computer_class, int size_particle_data>
class particles_system : public abstract_particles_system<partsize_t, real_number> {
MPI_Comm mpi_com;
@@ -47,6 +49,10 @@ class particles_system : public abstract_particles_system<partsize_t, real_numbe
int step_idx;
+ p2p_distr_mpi<partsize_t, real_number> distr_p2p;
+ p2p_computer_class computer_p2p;
+ std::unique_ptr<real_number[]> my_particles_data;
+
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,
@@ -57,7 +63,8 @@ public:
const field_class& in_field,
MPI_Comm in_mpi_com,
const partsize_t in_total_nb_particles,
- const int in_current_iteration = 1)
+ const int in_current_iteration = 1,
+ const real_number in_cutoff = 1.)
: mpi_com(in_mpi_com),
current_partition_interval({in_local_field_offset[IDX_Z], in_local_field_offset[IDX_Z] + in_local_field_dims[IDX_Z]}),
partition_interval_size(current_partition_interval.second - current_partition_interval.first),
@@ -69,7 +76,8 @@ public:
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){
+ 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){
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]);
@@ -85,6 +93,8 @@ public:
my_particles_positions_indexes = particles_input.getMyParticlesIndexes();
my_particles_rhs = particles_input.getMyRhs();
my_nb_particles = particles_input.getLocalNbParticles();
+ // TODO P2P get it from loader
+ my_particles_data.reset(new real_number[my_nb_particles*size_particle_data]());
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*3+IDX_Z], IDX_Z);
@@ -107,6 +117,10 @@ public:
my_particles_rhs[idx_rhs][idx2*size_particle_rhs + idx_val]);
}
}
+ for(int idx_val = 0 ; idx_val < size_particle_data ; ++idx_val){
+ std::swap(my_particles_data[idx1*size_particle_data + idx_val],
+ my_particles_data[idx2*size_particle_data + idx_val]);
+ }
});
{// TODO remove
@@ -131,6 +145,17 @@ public:
interp_neighbours);
}
+ void compute_p2p() final {
+ // TODO P2P
+ if(p2p_computer_class::isEmpty() == false){
+ TIMEZONE("particles_system::compute_p2p");
+ distr_p2p.template compute_distr<p2p_computer_class, 3, size_particle_data, size_particle_rhs>(
+ computer_p2p, current_my_nb_particles_per_partition.get(),
+ my_particles_positions.get(), my_particles_data.get(), my_particles_rhs.front().get(),
+ my_particles_positions_indexes.get());
+ }
+ }
+
template <class sample_field_class, int sample_size_particle_rhs>
void sample_compute(const sample_field_class& sample_field,
real_number sample_rhs[]) {
@@ -179,13 +204,14 @@ public:
void redistribute() final {
TIMEZONE("particles_system::redistribute");
- particles_distr.template redistribute<computer_class, 3, size_particle_rhs, 1>(
+ particles_distr.template redistribute<computer_class, 3, size_particle_data, 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);
+ &my_particles_positions_indexes,
+ &my_particles_data);
}
void inc_step_idx() final {
@@ -210,6 +236,7 @@ public:
void completeLoop(const real_number dt) final {
TIMEZONE("particles_system::completeLoop");
compute();
+ compute_p2p();
move(dt);
redistribute();
inc_step_idx();
diff --git a/bfps/cpp/particles/particles_system_builder.hpp b/bfps/cpp/particles/particles_system_builder.hpp
index a3bc689d..f9ce512d 100644
--- a/bfps/cpp/particles/particles_system_builder.hpp
+++ b/bfps/cpp/particles/particles_system_builder.hpp
@@ -8,6 +8,7 @@
#include "particles_system.hpp"
#include "particles_input_hdf5.hpp"
#include "particles_generic_interp.hpp"
+#include "p2p_computer_empty.hpp"
#include "field.hpp"
#include "kspace.hpp"
@@ -109,7 +110,7 @@ inline RetType evaluate(IterType1 value1, IterType2 value2, Args... args){
///
//////////////////////////////////////////////////////////////////////////////
-template <class partsize_t, class field_rnumber, field_backend be, field_components fc, class particles_rnumber>
+template <class partsize_t, class field_rnumber, field_backend be, field_components fc, class particles_rnumber, class p2p_computer_class>
struct particles_system_build_container {
template <const int interpolation_size, const int spline_mode>
static std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>> instanciate(
@@ -196,7 +197,8 @@ struct particles_system_build_container {
using particles_system_type = particles_system<partsize_t, particles_rnumber, field_rnumber,
field<field_rnumber, be, fc>,
particles_generic_interp<particles_rnumber, interpolation_size,spline_mode>,
- interpolation_size, ncomp(fc)>;
+ interpolation_size, ncomp(fc),
+ p2p_computer_class, p2p_computer_class::size_data>;
particles_system_type* part_sys = new particles_system_type(field_grid_dim,
spatial_box_width,
spatial_box_offset,
@@ -210,6 +212,7 @@ struct particles_system_build_container {
nparticles,
in_current_iteration);
+ // TODO P2P load particle data too
// Load particles from hdf5
particles_input_hdf5<partsize_t, particles_rnumber, 3,3> generator(mpi_comm, fname_input,
inDatanameState, inDatanameRhs, my_spatial_low_limit_z, my_spatial_up_limit_z);
@@ -251,7 +254,28 @@ inline std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>
return Template_double_for_if::evaluate<std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>,
int, 1, 11, 1, // interpolation_size
int, 0, 3, 1, // spline_mode
- particles_system_build_container<partsize_t, field_rnumber,be,fc,particles_rnumber>>(
+ particles_system_build_container<partsize_t, field_rnumber,be,fc,particles_rnumber, p2p_computer_empty<particles_rnumber,partsize_t>>>(
+ interpolation_size, // template iterator 1
+ spline_mode, // template iterator 2
+ fs_field,fs_kk, nsteps, nparticles, fname_input, inDatanameState, inDatanameRhs, mpi_comm, in_current_iteration);
+}
+
+template <class partsize_t, class field_rnumber, field_backend be, field_components fc, class p2p_computer_class, class particles_rnumber = double>
+inline std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>> particles_system_builder_with_p2p(
+ const field<field_rnumber, be, fc>* fs_field, // (field object)
+ const kspace<be, SMOOTH>* fs_kk, // (kspace object, contains dkx, dky, dkz)
+ const int nsteps, // to check coherency between parameters and hdf input file (nb rhs)
+ const partsize_t nparticles, // to check coherency between parameters and hdf input file
+ const std::string& fname_input, // particles input filename
+ const std::string& inDatanameState, const std::string& inDatanameRhs, // input dataset names
+ const int interpolation_size,
+ const int spline_mode,
+ MPI_Comm mpi_comm,
+ const int in_current_iteration){
+ return Template_double_for_if::evaluate<std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>,
+ int, 1, 11, 1, // interpolation_size
+ int, 0, 3, 1, // spline_mode
+ particles_system_build_container<partsize_t, field_rnumber,be,fc,particles_rnumber,p2p_computer_class>>(
interpolation_size, // template iterator 1
spline_mode, // template iterator 2
fs_field,fs_kk, nsteps, nparticles, fname_input, inDatanameState, inDatanameRhs, mpi_comm, in_current_iteration);
diff --git a/setup.py b/setup.py
index 9bba1701..f257fd7f 100644
--- a/setup.py
+++ b/setup.py
@@ -127,7 +127,8 @@ src_file_list = ['full_code/joint_acc_vel_stats',
'spline_n10',
'Lagrange_polys',
'scope_timer',
- 'full_code/NSVEparticles']
+ 'full_code/NSVEparticles',
+ 'full_code/NSVEparticlesP2P']
particle_headers = [
'cpp/particles/particles_distr_mpi.hpp',
--
GitLab