diff --git a/bfps/cpp/particles/NavierStokes-v2.0.1-single-particles.cpp b/bfps/cpp/particles/NavierStokes-v2.0.1-single-particles.cpp
index b6cd2a11d3a6521764538274faec23b31a65b53b..8bf93750d8676276f387c4caf4c0e28ad54fe4db 100644
--- a/bfps/cpp/particles/NavierStokes-v2.0.1-single-particles.cpp
+++ b/bfps/cpp/particles/NavierStokes-v2.0.1-single-particles.cpp
@@ -545,7 +545,7 @@ int main(int argc, char *argv[])
// ps0->read();
// Interpolator size
- const int InterpNbNeighbors = 5;
+ const int InterpNbNeighbors = 1; // start with one and then use 2 or 3
// Only the Z grid
const std::array<size_t,3> field_grid_dim{nx,ny,nz};
@@ -553,7 +553,7 @@ int main(int argc, char *argv[])
const int myPartitionInterval[2] = {fs->cvelocity->rlayout->start[0],
fs->cvelocity->rlayout->start[0] + fs->cvelocity->rlayout->subsizes[0]};
- const std::array<double,3> spatial_box_width{dkx, dky, dkz};
+ const std::array<double,3> spatial_box_width{2PI/dkx, 2PI/dky, 2PI/dkz};
const double spatial_partition_width = spatial_box_width/double(field_grid_dim);
const double my_spatial_low_limit = myPartitionInterval[0]*spatial_partition_width;
const double my_spatial_up_limit = myPartitionInterval[1]*spatial_partition_width;
diff --git a/bfps/cpp/particles/abstract_particle_system.hpp b/bfps/cpp/particles/abstract_particle_system.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7ffee82f17cd1200eba1b680be9f6bb3074b8bf8
--- /dev/null
+++ b/bfps/cpp/particles/abstract_particle_system.hpp
@@ -0,0 +1,27 @@
+#ifndef ABSTRACT_PARTICLE_SYSTEM_HPP
+#define ABSTRACT_PARTICLE_SYSTEM_HPP
+
+class abstract_particle_system {
+public:
+ virtual void compute() = 0;
+
+ virtual void move(const double dt) = 0;
+
+ virtual void redistribute() = 0;
+
+ virtual void inc_step_idx() = 0;
+
+ virtual void shift_rhs_vectors() = 0;
+
+ virtual void completeLoop(const double dt) = 0;
+
+ virtual const double* getParticlesPositions() const = 0;
+
+ virtual const double* getParticlesCurrentRhs() const = 0;
+
+ virtual const int* getParticlesIndexes() const = 0;
+
+ virtual int getLocalNbParticles() const = 0;
+};
+
+#endif
diff --git a/bfps/cpp/particles/abstract_particles_distr.hpp b/bfps/cpp/particles/abstract_particles_distr.hpp
index 1e6508cb359bcfa8fb7dd7a5bbbacd71f3498673..0cf7147c23c842b569b1dca7c66ae1ddf46db82f 100644
--- a/bfps/cpp/particles/abstract_particles_distr.hpp
+++ b/bfps/cpp/particles/abstract_particles_distr.hpp
@@ -431,7 +431,7 @@ public:
// Find particles outside my interval
const int nbOutLower = particles_utils::partition_extra<size_particle_positions>(&(*inout_positions_particles)[0], current_my_nb_particles_per_partition[0],
[&](const double val[]){
- const bool isLower = val[2] < mySpatialLowLimit;
+ const bool isLower = val[IDX_Z] < mySpatialLowLimit;
return isLower;
},
[&](const int idx1, const int idx2){
@@ -454,7 +454,7 @@ public:
&(*inout_positions_particles)[(current_offset_particles_for_partition[current_partition_size-1]+offesetOutLow)*size_particle_positions],
myTotalNbParticles - (current_offset_particles_for_partition[current_partition_size-1]+offesetOutLow),
[&](const double val[]){
- const bool isUpper = mySpatialUpLimit <= val[2];
+ const bool isUpper = mySpatialUpLimit <= val[IDX_Z];
return !isUpper;
},
[&](const int idx1, const int idx2){
@@ -578,7 +578,7 @@ public:
MPI_COMM_WORLD, &mpiRequests.back()));
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
- newParticlesLowRhs[idx_rhs].reset(new double[nbNewFromLow]);
+ newParticlesLowRhs[idx_rhs].reset(new double[nbNewFromLow*size_particle_rhs]);
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
AssertMpi(MPI_Irecv(&newParticlesLowRhs[idx_rhs][0], nbNewFromLow*size_particle_rhs, MPI_DOUBLE, (my_rank-1+nb_processes)%nb_processes, TAG_UP_LOW_MOVED_PARTICLES_RHS+idx_rhs,
@@ -605,7 +605,7 @@ public:
MPI_COMM_WORLD, &mpiRequests.back()));
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
- newParticlesUpRhs[idx_rhs].reset(new double[nbNewFromUp]);
+ newParticlesUpRhs[idx_rhs].reset(new double[nbNewFromUp*size_particle_rhs]);
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
AssertMpi(MPI_Irecv(&newParticlesUpRhs[idx_rhs][0], nbNewFromUp*size_particle_rhs, MPI_DOUBLE, (my_rank+1)%nb_processes, TAG_LOW_UP_MOVED_PARTICLES_RHS+idx_rhs,
@@ -650,31 +650,31 @@ public:
if(myTotalNewNbParticles > myTotalNbParticlesAllocated){
DEBUG_MSG("[%d] reuse array\n", my_rank);
std::unique_ptr<double[]> newArray(new double[myTotalNewNbParticles*size_particle_positions]);
- std::unique_ptr<int[]> newArrayIndexes(new int[myTotalNewNbParticles*size_particle_positions]);
+ std::unique_ptr<int[]> newArrayIndexes(new int[myTotalNewNbParticles]);
std::vector<std::unique_ptr<double[]>> newArrayRhs(in_nb_rhs);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
- newArrayRhs[idx_rhs].reset(new double[myTotalNewNbParticles*size_particle_positions]);
+ newArrayRhs[idx_rhs].reset(new double[myTotalNewNbParticles*size_particle_rhs]);
}
if(nbNewFromLow){
memcpy(&newArray[0], &newParticlesLow[0], sizeof(double)*nbNewFromLow*size_particle_positions);
memcpy(&newArrayIndexes[0], &newParticlesLowIndexes[0], sizeof(int)*nbNewFromLow);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
- memcpy(&newArrayRhs[idx_rhs][0], &newParticlesLowRhs[idx_rhs][0], sizeof(double)*nbNewFromLow);
+ memcpy(&newArrayRhs[idx_rhs][0], &newParticlesLowRhs[idx_rhs][0], sizeof(double)*nbNewFromLow*size_particle_rhs);
}
}
memcpy(&newArray[nbNewFromLow*size_particle_positions], &(*inout_positions_particles)[nbOutLower*size_particle_positions], sizeof(double)*nbOldParticlesInside*size_particle_positions);
memcpy(&newArrayIndexes[nbNewFromLow], &(*inout_positions_particles)[nbOutLower], sizeof(int)*nbOldParticlesInside);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
- memcpy(&newArrayRhs[idx_rhs][nbNewFromLow], &inout_positions_particles[idx_rhs][nbOutLower], sizeof(double)*nbOldParticlesInside);
+ memcpy(&newArrayRhs[idx_rhs][nbNewFromLow*size_particle_rhs], &inout_positions_particles[idx_rhs][nbOutLower*size_particle_rhs], sizeof(double)*nbOldParticlesInside*size_particle_rhs);
}
if(nbNewFromUp){
memcpy(&newArray[(nbNewFromLow+nbOldParticlesInside)*size_particle_positions], &newParticlesUp[0], sizeof(double)*nbNewFromUp*size_particle_positions);
memcpy(&newArrayIndexes[(nbNewFromLow+nbOldParticlesInside)], &newParticlesUpIndexes[0], sizeof(int)*nbNewFromUp);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
- memcpy(&newArrayRhs[idx_rhs][(nbNewFromLow+nbOldParticlesInside)], &newParticlesUpRhs[idx_rhs][0], sizeof(double)*nbNewFromUp);
+ memcpy(&newArrayRhs[idx_rhs][(nbNewFromLow+nbOldParticlesInside)*size_particle_rhs], &newParticlesUpRhs[idx_rhs][0], sizeof(double)*nbNewFromUp*size_particle_rhs);
}
}
@@ -695,21 +695,21 @@ public:
memcpy(&(*inout_positions_particles)[0], &newParticlesLow[0], sizeof(double)*nbOutLower*size_particle_positions);
memcpy(&(*inout_index_particles)[0], &newParticlesLowIndexes[0], sizeof(int)*nbOutLower);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
- memcpy(&inout_rhs_particles[idx_rhs][0], &newParticlesLowRhs[idx_rhs][0], sizeof(double)*nbOutLower);
+ memcpy(&inout_rhs_particles[idx_rhs][0], &newParticlesLowRhs[idx_rhs][0], sizeof(double)*nbOutLower*size_particle_rhs);
}
}
if(nbNewFromLow){
memcpy(&(*inout_positions_particles)[(nbOutLower+nbOldParticlesInside)*size_particle_positions], &newParticlesLow[nbOutLower*size_particle_positions], sizeof(double)*nbLowToMoveBack*size_particle_positions);
memcpy(&(*inout_index_particles)[(nbOutLower+nbOldParticlesInside)], &newParticlesLowIndexes[nbOutLower], sizeof(int)*nbLowToMoveBack);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
- memcpy(&inout_rhs_particles[idx_rhs][(nbOutLower+nbOldParticlesInside)], &newParticlesLowRhs[idx_rhs][nbOutLower], sizeof(double)*nbLowToMoveBack);
+ memcpy(&inout_rhs_particles[idx_rhs][(nbOutLower+nbOldParticlesInside)*size_particle_rhs], &newParticlesLowRhs[idx_rhs][nbOutLower*size_particle_rhs], sizeof(double)*nbLowToMoveBack*size_particle_rhs);
}
}
if(nbNewFromUp){
memcpy(&(*inout_positions_particles)[(nbNewFromLow+nbOldParticlesInside)*size_particle_positions], &newParticlesUp[0], sizeof(double)*nbNewFromUp*size_particle_positions);
memcpy(&(*inout_index_particles)[(nbNewFromLow+nbOldParticlesInside)], &newParticlesUpIndexes[0], sizeof(int)*nbNewFromUp);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
- memcpy(&inout_rhs_particles[idx_rhs][(nbNewFromLow+nbOldParticlesInside)], &newParticlesUpRhs[0], sizeof(double)*nbNewFromUp);
+ memcpy(&inout_rhs_particles[idx_rhs][(nbNewFromLow+nbOldParticlesInside)*size_particle_rhs], &newParticlesUpRhs[0], sizeof(double)*nbNewFromUp*size_particle_rhs);
}
}
}
@@ -721,24 +721,26 @@ public:
memcpy(&(*inout_positions_particles)[0], &newParticlesLow[0], sizeof(double)*nbNewFromLow*size_particle_positions);
memcpy(&(*inout_index_particles)[0], &newParticlesLowIndexes[0], sizeof(int)*nbNewFromLow);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
- memcpy(&inout_rhs_particles[idx_rhs][0], &newParticlesLowRhs[idx_rhs][0], sizeof(double)*nbNewFromLow);
+ memcpy(&inout_rhs_particles[idx_rhs][0], &newParticlesLowRhs[idx_rhs][0], sizeof(double)*nbNewFromLow*size_particle_rhs);
}
}
if(nbNewFromUp){
memcpy(&(*inout_positions_particles)[nbNewFromLow*size_particle_positions], &newParticlesUp[0], sizeof(double)*nbUpToMoveBegin*size_particle_positions);
memcpy(&(*inout_index_particles)[nbNewFromLow], &newParticlesUpIndexes[0], sizeof(int)*nbUpToMoveBegin);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
- memcpy(&inout_rhs_particles[idx_rhs][0], &newParticlesLowRhs[idx_rhs][0], sizeof(double)*nbNewFromLow);
+ memcpy(&inout_rhs_particles[idx_rhs][nbNewFromLow*size_particle_rhs], &newParticlesLowRhs[idx_rhs][0], sizeof(double)*nbUpToMoveBegin*size_particle_rhs);
}
}
if(nbNewFromUp){
- memcpy(&(*inout_positions_particles)[(nbOutLower+nbOldParticlesInside)*size_particle_positions], &newParticlesUp[nbUpToMoveBegin*size_particle_positions],
+ memcpy(&(*inout_positions_particles)[(nbOutLower+nbOldParticlesInside)*size_particle_positions],
+ &newParticlesUp[nbUpToMoveBegin*size_particle_positions],
sizeof(double)*(nbNewFromUp-nbUpToMoveBegin)*size_particle_positions);
memcpy(&(*inout_index_particles)[(nbOutLower+nbOldParticlesInside)], &newParticlesUpIndexes[nbUpToMoveBegin],
sizeof(int)*(nbNewFromUp-nbUpToMoveBegin));
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
- memcpy(&inout_rhs_particles[idx_rhs][(nbOutLower+nbOldParticlesInside)], &newParticlesUpRhs[idx_rhs][nbUpToMoveBegin],
- sizeof(double)*(nbNewFromUp-nbUpToMoveBegin));
+ memcpy(&inout_rhs_particles[idx_rhs][(nbOutLower+nbOldParticlesInside)*size_particle_rhs],
+ &newParticlesUpRhs[idx_rhs][nbUpToMoveBegin*size_particle_rhs],
+ sizeof(double)*(nbNewFromUp-nbUpToMoveBegin)*size_particle_rhs);
}
}
}
@@ -748,14 +750,14 @@ public:
memcpy(&(*inout_positions_particles)[0], &newParticlesLow[0], sizeof(double)*nbNewFromLow*size_particle_positions);
memcpy(&(*inout_index_particles)[0], &newParticlesLowIndexes[0], sizeof(int)*nbNewFromLow);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
- memcpy(&inout_rhs_particles[idx_rhs][0], &newParticlesLowRhs[idx_rhs][0], sizeof(double)*nbNewFromLow);
+ memcpy(&inout_rhs_particles[idx_rhs][0], &newParticlesLowRhs[idx_rhs][0], sizeof(double)*nbNewFromLow*size_particle_rhs);
}
}
if(nbNewFromUp){
memcpy(&(*inout_positions_particles)[0], &newParticlesUp[0], sizeof(double)*nbNewFromUp*size_particle_positions);
memcpy(&(*inout_index_particles)[0], &newParticlesUpIndexes[0], sizeof(int)*nbNewFromUp);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
- memcpy(&inout_rhs_particles[idx_rhs][0], &newParticlesUpRhs[idx_rhs][0], sizeof(double)*nbNewFromUp);
+ memcpy(&inout_rhs_particles[idx_rhs][0], &newParticlesUpRhs[idx_rhs][0], sizeof(double)*nbNewFromUp*size_particle_rhs);
}
}
const int padding = nbOutLower - nbNewFromLow+nbNewFromUp;
@@ -766,9 +768,9 @@ public:
&(*inout_index_particles)[(nbOutLower+nbOldParticlesInside-padding)],
sizeof(int)*padding);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
- memcpy(&inout_rhs_particles[idx_rhs][(nbNewFromLow+nbNewFromUp)],
- &inout_rhs_particles[idx_rhs][(nbOutLower+nbOldParticlesInside-padding)],
- sizeof(double)*padding);
+ memcpy(&inout_rhs_particles[idx_rhs][(nbNewFromLow+nbNewFromUp)*size_particle_rhs],
+ &inout_rhs_particles[idx_rhs][(nbOutLower+nbOldParticlesInside-padding)*size_particle_rhs],
+ sizeof(double)*padding*size_particle_rhs);
}
}
}
@@ -808,10 +810,10 @@ public:
current_offset_particles_for_partition[idxPartition+1] - current_offset_particles_for_partition[idxPartition]);
const double limitPartition = (idxPartition+1)*spatialPartitionWidth + mySpatialLowLimit;
for(int idx = 0 ; idx < current_offset_particles_for_partition[idxPartition+1] ; ++idx){
- assert((*inout_positions_particles)[idx*3+2] < limitPartition);
+ assert((*inout_positions_particles)[idx*3+IDX_Z] < limitPartition);
}
for(int idx = current_offset_particles_for_partition[idxPartition+1] ; idx < myTotalNbParticles ; ++idx){
- assert((*inout_positions_particles)[idx*3+2] >= limitPartition);
+ assert((*inout_positions_particles)[idx*3+IDX_Z] >= limitPartition);
}
}
}
diff --git a/bfps/cpp/particles/field_accessor.hpp b/bfps/cpp/particles/field_accessor.hpp
index e8937d31d79c6d7fc2ed65de3cf2d93b4468d7e1..cdf6c8560f86076548dc7b934a3f2440f4726a44 100644
--- a/bfps/cpp/particles/field_accessor.hpp
+++ b/bfps/cpp/particles/field_accessor.hpp
@@ -4,19 +4,23 @@
#include <algorithm>
#include <array>
+#include "particles_utils.hpp"
class field_accessor {
static const int nb_dim = 3;
const double* field_date;
std::array<size_t,3> local_field_dims;
- std::array<size_t,3> dim_offset;
+ std::array<size_t,3> local_field_offset;
+ std::array<size_t,3> field_memory_dims;
public:
field_accessor(const double* in_field_date, const std::array<size_t,3>& in_dims,
- const std::array<size_t,3>& in_dim_offset)
+ const std::array<size_t,3>& in_local_field_offset,
+ const std::array<size_t,3>& in_field_memory_dims)
: field_date(in_field_date), local_field_dims(in_dims),
- dim_offset(in_dim_offset){
+ local_field_offset(in_local_field_offset),
+ field_memory_dims(in_field_memory_dims){
}
~field_accessor(){}
@@ -27,17 +31,20 @@ public:
}
size_t getIndexFromGlobalPosition(const size_t in_global_x, const size_t in_global_y, const size_t in_global_z) const {
- return getIndexFromLocalPosition(in_global_x - dim_offset[0],
- in_global_y - dim_offset[1],
- in_global_z - dim_offset[2]);
+ return getIndexFromLocalPosition(in_global_x - local_field_offset[IDX_X],
+ in_global_y - local_field_offset[IDX_Y],
+ in_global_z - local_field_offset[IDX_Z]);
}
size_t getIndexFromLocalPosition(const size_t in_local_x, const size_t in_local_y, const size_t in_local_z) const {
- assert(in_local_x < local_field_dims[0]);
- assert(in_local_y < local_field_dims[1]);
- assert(in_local_z < local_field_dims[2]);
- return (((in_local_z)*local_field_dims[1] +
- in_local_y)*(local_field_dims[2]) + // TODO there was a +2 on local_field_dims[2]
+ assert(0 <= in_local_x && in_local_x < local_field_dims[IDX_X]);
+ assert(0 <= in_local_y && in_local_y < local_field_dims[IDX_Y]);
+ assert(0 <= in_local_z && in_local_z < local_field_dims[IDX_Z]);
+ static_assert(IDX_X == 2 && IDX_Y == 1 && IDX_Z == 0,
+ "Dimension idx does not match, please ensure getIndexFromLocalPosition"
+ "is correct before commenting this assert");
+ return (((in_local_z)*field_memory_dims[1] +
+ in_local_y)*(field_memory_dims[2]) +
in_local_x);
}
};
diff --git a/bfps/cpp/particles/main_tocompile.cpp b/bfps/cpp/particles/main_tocompile.cpp
index 38db2ec7bd4d9f1cf0d0fcc8b8cacfdd011652a6..9f6420d79b78e70a85232c77d693c616935c475f 100644
--- a/bfps/cpp/particles/main_tocompile.cpp
+++ b/bfps/cpp/particles/main_tocompile.cpp
@@ -9,6 +9,7 @@
#include "particles_interp_spline.hpp"
#include "abstract_particles_input.hpp"
#include "particles_input_hdf5.hpp"
+#include "particles_utils.hpp"
class random_particles : public abstract_particles_input {
const int nb_particles;
@@ -43,9 +44,9 @@ public:
std::unique_ptr<double[]> particles(new double[nb_particles*3]);
for(int idx_part = 0 ; idx_part < nb_particles ; ++idx_part){
- particles[idx_part*3+0] = drand48() * box_width;
- particles[idx_part*3+1] = drand48() * box_width;
- particles[idx_part*3+2] = (drand48() * (upper_limit-lower_limit))
+ particles[idx_part*3+IDX_X] = drand48() * box_width;
+ particles[idx_part*3+IDX_Y] = drand48() * box_width;
+ particles[idx_part*3+IDX_Z] = (drand48() * (upper_limit-lower_limit))
+ lower_limit;
}
@@ -91,27 +92,33 @@ int main(int argc, char** argv){
assert(my_rank < field_grid_dim[1]);
assert(my_rank < field_grid_dim[2]);
- const double partitionIntervalSize = double(field_grid_dim[2])/double(nb_processes);
- const int myPartitionInterval[2] = { int(partitionIntervalSize*my_rank), (my_rank==nb_processes-1?field_grid_dim[2]:int(partitionIntervalSize*(my_rank+1)))};
+ const double partitionIntervalSize = double(field_grid_dim[IDX_Z])/double(nb_processes);
+ const int myPartitionInterval[2] = { int(partitionIntervalSize*my_rank), (my_rank==nb_processes-1?field_grid_dim[IDX_Z]:int(partitionIntervalSize*(my_rank+1)))};
const std::array<double,3> spatial_box_width{10., 10., 10.};
- const double spatial_partition_width = spatial_box_width[2]/double(field_grid_dim[2]);
+ const double spatial_partition_width = spatial_box_width[IDX_Z]/double(field_grid_dim[IDX_Z]);
const double my_spatial_low_limit = myPartitionInterval[0]*spatial_partition_width;
const double my_spatial_up_limit = myPartitionInterval[1]*spatial_partition_width;
if(my_rank == 0){
- std::cout << "spatial_box_width = " << spatial_box_width[0] << " " << spatial_box_width[1] << " " << spatial_box_width[2] << std::endl;
+ std::cout << "spatial_box_width = " << spatial_box_width[IDX_X] << " " << spatial_box_width[IDX_Y] << " " << spatial_box_width[IDX_Z] << std::endl;
std::cout << "spatial_partition_width = " << spatial_partition_width << std::endl;
std::cout << "my_spatial_low_limit = " << my_spatial_low_limit << std::endl;
std::cout << "my_spatial_up_limit = " << my_spatial_up_limit << std::endl;
}
- std::array<size_t,3> local_field_dims{ field_grid_dim[0], field_grid_dim[1], myPartitionInterval[1]-myPartitionInterval[0]};
- std::array<size_t,3> local_field_offset{ 0, 0, myPartitionInterval[0]};
+ std::array<size_t,3> local_field_dims;
+ local_field_dims[IDX_X] = field_grid_dim[IDX_X];
+ local_field_dims[IDX_Y] = field_grid_dim[IDX_Y];
+ local_field_dims[IDX_Z] = myPartitionInterval[1]-myPartitionInterval[0];
+ std::array<size_t,3> local_field_offset;
+ local_field_offset[IDX_X] = 0;
+ local_field_offset[IDX_Y] = 0;
+ local_field_offset[IDX_Z] = myPartitionInterval[0];
- std::unique_ptr<double[]> field_data(new double[local_field_dims[0]*local_field_dims[1]*local_field_dims[2]]);
- particles_utils::memzero(field_data.get(), local_field_dims[0]*local_field_dims[1]*local_field_dims[2]);
+ std::unique_ptr<double[]> field_data(new double[local_field_dims[IDX_X]*local_field_dims[IDX_Y]*local_field_dims[IDX_Z]*3]);
+ particles_utils::memzero(field_data.get(), local_field_dims[IDX_X]*local_field_dims[IDX_Y]*local_field_dims[IDX_Z]*3);
particles_system<particles_interp_spline<InterpNbNeighbors,0>, InterpNbNeighbors> part_sys(field_grid_dim,
@@ -122,6 +129,7 @@ int main(int argc, char** argv){
field_data.get(),
local_field_dims,
local_field_offset,
+ local_field_dims,
MPI_COMM_WORLD);
int total_nb_particles;
@@ -134,7 +142,7 @@ int main(int argc, char** argv){
spatial_interval_per_proc[idx_proc] = partitionIntervalSize*spatial_partition_width*idx_proc;
std::cout << "spatial_interval_per_proc[idx_proc] " << spatial_interval_per_proc[idx_proc] << std::endl;
}
- spatial_interval_per_proc[nb_processes] = spatial_box_width[2];
+ spatial_interval_per_proc[nb_processes] = spatial_box_width[IDX_Z];
assert(my_spatial_low_limit == spatial_interval_per_proc[my_rank] || fabs((spatial_interval_per_proc[my_rank]-my_spatial_low_limit)/my_spatial_low_limit) < 1e-13);
assert(my_spatial_up_limit == spatial_interval_per_proc[my_rank+1] || fabs((spatial_interval_per_proc[my_rank+1]-my_spatial_up_limit)/my_spatial_up_limit) < 1e-13);
diff --git a/bfps/cpp/particles/particles_field_computer.hpp b/bfps/cpp/particles/particles_field_computer.hpp
index 1baac156b20005fb6e6c91a1bc0523bfd31c9550..0ca83ce40246c578a2d2afdf4a06565e55af1d78 100644
--- a/bfps/cpp/particles/particles_field_computer.hpp
+++ b/bfps/cpp/particles/particles_field_computer.hpp
@@ -6,6 +6,7 @@
#include "abstract_particles_distr.hpp"
#include "scope_timer.hpp"
+#include "particles_utils.hpp"
template <class interpolator_class, class field_class, int interp_neighbours, class positions_updater_class >
class particles_field_computer : public abstract_particles_distr<3,3,1> {
@@ -40,17 +41,17 @@ class particles_field_computer : public abstract_particles_distr<3,3,1> {
TIMEZONE("particles_field_computer::apply_computation");
for(int idxPart = 0 ; idxPart < nb_particles ; ++idxPart){
double bx[interp_neighbours*2+2], by[interp_neighbours*2+2], bz[interp_neighbours*2+2];
- interpolator.compute_beta(deriv[0], particles_positions[idxPart*3], bx);
- interpolator.compute_beta(deriv[1], particles_positions[idxPart*3+1], by);
- interpolator.compute_beta(deriv[2], particles_positions[idxPart*3+2], bz);
+ interpolator.compute_beta(deriv[IDX_X], particles_positions[idxPart*3+IDX_X], bx);
+ interpolator.compute_beta(deriv[IDX_Y], particles_positions[idxPart*3+IDX_Y], by);
+ interpolator.compute_beta(deriv[IDX_Z], particles_positions[idxPart*3+IDX_Z], bz);
- const int partGridIdx_x = int(particles_positions[idxPart*3]/box_step_width);
- const int partGridIdx_y = int(particles_positions[idxPart*3+1]/box_step_width);
- const int partGridIdx_z = int(particles_positions[idxPart*3+2]/box_step_width);
+ const int partGridIdx_x = int(particles_positions[idxPart*3+IDX_X]/box_step_width);
+ const int partGridIdx_y = int(particles_positions[idxPart*3+IDX_Y]/box_step_width);
+ const int partGridIdx_z = int(particles_positions[idxPart*3+IDX_Z]/box_step_width);
- assert(0 <= partGridIdx_z && partGridIdx_z < field_grid_dim[2]);
- assert(0 <= partGridIdx_x && partGridIdx_x < field_grid_dim[0]);
- assert(0 <= partGridIdx_y && partGridIdx_y < field_grid_dim[1]);
+ assert(0 <= partGridIdx_x && partGridIdx_x < field_grid_dim[IDX_X]);
+ assert(0 <= partGridIdx_y && partGridIdx_y < field_grid_dim[IDX_Y]);
+ assert(0 <= partGridIdx_z && partGridIdx_z < field_grid_dim[IDX_Z]);
const int interp_limit_mx = partGridIdx_x-interp_neighbours;
const int interp_limit_x = partGridIdx_x+interp_neighbours+1;
@@ -62,8 +63,8 @@ class particles_field_computer : public abstract_particles_distr<3,3,1> {
int nb_z_intervals;
if((partGridIdx_z-interp_neighbours) < 0){
- assert(partGridIdx_z+interp_neighbours+1 < field_grid_dim[2]);
- interp_limit_mz[0] = ((partGridIdx_z-interp_neighbours)+field_grid_dim[2])%field_grid_dim[2];
+ assert(partGridIdx_z+interp_neighbours+1 < field_grid_dim[IDX_Z]);
+ interp_limit_mz[0] = ((partGridIdx_z-interp_neighbours)+field_grid_dim[IDX_Z])%field_grid_dim[IDX_Z];
interp_limit_z[0] = current_partition_interval.second-1;
interp_limit_mz[1] = std::max(0, current_partition_interval.first);// max is not really needed here
@@ -73,10 +74,10 @@ class particles_field_computer : public abstract_particles_distr<3,3,1> {
}
else if(field_grid_dim[2] <= (partGridIdx_z+interp_neighbours+1)){
interp_limit_mz[0] = std::max(current_partition_interval.first, partGridIdx_z-interp_neighbours);
- interp_limit_z[0] = std::min(int(field_grid_dim[2])-1,current_partition_interval.second-1);// max is not really needed here
+ interp_limit_z[0] = std::min(int(field_grid_dim[IDX_Z])-1,current_partition_interval.second-1);// max is not really needed here
interp_limit_mz[1] = std::max(0, current_partition_interval.first);
- interp_limit_z[1] = std::min(int((partGridIdx_z+interp_neighbours+1+field_grid_dim[2])%field_grid_dim[2]), current_partition_interval.second-1);
+ interp_limit_z[1] = std::min(int((partGridIdx_z+interp_neighbours+1+field_grid_dim[IDX_Z])%field_grid_dim[IDX_Z]), current_partition_interval.second-1);
nb_z_intervals = 2;
}
@@ -88,16 +89,16 @@ class particles_field_computer : public abstract_particles_distr<3,3,1> {
for(int idx_inter = 0 ; idx_inter < nb_z_intervals ; ++idx_inter){
for(int idx_z = interp_limit_mz[idx_inter] ; idx_z <= interp_limit_z[idx_inter] ; ++idx_z ){
- const int idx_z_pbc = (idx_z + field_grid_dim[2])%field_grid_dim[2];
+ const int idx_z_pbc = (idx_z + field_grid_dim[IDX_Z])%field_grid_dim[IDX_Z];
assert(current_partition_interval.first <= idx_z_pbc && idx_z_pbc < current_partition_interval.second);
assert(idx_z-interp_limit_mz[idx_inter] < interp_neighbours*2+2);
for(int idx_x = interp_limit_mx ; idx_x <= interp_limit_x ; ++idx_x ){
- const int idx_x_pbc = (idx_x + field_grid_dim[0])%field_grid_dim[0];
+ const int idx_x_pbc = (idx_x + field_grid_dim[IDX_X])%field_grid_dim[IDX_X];
assert(idx_x-interp_limit_mx < interp_neighbours*2+2);
for(int idx_y = interp_limit_my ; idx_y <= interp_limit_y ; ++idx_y ){
- const int idx_y_pbc = (idx_y + field_grid_dim[1])%field_grid_dim[1];
+ const int idx_y_pbc = (idx_y + field_grid_dim[IDX_Y])%field_grid_dim[IDX_Y];
assert(idx_y-interp_limit_my < interp_neighbours*2+2);
const double coef = (bz[idx_z-interp_limit_mz[idx_inter]]
@@ -106,9 +107,9 @@ class particles_field_computer : public abstract_particles_distr<3,3,1> {
const ptrdiff_t tindex = field.getIndexFromGlobalPosition(idx_x_pbc, idx_y_pbc, idx_z_pbc);
- particles_current_rhs[idxPart*3+0] += field.getValue(tindex,0)*coef;
- particles_current_rhs[idxPart*3+1] += field.getValue(tindex,1)*coef;
- particles_current_rhs[idxPart*3+2] += field.getValue(tindex,2)*coef;
+ particles_current_rhs[idxPart*3+IDX_X] += field.getValue(tindex,IDX_X)*coef;
+ particles_current_rhs[idxPart*3+IDX_Y] += field.getValue(tindex,IDX_Y)*coef;
+ particles_current_rhs[idxPart*3+IDX_Z] += field.getValue(tindex,IDX_Z)*coef;
}
}
}
@@ -123,9 +124,9 @@ class particles_field_computer : public abstract_particles_distr<3,3,1> {
TIMEZONE("particles_field_computer::reduce_particles");
// Simply sum values
for(int idxPart = 0 ; idxPart < nb_particles ; ++idxPart){
- particles_current_rhs[idxPart*3+0] += extra_particles_current_rhs[idxPart*3+0];
- particles_current_rhs[idxPart*3+1] += extra_particles_current_rhs[idxPart*3+1];
- particles_current_rhs[idxPart*3+2] += extra_particles_current_rhs[idxPart*3+2];
+ particles_current_rhs[idxPart*3+IDX_X] += extra_particles_current_rhs[idxPart*3+IDX_X];
+ particles_current_rhs[idxPart*3+IDX_Y] += extra_particles_current_rhs[idxPart*3+IDX_Y];
+ particles_current_rhs[idxPart*3+IDX_Z] += extra_particles_current_rhs[idxPart*3+IDX_Z];
}
}
@@ -136,9 +137,10 @@ class particles_field_computer : public abstract_particles_distr<3,3,1> {
void apply_pbc_xy(double* inout_particles, const int size) const final {
TIMEZONE("particles_field_computer::apply_pbc_xy");
+ const std::array<int, 2> dims_xy={IDX_X, IDX_Y};
for(int idxPart = 0 ; idxPart < size ; ++idxPart){
// Consider it will never move for more than one box repeatition
- for(int idxDim = 0 ; idxDim < 2 ; ++idxDim){
+ for(const int idxDim : dims_xy){
if(inout_particles[idxPart*3+idxDim] < 0) inout_particles[idxPart*3+idxDim] += spatial_box_width[idxDim];
else if(spatial_box_width[idxDim] <= inout_particles[idxPart*3+idxDim]) inout_particles[idxPart*3+idxDim] -= spatial_box_width[idxDim];
assert(0 <= inout_particles[idxPart*3+idxDim] && inout_particles[idxPart*3+idxDim] < spatial_box_width[idxDim]);
@@ -149,7 +151,7 @@ class particles_field_computer : public abstract_particles_distr<3,3,1> {
void apply_pbc_z_new_particles(double* values, const int size) const final {
TIMEZONE("particles_field_computer::apply_pbc_z_new_particles");
if(my_rank == 0){
- const int idxDim = 2;
+ const int idxDim = IDX_Z;
for(int idxPart = 0 ; idxPart < size ; ++idxPart){
assert(values[idxPart*3+idxDim] < my_spatial_up_limit || spatial_box_width[idxDim] <= values[idxPart*3+idxDim]);
assert(my_spatial_low_limit <= values[idxPart*3+idxDim]);
@@ -161,7 +163,7 @@ class particles_field_computer : public abstract_particles_distr<3,3,1> {
}
}
else if(my_rank == nb_processes - 1){
- const int idxDim = 2;
+ const int idxDim = IDX_Z;
for(int idxPart = 0 ; idxPart < size ; ++idxPart){
assert(my_spatial_low_limit <= values[idxPart*3+idxDim] || values[idxPart*3+idxDim] < 0);
assert(values[idxPart*3+idxDim] < spatial_box_width[idxDim]);
@@ -173,7 +175,7 @@ class particles_field_computer : public abstract_particles_distr<3,3,1> {
}
}
else{
- const int idxDim = 2;
+ const int idxDim = IDX_Z;
for(int idxPart = 0 ; idxPart < size ; ++idxPart){
assert(my_spatial_low_limit <= values[idxPart*3+idxDim] && values[idxPart*3+idxDim] < my_spatial_up_limit);
}
@@ -194,9 +196,9 @@ public:
interpolator(in_interpolator), field(in_field), positions_updater(),
spatial_box_width(in_spatial_box_width), box_step_width(in_box_step_width),
my_spatial_low_limit(in_my_spatial_low_limit), my_spatial_up_limit(in_my_spatial_up_limit){
- deriv[0] = 0;
- deriv[1] = 0;
- deriv[2] = 0;
+ deriv[IDX_X] = 0;
+ deriv[IDX_Y] = 0;
+ deriv[IDX_Z] = 0;
}
////////////////////////////////////////////////////////////////////////
diff --git a/bfps/cpp/particles/particles_input_hdf5.hpp b/bfps/cpp/particles/particles_input_hdf5.hpp
index c3b377a300e610d369938fdc8fd642f4ec7429c3..1ca9ae37d79a0018f6bc28486d680e7afeced69e 100644
--- a/bfps/cpp/particles/particles_input_hdf5.hpp
+++ b/bfps/cpp/particles/particles_input_hdf5.hpp
@@ -224,7 +224,7 @@ public:
&split_particles_positions[previousOffset*size_particle_positions],
load_splitter.getMySize()-previousOffset,
[&](const double val[]){
- return val[2] < limitPartition;
+ return val[IDX_Z] < limitPartition;
},
[&](const int idx1, const int idx2){
std::swap(split_particles_indexes[idx1], split_particles_indexes[idx2]);
diff --git a/bfps/cpp/particles/particles_system.hpp b/bfps/cpp/particles/particles_system.hpp
index 52240591d1de54330b8f9a2ed95b020b48dfa70c..72c460993889e7abfed2ce7f1b529fc661b2ddd2 100644
--- a/bfps/cpp/particles/particles_system.hpp
+++ b/bfps/cpp/particles/particles_system.hpp
@@ -3,6 +3,7 @@
#include <array>
+#include "abstract_particle_system.hpp"
#include "particles_output_hdf5.hpp"
#include "particles_output_mpiio.hpp"
#include "particles_field_computer.hpp"
@@ -12,7 +13,7 @@
#include "scope_timer.hpp"
template <class interpolator_class, int interp_neighbours>
-class particles_system {
+class particles_system : public abstract_particle_system {
MPI_Comm mpi_com;
const std::pair<int,int> current_partition_interval;
@@ -45,11 +46,12 @@ public:
const double in_my_spatial_low_limit, const double in_my_spatial_up_limit,
const double* in_field_data, const std::array<size_t,3>& in_local_field_dims,
const std::array<size_t,3>& in_local_field_offset,
+ const std::array<size_t,3>& in_field_memory_dims,
MPI_Comm in_mpi_com)
: mpi_com(in_mpi_com),
- current_partition_interval({in_local_field_offset[2], in_local_field_offset[2] + in_local_field_dims[2]}),
+ 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),
- field(in_field_data, in_local_field_dims, in_local_field_offset),
+ field(in_field_data, in_local_field_dims, in_local_field_offset, in_field_memory_dims),
interpolator(),
computer(in_mpi_com, field_grid_dim, current_partition_interval,
interpolator, field, in_spatial_box_width, in_spatial_partition_width,
@@ -74,8 +76,8 @@ public:
my_nb_particles = particles_input.getLocalNbParticles();
for(int idx_part = 0 ; idx_part < my_nb_particles ; ++idx_part){ // TODO remove me
- assert(my_particles_positions[idx_part*3+2] >= my_spatial_low_limit);
- assert(my_particles_positions[idx_part*3+2] < my_spatial_up_limit);
+ assert(my_particles_positions[idx_part*3+IDX_Z] >= my_spatial_low_limit);
+ assert(my_particles_positions[idx_part*3+IDX_Z] < my_spatial_up_limit);
}
particles_utils::partition_extra_z<3>(&my_particles_positions[0], my_nb_particles, partition_interval_size,
@@ -100,17 +102,17 @@ public:
current_offset_particles_for_partition[idxPartition+1] - current_offset_particles_for_partition[idxPartition]);
const double limitPartition = (idxPartition+1)*spatial_partition_width + my_spatial_low_limit;
for(int idx = 0 ; idx < current_offset_particles_for_partition[idxPartition+1] ; ++idx){
- assert(my_particles_positions[idx*3+2] < limitPartition);
+ assert(my_particles_positions[idx*3+IDX_Z] < limitPartition);
}
for(int idx = current_offset_particles_for_partition[idxPartition+1] ; idx < my_nb_particles ; ++idx){
- assert(my_particles_positions[idx*3+2] >= limitPartition);
+ assert(my_particles_positions[idx*3+IDX_Z] >= limitPartition);
}
}
}
}
- void compute(){
+ void compute() final {
TIMEZONE("particles_system::compute");
computer.compute_distr(current_my_nb_particles_per_partition.get(),
my_particles_positions.get(),
@@ -118,14 +120,14 @@ public:
interp_neighbours);
}
- void move(const double dt){
+ void move(const double dt) final {
TIMEZONE("particles_system::move");
computer.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(){
+ void redistribute() final {
TIMEZONE("particles_system::redistribute");
computer.redistribute(current_my_nb_particles_per_partition.get(),
&my_nb_particles,
@@ -137,11 +139,11 @@ public:
spatial_partition_width);
}
- void inc_step_idx(){
+ void inc_step_idx() final {
step_idx += 1;
}
- void shift_rhs_vectors(){
+ void shift_rhs_vectors() final {
if(my_particles_rhs.size()){
std::unique_ptr<double[]> next_current(std::move(my_particles_rhs.back()));
for(int idx_rhs = my_particles_rhs.size()-1 ; idx_rhs > 0 ; --idx_rhs){
@@ -152,7 +154,7 @@ public:
}
}
- void completeLoop(const double dt){
+ void completeLoop(const double dt) final {
TIMEZONE("particles_system::completeLoop");
compute();
move(dt);
@@ -161,32 +163,32 @@ public:
shift_rhs_vectors();
}
- const double* getParticlesPositions() const{
+ const double* getParticlesPositions() const final {
return my_particles_positions.get();
}
- const double* getParticlesCurrentRhs() const{
+ const double* getParticlesCurrentRhs() const final {
return my_particles_rhs.front().get();
}
- const int* getParticlesIndexes() const{
+ const int* getParticlesIndexes() const final {
return my_particles_positions_indexes.get();
}
- int getLocalNbParticles() const{
+ int getLocalNbParticles() const final {
return my_nb_particles;
}
void checkNan() const { // TODO remove
for(int idx_part = 0 ; idx_part < my_nb_particles ; ++idx_part){ // TODO remove me
- assert(std::isnan(my_particles_positions[idx_part*3+0]) == false);
- assert(std::isnan(my_particles_positions[idx_part*3+1]) == false);
- assert(std::isnan(my_particles_positions[idx_part*3+2]) == false);
+ assert(std::isnan(my_particles_positions[idx_part*3+IDX_X]) == false);
+ assert(std::isnan(my_particles_positions[idx_part*3+IDX_Y]) == false);
+ assert(std::isnan(my_particles_positions[idx_part*3+IDX_Z]) == false);
for(int idx_rhs = 0 ; idx_rhs < my_particles_rhs.size() ; ++idx_rhs){
- assert(std::isnan(my_particles_rhs[idx_rhs][idx_part*3+0]) == false);
- assert(std::isnan(my_particles_rhs[idx_rhs][idx_part*3+1]) == false);
- assert(std::isnan(my_particles_rhs[idx_rhs][idx_part*3+2]) == false);
+ assert(std::isnan(my_particles_rhs[idx_rhs][idx_part*3+IDX_X]) == false);
+ assert(std::isnan(my_particles_rhs[idx_rhs][idx_part*3+IDX_Y]) == false);
+ assert(std::isnan(my_particles_rhs[idx_rhs][idx_part*3+IDX_Z]) == false);
}
}
}
diff --git a/bfps/cpp/particles/particles_utils.hpp b/bfps/cpp/particles/particles_utils.hpp
index eb3d22ab15e718026f53cf78bb1ea659e60bc6d1..3302a9f910b2e82a7236245c63899834cc9e8414 100644
--- a/bfps/cpp/particles/particles_utils.hpp
+++ b/bfps/cpp/particles/particles_utils.hpp
@@ -4,6 +4,11 @@
#include <cassert>
#include <stack>
+enum IDXS_3D {
+ IDX_X = 2,
+ IDX_Y = 1,
+ IDX_Z = 0
+};
namespace particles_utils {
@@ -78,7 +83,7 @@ inline void partition_extra_z(double* array, const int size, const int nb_partit
const double limit = partitions_limits(0);
const int size_current = partition_extra<nb_values>(array, size,
[&](const double inval[]){
- return inval[nb_values-1] < limit;
+ return inval[IDX_Z] < limit;
}, pdcswap);
partitions_size[0] = size_current;
partitions_size[1] = size-size_current;
@@ -108,7 +113,7 @@ inline void partition_extra_z(double* array, const int size, const int nb_partit
const int size_current = partition_extra<nb_values>(&array[partitions_offset[current_part.first]*nb_values],
size_unpart,
[&](const double inval[]){
- return inval[nb_values-1] < limit;
+ return inval[IDX_Z] < limit;
}, pdcswap, partitions_offset[current_part.first]);
partitions_offset[idx_middle+1] = size_current + partitions_offset[current_part.first];
@@ -185,7 +190,9 @@ public:
}
NumType getOwner(const NumType in_item_idx) const {
- return NumType(double(in_item_idx)/step_split);
+ const NumType owner = NumType(double(in_item_idx)/step_split);
+ assert(owner < nb_intervals);
+ return owner;
}
};