Merge branch 'bugfix/issue38' into bugfix/particle_subsample

cpp/particles/p2p/p2p_distr_mpi.hpp

@@ -96,9 +96,9 @@ protected:
     std::array<real_number,3> spatial_box_width;
     std::array<real_number,3> spatial_box_offset;
 
-    const real_number cutoff_radius_compute;
+    const long double cutoff_radius_compute;
     const int nb_cells_factor;
-    const real_number cutoff_radius;
+    const long double cutoff_radius;
     std::array<long int,3> nb_cell_levels;
 
     template <class DataType, int sizeElement>
@@ -291,6 +291,7 @@ public:
 
         // Some processes might not be involved
         if(nb_processes_involved <= my_rank){
+            DEBUG_MSG("warning: nb_processes_involved <= my_rank, and this process is exiting p2p_distr_mpi::compute_distr now.\nHowever, there is a check below which calls an MPI_Gather over MPI_COMM_WORLD.\n");
             return;
         }
 
@@ -492,10 +493,14 @@ public:
             if(descriptor.isRecv == false){
                 whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
                 mpiRequests.emplace_back();
-                AssertMpi(MPI_Isend(const_cast<partsize_t*>(&descriptor.nbParticlesToExchange),
-                                    1, particles_utils::GetMpiType(partsize_t()),
-                                    descriptor.destProc, TAG_NB_PARTICLES,
-                                    current_com, &mpiRequests.back()));
+                AssertMpi(MPI_Isend(
+                            const_cast<partsize_t*>(&descriptor.nbParticlesToExchange),
+                            1,
+                            particles_utils::GetMpiType(partsize_t()),
+                            descriptor.destProc,
+                            TAG_NB_PARTICLES,
+                            current_com,
+                            &mpiRequests.back()));
 #ifndef NDEBUG // Just for assertion
                 willsend[descriptor.destProc] += 1;
 #endif
@@ -570,21 +575,42 @@ public:
                            MPI_INT,
                            0,
                            MPI_COMM_WORLD);
-                MPI_Gather(willsend.data(), nb_processes_involved, MPI_INT, willsendall.data(),
-                            nb_processes_involved, MPI_INT, 0, MPI_COMM_WORLD);
+                MPI_Gather(willsend.data(),
+                           nb_processes_involved,
+                           MPI_INT,
+                           willsendall.data(),
+                           nb_processes_involved,
+                           MPI_INT,
+                           0,
+                           MPI_COMM_WORLD);
 
                 for(int idxproc = 0 ; idxproc < nb_processes_involved ; ++idxproc){
                     for(int idxtest = 0 ; idxtest < nb_processes_involved ; ++idxtest){
+                        DEBUG_MSG("p2p_distr_mpi::compute_distribution, comparing willsendall(%d, %d)=%d with willrecvall(%d, %d) = %d\n",
+                                idxproc, idxtest, willsendall[idxproc*nb_processes_involved + idxtest],
+                                idxtest, idxproc, willrecvall[idxtest*nb_processes_involved + idxproc]);
                         assert(willsendall[idxproc*nb_processes_involved + idxtest]
                                 == willrecvall[idxtest*nb_processes_involved + idxproc]);
                     }
                 }
             }
             else{
-                MPI_Gather(willrecv.data(), nb_processes_involved, MPI_INT, nullptr,
-                            0, MPI_INT, 0, MPI_COMM_WORLD);
-                MPI_Gather(willsend.data(), nb_processes_involved, MPI_INT, nullptr,
-                            0, MPI_INT, 0, MPI_COMM_WORLD);
+                MPI_Gather(willrecv.data(),
+                           nb_processes_involved,
+                           MPI_INT,
+                           nullptr,
+                           0,
+                           MPI_INT,
+                           0,
+                           MPI_COMM_WORLD);
+                MPI_Gather(willsend.data(),
+                           nb_processes_involved,
+                           MPI_INT,
+                           nullptr,
+                           0,
+                           MPI_INT,
+                           0,
+                           MPI_COMM_WORLD);
             }
         }
 #endif

cpp/particles/particles_utils.hpp

@@ -400,11 +400,11 @@ std::vector<real_number> BuildLimitsAllProcesses(
 
 template <typename partsize_t, typename rnumber, int size_of_particle>
 int set_particle_data_to_zero(
-        rnumber __restrict__ *data,
+        rnumber *data, // TODO: this should be "restrict", but intel can't handle it.
         const partsize_t numberParticles)
 {
     // TODO: ensure simd.
-    // don't use openmp here, as this function may be called from within openmp parallel regions
+    // don't use openmp here, as this function WILL be called from within openmp parallel regions
     std::fill_n(
             data,
             numberParticles*size_of_particle,