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