From e2e52fa062475eeb23073e79fc8beec43bda5971 Mon Sep 17 00:00:00 2001
From: Berenger Bramas <berenger.bramas@mpcdf.mpg.de>
Date: Wed, 20 Sep 2017 17:25:01 +0200
Subject: [PATCH] Debug
---
bfps/cpp/particles/p2p_distr_mpi.hpp | 240 +++++++++++++++------------
bfps/cpp/particles/p2p_tree.hpp | 1 -
2 files changed, 132 insertions(+), 109 deletions(-)
diff --git a/bfps/cpp/particles/p2p_distr_mpi.hpp b/bfps/cpp/particles/p2p_distr_mpi.hpp
index a32c3375..45f8bba3 100644
--- a/bfps/cpp/particles/p2p_distr_mpi.hpp
+++ b/bfps/cpp/particles/p2p_distr_mpi.hpp
@@ -25,6 +25,50 @@
- update particles property
*/
+
+
+template <class partsize_t, class real_number, int size_particle_positions, int size_particle_rhs>
+struct ParticleView{
+ partsize_t p_index;
+ real_number* ptr_particles_positions;
+ real_number* ptr_particles_current_rhs;
+ partsize_t* ptr_global_idx;
+ long int* ptr_cell_idx;
+
+ ParticleView()
+ : p_index(-1), ptr_particles_positions(nullptr),
+ ptr_particles_current_rhs(nullptr), ptr_global_idx(nullptr),
+ ptr_cell_idx(nullptr){}
+};
+
+template <class partsize_t, class real_number, int size_particle_positions, int size_particle_rhs>
+void swap(ParticleView<partsize_t, real_number, size_particle_positions,size_particle_rhs>& p1,
+ ParticleView<partsize_t, real_number, size_particle_positions,size_particle_rhs>& p2){
+ if(p1.p_index != -1 && p2.p_index != -1){
+ for(int idx_pos = 0 ; idx_pos < size_particle_positions ; ++idx_pos){
+ std::swap(p1.ptr_particles_positions[p1.p_index*size_particle_positions+idx_pos],
+ p2.ptr_particles_positions[p2.p_index*size_particle_positions+idx_pos]);
+ }
+ for(int idx_rhs = 0 ; idx_rhs < size_particle_rhs ; ++idx_rhs){
+ std::swap(p1.ptr_particles_current_rhs[p1.p_index*size_particle_rhs+idx_rhs],
+ p2.ptr_particles_current_rhs[p2.p_index*size_particle_rhs+idx_rhs]);
+ }
+ std::swap(p1.ptr_cell_idx[p1.p_index],p2.ptr_cell_idx[p2.p_index]);
+ std::swap(p1.ptr_global_idx[p1.p_index],p2.ptr_global_idx[p2.p_index]);
+ std::swap(p1.p_index,p2.p_index);
+ }
+ else if(p1.p_index != -1){
+ p2 = p1;
+ p1 = ParticleView<partsize_t, real_number, size_particle_positions,size_particle_rhs>();
+ }
+ else if(p2.p_index != -1){
+ p1 = p2;
+ p2 = ParticleView<partsize_t, real_number, size_particle_positions,size_particle_rhs>();
+ }
+}
+
+
+
template <class partsize_t, class real_number>
class p2p_distr_mpi {
protected:
@@ -183,32 +227,23 @@ public:
real_number compute_distance_r2(const real_number x1, const real_number y1, const real_number z1,
const real_number x2, const real_number y2, const real_number z2) const {
- return ((x1-x2)*(x1-x2)) + ((y1-y2)*(y1-y2)) + ((z1-z2)*(z1-z2));
- }
-
+ real_number diff_x = std::abs(x1-x2);
+ while(diff_x > spatial_box_width[IDX_X]/2){
+ diff_x = std::abs(diff_x - spatial_box_width[IDX_X]);
+ }
- template <int size_particle_positions, int size_particle_rhs>
- struct ParticleView{
- partsize_t p_index;
- real_number* ptr_particles_positions;
- real_number* ptr_particles_current_rhs;
- partsize_t* ptr_global_idx;
- long int* ptr_cell_idx;
+ real_number diff_y = std::abs(y1-y2);
+ while(diff_y > spatial_box_width[IDX_Y]/2){
+ diff_y = std::abs(diff_y - spatial_box_width[IDX_Y]);
+ }
- void swap(ParticleView& p1, ParticleView& p2){
- for(int idx_pos = 0 ; idx_pos < size_particle_positions ; ++idx_pos){
- std::swap(p1.ptr_particles_positions[p1.p_index*size_particle_positions+idx_pos],
- p2.ptr_particles_positions[p2.p_index*size_particle_positions+idx_pos]);
- }
- for(int idx_rhs = 0 ; idx_rhs < size_particle_rhs ; ++idx_rhs){
- std::swap(p1.ptr_particles_current_rhs[p1.p_index*size_particle_rhs+idx_rhs],
- p2.ptr_particles_current_rhs[p2.p_index*size_particle_rhs+idx_rhs]);
- }
- std::swap(p1.ptr_cell_idx[p1.p_index],p2.ptr_cell_idx[p2.p_index]);
- std::swap(p1.ptr_global_idx[p1.p_index],p2.ptr_global_idx[p2.p_index]);
- std::swap(p1.p_index,p2.p_index);
+ real_number diff_z = std::abs(z1-z2);
+ while(diff_z > spatial_box_width[IDX_Z]/2){
+ diff_z = std::abs(diff_z - spatial_box_width[IDX_Z]);
}
- };
+
+ return (diff_x*diff_x) + (diff_y*diff_y) + (diff_z*diff_z);
+ }
template <class computer_class, int size_particle_positions, int size_particle_rhs>
void compute_distr(computer_class& in_computer,
@@ -246,16 +281,18 @@ public:
particles_positions[(idxPart)*size_particle_positions + IDX_Z]);
assert(my_down_z_cell_level <= get_cell_coord_z_from_index(particles_coord[idxPart]));
assert(get_cell_coord_z_from_index(particles_coord[idxPart]) <= my_top_z_cell_level);
- if(inout_index_particles[idxPart] == 3 || inout_index_particles[idxPart] == 52){// TODO
- printf("Coord index %ld (tree index %ld)\n", inout_index_particles[idxPart],particles_coord[idxPart]);
+ if(inout_index_particles[idxPart] == 58576 || inout_index_particles[idxPart] == 0){// TODO
+ printf("Coord index %ld - %ld (tree index %ld)\n", idxPart, inout_index_particles[idxPart],particles_coord[idxPart]);
printf(">> Box index %ld - %ld - %ld\n", get_cell_coord_x_from_index(particles_coord[idxPart]),
get_cell_coord_y_from_index(particles_coord[idxPart]),
get_cell_coord_z_from_index(particles_coord[idxPart]));
+ printf("idxPartition %d\n", idxPartition);
}
}
}
- std::vector<ParticleView<size_particle_positions,size_particle_rhs>> part_to_sort;
+ using ParticleView_t = ParticleView<partsize_t, real_number, size_particle_positions,size_particle_rhs>;
+ std::vector<ParticleView_t> part_to_sort;
// Sort each partition in cells
for(int idxPartition = 0 ; idxPartition < current_partition_size ; ++idxPartition){
@@ -273,8 +310,10 @@ public:
assert(part_to_sort.size() == (current_offset_particles_for_partition[idxPartition+1]-current_offset_particles_for_partition[idxPartition]));
std::sort(part_to_sort.begin(), part_to_sort.end(),
- [](const ParticleView<size_particle_positions,size_particle_rhs>& p1,
- const ParticleView<size_particle_positions,size_particle_rhs>& p2){
+ [](const ParticleView_t& p1,
+ const ParticleView_t& p2){
+ assert(p1.p_index != -1 && p1.ptr_cell_idx);
+ assert(p2.p_index != -1 && p2.ptr_cell_idx);
return p1.ptr_cell_idx[p1.p_index] < p2.ptr_cell_idx[p2.p_index];
});
}
@@ -297,6 +336,13 @@ public:
current_cell_idx = particles_coord[idx_part];
current_nb_particles_in_cell = 1;
current_cell_offset = idx_part;
+ if(inout_index_particles[idx_part] == 58576 || inout_index_particles[idx_part] == 0){// TODO
+ printf("Coord index %ld - %ld (tree index %ld)\n", idx_part, inout_index_particles[idx_part],particles_coord[idx_part]);
+ printf(">> Box index %ld - %ld - %ld\n", get_cell_coord_x_from_index(particles_coord[idx_part]),
+ get_cell_coord_y_from_index(particles_coord[idx_part]),
+ get_cell_coord_z_from_index(particles_coord[idx_part]));
+ printf("current_cell_offset %ld current_nb_particles_in_cell %ld\n", current_cell_offset, current_nb_particles_in_cell);
+ }
}
}
if(current_nb_particles_in_cell){
@@ -310,13 +356,17 @@ public:
fflush(stdout); // TODO
// Offset per cell layers
+ long int previous_index = 0;
std::unique_ptr<partsize_t[]> particles_offset_layers(new partsize_t[my_nb_cell_levels+1]());
for(int idxPartition = 0 ; idxPartition < current_partition_size ; ++idxPartition){
for(partsize_t idx_part = current_offset_particles_for_partition[idxPartition] ;
idx_part != current_offset_particles_for_partition[idxPartition+1]; ++idx_part){
- assert(my_down_z_cell_level <= get_cell_coord_z_from_index(particles_coord[idx_part]));
- assert(get_cell_coord_z_from_index(particles_coord[idx_part]) <= my_top_z_cell_level);
- particles_offset_layers[get_cell_coord_z_from_index(particles_coord[idx_part])+1-my_down_z_cell_level] += 1;
+ const long int part_box_z_index = get_cell_coord_z_from_index(particles_coord[idx_part]);
+ assert(my_down_z_cell_level <= part_box_z_index);
+ assert(part_box_z_index <= my_top_z_cell_level);
+ particles_offset_layers[part_box_z_index+1-my_down_z_cell_level] += 1;
+ assert(previous_index <= part_box_z_index);
+ previous_index = part_box_z_index;
}
}
for(size_t idx_layer = 0 ; idx_layer < my_nb_cell_levels ; ++idx_layer){
@@ -363,7 +413,6 @@ public:
std::cout << "descriptor.destProc " << descriptor.destProc << std::endl;
std::cout << "descriptor.nbLevelsToExchange " << descriptor.nbLevelsToExchange << std::endl;
std::cout << "descriptor.nbParticlesToExchange " << descriptor.nbParticlesToExchange << std::endl;
- std::cout << "descriptor.nbParticlesToExchange " << descriptor.nbParticlesToExchange << std::endl;
std::cout << "descriptor.isRecv " << descriptor.isRecv << std::endl;
std::cout << "neigDescriptors.size() " << neigDescriptors.size() << std::endl;
std::cout.flush();// TODO
@@ -440,18 +489,19 @@ public:
std::cout << "descriptor.nbParticlesToExchange " << descriptor.nbParticlesToExchange << std::endl;
std::cout << "descriptor.destProc " << descriptor.destProc << std::endl;
std::cout << "idxDescr " << idxDescr << std::endl;
+ std::cout << "send from part " << particles_offset_layers[my_nb_cell_levels-descriptor.nbLevelsToExchange] << std::endl;
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
assert(descriptor.nbParticlesToExchange*size_particle_positions < std::numeric_limits<int>::max());
- AssertMpi(MPI_Isend(const_cast<real_number*>(&particles_positions[particles_offset_layers[my_nb_cell_levels-descriptor.nbLevelsToExchange]]),
+ AssertMpi(MPI_Isend(const_cast<real_number*>(&particles_positions[particles_offset_layers[my_nb_cell_levels-descriptor.nbLevelsToExchange]*size_particle_positions]),
int(descriptor.nbParticlesToExchange*size_particle_positions), particles_utils::GetMpiType(real_number()),
descriptor.destProc, TAG_POSITION_PARTICLES,
current_com, &mpiRequests.back()));
assert(descriptor.toRecvAndMerge == nullptr);
descriptor.toRecvAndMerge.reset(new real_number[descriptor.nbParticlesToExchange*size_particle_rhs]);
- whatNext.emplace_back(std::pair<Action,int>{MERGE_PARTICLES, int(neigDescriptors.size())});
+ whatNext.emplace_back(std::pair<Action,int>{MERGE_PARTICLES, idxDescr});
mpiRequests.emplace_back();
assert(descriptor.nbParticlesToExchange*size_particle_rhs < std::numeric_limits<int>::max());
AssertMpi(MPI_Irecv(descriptor.toRecvAndMerge.get(), int(descriptor.nbParticlesToExchange*size_particle_rhs),
@@ -470,8 +520,6 @@ public:
}
}
- const bool more_than_one_thread = (omp_get_max_threads() > 1);
-
TIMEZONE_OMP_INIT_PREPARALLEL(omp_get_max_threads())
#pragma omp parallel default(shared)
{
@@ -496,8 +544,7 @@ public:
//////////////////////////////////////////////////////////////////////
if(releasedAction.first == RECV_PARTICLES){
NeighborDescriptor& descriptor = neigDescriptors[releasedAction.second];
-
- assert(descriptor.isRecv == true);
+ assert(descriptor.isRecv);
const int destProc = descriptor.destProc;
const partsize_t NbParticlesToReceive = descriptor.nbParticlesToExchange;
assert(NbParticlesToReceive != -1);
@@ -525,11 +572,12 @@ public:
//////////////////////////////////////////////////////////////////////
if(releasedAction.first == COMPUTE_PARTICLES){
NeighborDescriptor& descriptor = neigDescriptors[releasedAction.second];
+ assert(descriptor.isRecv);
const partsize_t NbParticlesToReceive = descriptor.nbParticlesToExchange;
assert(descriptor.toCompute != nullptr);
descriptor.results.reset(new real_number[NbParticlesToReceive*size_particle_rhs]);
- // TODO in_computer.template init_result_array<size_particle_rhs>(descriptor.results.get(), NbParticlesToReceive);
+ in_computer.template init_result_array<size_particle_rhs>(descriptor.results.get(), NbParticlesToReceive);
// Compute
partsize_t idxPart = 0;
@@ -537,7 +585,7 @@ public:
const long int current_cell_idx = get_cell_idx(descriptor.toCompute[idxPart*size_particle_positions + IDX_X],
descriptor.toCompute[idxPart*size_particle_positions + IDX_Y],
descriptor.toCompute[idxPart*size_particle_positions + IDX_Z]);
- partsize_t nb_parts_in_cell = 0;
+ partsize_t nb_parts_in_cell = 1;
while(idxPart+nb_parts_in_cell != NbParticlesToReceive
&& current_cell_idx == get_cell_idx(descriptor.toCompute[(idxPart+nb_parts_in_cell)*size_particle_positions + IDX_X],
descriptor.toCompute[(idxPart+nb_parts_in_cell)*size_particle_positions + IDX_Y],
@@ -549,6 +597,23 @@ public:
long int neighbors_indexes[27];
const int nbNeighbors = my_tree.getNeighbors(current_cell_idx, neighbors, neighbors_indexes, true);
+ for(int idx_test = 0 ; idx_test < nb_parts_in_cell ; ++idx_test){ // TODO
+ if(int(descriptor.toCompute[(idxPart+idx_test)*size_particle_positions + IDX_X]*1000) == int(1.685800e-01*1000)
+ && int(descriptor.toCompute[(idxPart+idx_test)*size_particle_positions + IDX_Y]*1000) == int(7.524981e-01*1000)
+ && int(descriptor.toCompute[(idxPart+idx_test)*size_particle_positions + IDX_Z]*1000) == int(9.999596e-01*1000)){
+ printf("Found a pos %ld\n", idxPart+idx_test);
+ printf("pos %e %e %e\n",
+ descriptor.toCompute[(idxPart+idx_test)*size_particle_positions + IDX_X],
+ descriptor.toCompute[(idxPart+idx_test)*size_particle_positions + IDX_Y],
+ descriptor.toCompute[(idxPart+idx_test)*size_particle_positions + IDX_Z]);
+ }
+ }
+ printf("Remote part from %ld for %ld at idx %ld\n", idxPart, nb_parts_in_cell, current_cell_idx); // TODO
+ printf("pos of first %e %e %e\n", descriptor.toCompute[idxPart*size_particle_positions + IDX_X],
+ descriptor.toCompute[idxPart*size_particle_positions + IDX_Y],
+ descriptor.toCompute[idxPart*size_particle_positions + IDX_Z]); // TODO
+ printf("nbNeighbors %d\n", nbNeighbors); // TODO
+
// with other interval
for(size_t idx_neighbor = 0 ; idx_neighbor < nbNeighbors ; ++idx_neighbor){
for(size_t idx_2 = 0 ; idx_2 < (*neighbors[idx_neighbor]).size() ; ++idx_2){
@@ -561,12 +626,28 @@ public:
particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Y],
particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Z]);
if(dist_r2 < cutoff_radius*cutoff_radius){
- // TODO in_computer.template compute_interaction<size_particle_positions,size_particle_rhs>(
-// &descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions],
-// &descriptor.results[(idxPart+idx_p1)*size_particle_rhs],
-// &particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions],
-// &particles_current_rhs[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_rhs],
-// dist_r2);
+ in_computer.template compute_interaction<size_particle_positions,size_particle_rhs>(
+ &descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions],
+ &descriptor.results[(idxPart+idx_p1)*size_particle_rhs],
+ &particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions],
+ &particles_current_rhs[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_rhs],
+ dist_r2);
+ }
+
+ if(inout_index_particles[(*neighbors[idx_neighbor])[idx_2].first+idx_p2] == 132){// TODO
+ printf("test interaction between :\n");
+ printf("index %ld (%ld) pos %e %e %e\n",
+ (idxPart+idx_p1), -1,
+ descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions + IDX_X],
+ descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions + IDX_Y],
+ descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions + IDX_Z]);
+ printf("index %ld (%ld) pos %e %e %e\n",
+ ((*neighbors[idx_neighbor])[idx_2].first+idx_p2),
+ inout_index_particles[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)],
+ particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_X],
+ particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Y],
+ particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Z]);
+ printf("Radius = %e (%e)\n", sqrt(dist_r2), dist_r2);
}
}
}
@@ -591,16 +672,18 @@ public:
if(releasedAction.first == RELEASE_BUFFER_PARTICLES){
NeighborDescriptor& descriptor = neigDescriptors[releasedAction.second];
assert(descriptor.toCompute != nullptr);
+ assert(descriptor.isRecv);
descriptor.toCompute.release();
}
//////////////////////////////////////////////////////////////////////
/// Merge
//////////////////////////////////////////////////////////////////////
- if(releasedAction.first == MERGE_PARTICLES && more_than_one_thread == false){
+ if(releasedAction.first == MERGE_PARTICLES){
NeighborDescriptor& descriptor = neigDescriptors[releasedAction.second];
- assert(descriptor.isRecv);
+ assert(descriptor.isRecv == false);
assert(descriptor.toRecvAndMerge != nullptr);
- // TODO in_computer.template reduce_particles_rhs<size_particle_rhs>(&particles_current_rhs[0], descriptor.toRecvAndMerge.get(), descriptor.nbParticlesToExchange);
+ in_computer.template reduce_particles_rhs<size_particle_rhs>(&particles_current_rhs[particles_offset_layers[my_nb_cell_levels-descriptor.nbLevelsToExchange]*size_particle_rhs],
+ descriptor.toRecvAndMerge.get(), descriptor.nbParticlesToExchange);
descriptor.toRecvAndMerge.release();
}
}
@@ -632,24 +715,6 @@ public:
&particles_current_rhs[(intervals[idx_1].first+idx_p2)*size_particle_rhs],
dist_r2);
}
- if((inout_index_particles[(intervals[idx_1].first+idx_p1)] == 3
- && inout_index_particles[(intervals[idx_1].first+idx_p2)] == 52)
- || (inout_index_particles[(intervals[idx_1].first+idx_p1)] == 52
- && inout_index_particles[(intervals[idx_1].first+idx_p2)] == 3)){// TODO
- printf("test interaction between :\n");
- printf("index %ld (%ld) pos %e %e %e\n",
- (intervals[idx_1].first+idx_p1), inout_index_particles[(intervals[idx_1].first+idx_p1)],
- particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_X],
- particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Y],
- particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Z]);
- printf("index %ld (%ld) pos %e %e %e\n",
- (intervals[idx_1].first+idx_p2),
- inout_index_particles[(intervals[idx_1].first+idx_p2)],
- particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions + IDX_X],
- particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions + IDX_Y],
- particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions + IDX_Z]);
- printf("Radius = %e (%e)\n", sqrt(dist_r2), dist_r2);
- }
}
}
@@ -671,25 +736,6 @@ public:
&particles_current_rhs[(intervals[idx_2].first+idx_p2)*size_particle_rhs],
dist_r2);
}
-
- if((inout_index_particles[(intervals[idx_1].first+idx_p1)] == 3
- && inout_index_particles[(intervals[idx_2].first+idx_p2)] == 52)
- || (inout_index_particles[(intervals[idx_1].first+idx_p1)] == 52
- && inout_index_particles[(intervals[idx_2].first+idx_p2)] == 3)){// TODO
- printf("interaction between :\n");
- printf("index %ld (%ld) pos %e %e %e\n",
- (intervals[idx_1].first+idx_p1), inout_index_particles[(intervals[idx_1].first+idx_p1)],
- particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_X],
- particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Y],
- particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Z]);
- printf("index %ld (%ld) pos %e %e %e\n",
- (intervals[idx_2].first+idx_p2),
- inout_index_particles[(intervals[idx_2].first+idx_p2)],
- particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions + IDX_X],
- particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions + IDX_Y],
- particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions + IDX_Z]);
- printf("Radius = %e (%e)\n", sqrt(dist_r2), dist_r2);
- }
}
}
}
@@ -701,10 +747,6 @@ public:
long int neighbors_indexes[27];
const int nbNeighbors = my_tree.getNeighbors(currenct_cell_idx, neighbors, neighbors_indexes, false);
- if(iter_cell.first == 3669){ // TODO
- printf("Box %ld has %d neighbors\n", iter_cell.first, nbNeighbors);
- }
-
for(size_t idx_1 = 0 ; idx_1 < intervals.size() ; ++idx_1){
// with other interval
for(size_t idx_neighbor = 0 ; idx_neighbor < nbNeighbors ; ++idx_neighbor){
@@ -726,24 +768,6 @@ public:
&particles_current_rhs[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_rhs],
dist_r2);
}
- if((inout_index_particles[(intervals[idx_1].first+idx_p1)] == 3
- && inout_index_particles[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)] == 52)
- || (inout_index_particles[(intervals[idx_1].first+idx_p1)] == 52
- && inout_index_particles[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)] == 3)){// TODO
- printf("interaction between :\n");
- printf("index %ld (%ld) pos %e %e %e\n",
- (intervals[idx_1].first+idx_p1), inout_index_particles[(intervals[idx_1].first+idx_p1)],
- particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_X],
- particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Y],
- particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Z]);
- printf("index %ld (%ld) pos %e %e %e\n",
- ((*neighbors[idx_neighbor])[idx_2].first+idx_p2),
- inout_index_particles[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)],
- particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_X],
- particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Y],
- particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Z]);
- printf("Radius = %e (%e)\n", sqrt(dist_r2), dist_r2);
- }
}
}
}
diff --git a/bfps/cpp/particles/p2p_tree.hpp b/bfps/cpp/particles/p2p_tree.hpp
index 4715dfad..9f146644 100644
--- a/bfps/cpp/particles/p2p_tree.hpp
+++ b/bfps/cpp/particles/p2p_tree.hpp
@@ -82,7 +82,6 @@ public:
neigh_z_pbc -= nb_cell_levels[IDX_Z];
}
- // Not the current cell
if(include_target || neigh_x_pbc != idx_x || neigh_y_pbc != idx_y || neigh_z_pbc != idx_z){
const long int idx_neigh = get_cell_idx(neigh_x_pbc,
neigh_y_pbc,
--
GitLab