Skip to content
Snippets Groups Projects
Commit 93462d10 authored by Berenger Bramas's avatar Berenger Bramas
Browse files

Call the spline using the position of the particle inside the cell (normalize... 

Call the spline using the position of the particle inside the cell (normalize to 0-1) -- pass an array of box width to the constructor
parent c2b2cac8
Branches
Tags
2 merge requests!21Bugfix/nansampling,!3Bugfix/event manager show html
Show whitespace changes
Inline Side-by-side
bfps/cpp/particles/main_tocompile.cpp
+ 9
− 5
View file @ 93462d10
......@@ -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];
......
......
bfps/cpp/particles/particles_field_computer.hpp
+ 29
− 18
View file @ 93462d10
......@@ -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;
......
......
bfps/cpp/particles/particles_system.hpp
+ 5
− 5
View file @ 93462d10
......@@ -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 {
......
......
bfps/cpp/particles/particles_system_builder.hpp
+ 7
− 4
View file @ 93462d10
......@@ -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(),
......
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment