diff --git a/bfps/cpp/particles/main_tocompile.cpp b/bfps/cpp/particles/main_tocompile.cpp
index 25a72e4905e597af43ac7a5030f082745b961360..fb62c8f6868da503bb928d898d37c1f8c7a26f66 100644
--- a/bfps/cpp/particles/main_tocompile.cpp
+++ b/bfps/cpp/particles/main_tocompile.cpp
@@ -98,13 +98,17 @@ int main(int argc, char** argv){
const std::array<double,3> spatial_box_width{10., 10., 10.};
- 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;
+ std::array<double,3> spatial_partition_width;
+ spatial_partition_width[IDX_X] = spatial_box_width[IDX_X]/double(field_grid_dim[IDX_X]);
+ spatial_partition_width[IDX_Y] = spatial_box_width[IDX_Y]/double(field_grid_dim[IDX_Y]);
+ spatial_partition_width[IDX_Z] = spatial_box_width[IDX_Z]/double(field_grid_dim[IDX_Z]);
+
+ const double my_spatial_low_limit = myPartitionInterval[0]*spatial_partition_width[IDX_Z];
+ const double my_spatial_up_limit = myPartitionInterval[1]*spatial_partition_width[IDX_Z];
if(my_rank == 0){
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 << "spatial_partition_width = " << spatial_partition_width[IDX_Z] << 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;
}
@@ -140,7 +144,7 @@ int main(int argc, char** argv){
// my_spatial_up_limit, my_rank, nb_processes);
std::vector<double> spatial_interval_per_proc(nb_processes+1);
for(int idx_proc = 0 ; idx_proc < nb_processes ; ++idx_proc){
- spatial_interval_per_proc[idx_proc] = partitionIntervalSize*spatial_partition_width*idx_proc;
+ spatial_interval_per_proc[idx_proc] = partitionIntervalSize*spatial_partition_width[IDX_Z]*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[IDX_Z];
diff --git a/bfps/cpp/particles/particles_field_computer.hpp b/bfps/cpp/particles/particles_field_computer.hpp
index b0f9077a5ab80f451b670de5556e721d43fde3a0..1e39d63a62049c40d7f58fa6b19d46b845e78d52 100644
--- a/bfps/cpp/particles/particles_field_computer.hpp
+++ b/bfps/cpp/particles/particles_field_computer.hpp
@@ -21,9 +21,9 @@ class particles_field_computer : public abstract_particles_distr<real_number, 3,
const positions_updater_class positions_updater;
const std::array<real_number,3> spatial_box_width;
- const real_number box_step_width;
- const real_number my_spatial_low_limit;
- const real_number my_spatial_up_limit;
+ const std::array<real_number,3> box_step_width;
+ const real_number my_spatial_low_limit_z;
+ const real_number my_spatial_up_limit_z;
int deriv[3];
@@ -37,19 +37,30 @@ class particles_field_computer : public abstract_particles_distr<real_number, 3,
std::fill(particles_current_rhs, particles_current_rhs+nb_particles*3, 0);
}
+ real_number get_norm_pos_in_cell(const real_number in_pos, const int idx_pos) const {
+ const real_number cell_idx = floor(in_pos/box_step_width[idx_pos]);
+ const real_number pos_in_cell = (in_pos - cell_idx*box_step_width[idx_pos]) / box_step_width[idx_pos];
+ assert(0 <= pos_in_cell && pos_in_cell < box_step_width[idx_pos]);
+ return pos_in_cell;
+ }
+
virtual void apply_computation(const real_number particles_positions[],
real_number particles_current_rhs[],
const int nb_particles) const final{
TIMEZONE("particles_field_computer::apply_computation");
for(int idxPart = 0 ; idxPart < nb_particles ; ++idxPart){
+ const real_number reltv_x = get_norm_pos_in_cell(particles_positions[idxPart*3+IDX_X], IDX_X);
+ const real_number reltv_y = get_norm_pos_in_cell(particles_positions[idxPart*3+IDX_Y], IDX_Y);
+ const real_number reltv_z = get_norm_pos_in_cell(particles_positions[idxPart*3+IDX_Z], IDX_Z);
+
typename interpolator_class::real_number bx[interp_neighbours*2+2], by[interp_neighbours*2+2], bz[interp_neighbours*2+2];
- 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);
+ interpolator.compute_beta(deriv[IDX_X], reltv_x, bx);
+ interpolator.compute_beta(deriv[IDX_Y], reltv_y, by);
+ interpolator.compute_beta(deriv[IDX_Z], reltv_z, bz);
- 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);
+ const int partGridIdx_x = int(particles_positions[idxPart*3+IDX_X]/box_step_width[IDX_X]);
+ const int partGridIdx_y = int(particles_positions[idxPart*3+IDX_Y]/box_step_width[IDX_Y]);
+ const int partGridIdx_z = int(particles_positions[idxPart*3+IDX_Z]/box_step_width[IDX_Z]);
assert(0 <= partGridIdx_x && partGridIdx_x < field_grid_dim[IDX_X]);
assert(0 <= partGridIdx_y && partGridIdx_y < field_grid_dim[IDX_Y]);
@@ -155,31 +166,31 @@ class particles_field_computer : public abstract_particles_distr<real_number, 3,
if(Parent::my_rank == 0){
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]);
+ assert(values[idxPart*3+idxDim] < my_spatial_up_limit_z || spatial_box_width[idxDim] <= values[idxPart*3+idxDim]);
+ assert(my_spatial_low_limit_z <= values[idxPart*3+idxDim]);
if(spatial_box_width[idxDim] <= values[idxPart*3+idxDim]) values[idxPart*3+idxDim] -= spatial_box_width[idxDim];
assert(0 <= values[idxPart*3+idxDim] && values[idxPart*3+idxDim] < spatial_box_width[idxDim]);
- assert(my_spatial_low_limit <= values[idxPart*3+idxDim] && values[idxPart*3+idxDim] < my_spatial_up_limit);
+ assert(my_spatial_low_limit_z <= values[idxPart*3+idxDim] && values[idxPart*3+idxDim] < my_spatial_up_limit_z);
}
}
else if(Parent::my_rank == Parent::nb_processes - 1){
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(my_spatial_low_limit_z <= values[idxPart*3+idxDim] || values[idxPart*3+idxDim] < 0);
assert(values[idxPart*3+idxDim] < spatial_box_width[idxDim]);
if(values[idxPart*3+idxDim] < 0) values[idxPart*3+idxDim] += spatial_box_width[idxDim];
assert(0 <= values[idxPart*3+idxDim] && values[idxPart*3+idxDim] < spatial_box_width[idxDim]);
- assert(my_spatial_low_limit <= values[idxPart*3+idxDim] && values[idxPart*3+idxDim] < my_spatial_up_limit);
+ assert(my_spatial_low_limit_z <= values[idxPart*3+idxDim] && values[idxPart*3+idxDim] < my_spatial_up_limit_z);
}
}
else{
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);
+ assert(my_spatial_low_limit_z <= values[idxPart*3+idxDim] && values[idxPart*3+idxDim] < my_spatial_up_limit_z);
}
}
}
@@ -191,13 +202,13 @@ public:
const interpolator_class& in_interpolator,
const field_class& in_field,
const std::array<real_number,3>& in_spatial_box_width,
- const real_number in_box_step_width, const real_number in_my_spatial_low_limit,
- const real_number in_my_spatial_up_limit)
+ const std::array<real_number,3>& in_box_step_width, const real_number in_my_spatial_low_limit_z,
+ const real_number in_my_spatial_up_limit_z)
: abstract_particles_distr<real_number, 3,3,1>(in_current_com, in_current_partitions),
field_grid_dim(in_field_grid_dim), current_partition_interval(in_current_partitions),
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){
+ my_spatial_low_limit_z(in_my_spatial_low_limit_z), my_spatial_up_limit_z(in_my_spatial_up_limit_z){
deriv[IDX_X] = 0;
deriv[IDX_Y] = 0;
deriv[IDX_Z] = 0;
diff --git a/bfps/cpp/particles/particles_system.hpp b/bfps/cpp/particles/particles_system.hpp
index dfc5683c315c3d9a404297d143192863275c5411..7a9c3f6ffe0c2a28490a197527c72d279bf08374 100644
--- a/bfps/cpp/particles/particles_system.hpp
+++ b/bfps/cpp/particles/particles_system.hpp
@@ -29,7 +29,7 @@ class particles_system : public abstract_particles_system<real_number> {
std::unique_ptr<int[]> current_offset_particles_for_partition;
const std::array<real_number,3> spatial_box_width;
- const real_number spatial_partition_width;
+ const std::array<real_number,3> spatial_partition_width;
const real_number my_spatial_low_limit;
const real_number my_spatial_up_limit;
@@ -42,7 +42,7 @@ class particles_system : public abstract_particles_system<real_number> {
public:
particles_system(const std::array<size_t,3>& field_grid_dim, const std::array<real_number,3>& in_spatial_box_width,
- const real_number in_spatial_partition_width,
+ 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 real_number* in_field_data, const std::array<size_t,3>& in_local_field_dims,
const std::array<size_t,3>& in_local_field_offset,
@@ -83,7 +83,7 @@ public:
particles_utils::partition_extra_z<3>(&my_particles_positions[0], my_nb_particles, partition_interval_size,
current_my_nb_particles_per_partition.get(), current_offset_particles_for_partition.get(),
[&](const int idxPartition){
- const real_number limitPartition = (idxPartition+1)*spatial_partition_width + my_spatial_low_limit;
+ const real_number limitPartition = (idxPartition+1)*spatial_partition_width[IDX_Z] + my_spatial_low_limit;
return limitPartition;
},
[&](const int idx1, const int idx2){
@@ -100,7 +100,7 @@ public:
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]);
- const real_number limitPartition = (idxPartition+1)*spatial_partition_width + my_spatial_low_limit;
+ const real_number limitPartition = (idxPartition+1)*spatial_partition_width[IDX_Z] + my_spatial_low_limit;
for(int idx = 0 ; idx < current_offset_particles_for_partition[idxPartition+1] ; ++idx){
assert(my_particles_positions[idx*3+IDX_Z] < limitPartition);
}
@@ -136,7 +136,7 @@ public:
&my_particles_positions_indexes,
my_spatial_low_limit,
my_spatial_up_limit,
- spatial_partition_width);
+ spatial_partition_width[IDX_Z]);
}
void inc_step_idx() final {
diff --git a/bfps/cpp/particles/particles_system_builder.hpp b/bfps/cpp/particles/particles_system_builder.hpp
index dd463c0a8820d657ed7e8621d58a6ccbfac1dadf..fbde1fa57317a07445f193509e51cc15354f7aab 100644
--- a/bfps/cpp/particles/particles_system_builder.hpp
+++ b/bfps/cpp/particles/particles_system_builder.hpp
@@ -166,16 +166,19 @@ struct particles_system_build_container {
spatial_box_width[IDX_Z] = fs_kk->dkz;
// The distance between two field nodes in z
- const rnumber spatial_partition_width_z = spatial_box_width[IDX_Z]/rnumber(field_grid_dim[IDX_Z]);
+ std::array<rnumber,3> spatial_partition_width;
+ spatial_partition_width[IDX_X] = spatial_box_width[IDX_X]/rnumber(field_grid_dim[IDX_X]);
+ spatial_partition_width[IDX_Y] = spatial_box_width[IDX_Y]/rnumber(field_grid_dim[IDX_Y]);
+ spatial_partition_width[IDX_Z] = spatial_box_width[IDX_Z]/rnumber(field_grid_dim[IDX_Z]);
// The spatial interval of the current process
- const rnumber my_spatial_low_limit_z = rnumber(local_field_offset[IDX_Z])*spatial_partition_width_z;
- const rnumber my_spatial_up_limit_z = rnumber(local_field_offset[IDX_Z]+local_field_dims[IDX_Z])*spatial_partition_width_z;
+ const rnumber my_spatial_low_limit_z = rnumber(local_field_offset[IDX_Z])*spatial_partition_width[IDX_Z];
+ const rnumber my_spatial_up_limit_z = rnumber(local_field_offset[IDX_Z]+local_field_dims[IDX_Z])*spatial_partition_width[IDX_Z];
// Create the particles system
particles_system<rnumber, particles_interp_spline<double, interpolation_size,spline_mode>, interpolation_size>* part_sys
= new particles_system<rnumber, particles_interp_spline<double, interpolation_size,spline_mode>, interpolation_size>(field_grid_dim,
spatial_box_width,
- spatial_partition_width_z,
+ spatial_partition_width,
my_spatial_low_limit_z,
my_spatial_up_limit_z,
fs_cvorticity->get_rdata(),