particles_system.hpp 20.36 KiB
/******************************************************************************
* *
* 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 PARTICLES_SYSTEM_HPP
#define PARTICLES_SYSTEM_HPP
#include <array>
#include "abstract_particles_system.hpp"
#include "abstract_particles_system_with_p2p.hpp"
#include "particles_distr_mpi.hpp"
#include "particles_output_hdf5.hpp"
#include "particles_output_mpiio.hpp"
#include "particles_field_computer.hpp"
#include "abstract_particles_input.hpp"
#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_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;
const 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 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.front().get(),
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 {
TIMEZONE("particles_system::completeLoop");
compute();
compute_p2p();
compute_particles_inner();
move(dt);
enforce_unit_orientation();
redistribute();
inc_step_idx();
shift_rhs_vectors();
}
void complete2ndOrderLoop(const real_number dt) final {
assert(size_particle_positions == 6);
assert(size_particle_rhs == 6);
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();
}
void completeLoopWithVorticity(
const real_number dt,
const real_number particle_extra_field[]) final {
TIMEZONE("particles_system::completeLoopWithVorticity");
compute();
compute_p2p();
compute_sphere_particles_inner(particle_extra_field);
move(dt);
enforce_unit_orientation();
redistribute();
inc_step_idx();
shift_rhs_vectors();
}
void completeLoopWithVelocityGradient(
const real_number dt,
const real_number particle_extra_field[]) final {
TIMEZONE("particles_system::completeLoopWithVelocityGradient");
compute();
compute_p2p();
compute_ellipsoid_particles_inner(particle_extra_field);
move(dt);
enforce_unit_orientation();
redistribute();
inc_step_idx();
shift_rhs_vectors();
}
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;
}
};
#endif