Merge branch 'bugfix/issue38' into develop

TurTLE/test/test_Heun_p2p.py

@@ -95,7 +95,7 @@ class ADNS(TurTLE.DNS):
         self.parameters['tracers0_enable_p2p'] = int(0)
         self.parameters['tracers0_enable_inner'] = int(0)
         self.parameters['tracers0_enable_vorticity_omega'] = int(0)
-        self.parameters['tracers0_cutoff'] = float(0.314)
+        self.parameters['tracers0_cutoff'] = float(0.920272)
         self.parameters['tracers0_inner_v0'] = float(1)
         self.parameters['tracers0_lambda'] = float(1)
 

cpp/particles/lock_free_bool_array.hpp

@@ -27,31 +27,61 @@
 #define LOCK_FREE_BOOL_ARRAY_HPP
 
 #include <vector>
-#include <memory>
+#include <atomic>
+#include <unistd.h>
+#include <cstdio>
+#include <omp.h>
+
 
 class lock_free_bool_array{
-    std::vector<std::unique_ptr<long int>> keys;
+    static const int Available = 0;
+    static const int Busy = 1;
+    static const int NoOwner = -1;
+
+    struct Locker {
+        Locker() : lock(Available), ownerId(NoOwner), counter(0) {}
+        std::atomic_int lock;
+        std::atomic_int ownerId;
+        int counter;
+    };
+
+    std::vector<std::unique_ptr<Locker>> keys;
+
 
 public:
-    explicit lock_free_bool_array(const long int inNbKeys = 512){
+    explicit lock_free_bool_array(const long int inNbKeys = 1024){
         keys.resize(inNbKeys);
-        for(std::unique_ptr<long int>& k : keys){
-            k.reset(new long int(0));
+        for(auto& k : keys){
+            k.reset(new Locker());
         }
     }
 
     void lock(const long int inKey){
-        volatile long int* k = keys[inKey%keys.size()].get();
-        long int res = 1;
-        while(res == 1){
-            res = __sync_val_compare_and_swap(k, 0, res);
+        Locker* k = keys[inKey%keys.size()].get();
+        if(k->ownerId.load() != omp_get_thread_num()){
+            int expected = Available;
+            while(!std::atomic_compare_exchange_strong(&k->lock, &expected, Busy)){
+                usleep(1);
+            }
+            k->ownerId.store(omp_get_thread_num());
+            k->counter = 0; // must remain
         }
+        k->counter += 1;
+        assert(k->lock.load() == Busy);
+        assert(k->counter >= 1);
+        assert(k->ownerId.load() == omp_get_thread_num());
     }
 
     void unlock(const long int inKey){
-        volatile long int* k = keys[inKey%keys.size()].get();
-        assert(k && *k);
-        (*k) = 0;
+        Locker* k = keys[inKey%keys.size()].get();
+        assert(k->lock.load() == Busy);
+        assert(k->counter >= 1);
+        assert(k->ownerId.load() == omp_get_thread_num());
+        k->counter -= 1;
+        if(k->counter == 0){
+            k->ownerId.store(NoOwner);
+            k->lock.store(Available);
+        }
     }
 };
 

cpp/particles/p2p/p2p_distr_mpi.hpp

@@ -138,7 +138,8 @@ protected:
 
     static int foundGridFactor(const real_number in_cutoff_radius, const std::array<real_number,3>& in_spatial_box_width){
         int idx_factor = 1;
-        while(in_cutoff_radius <= in_spatial_box_width[IDXC_Z]/real_number(idx_factor+1)){
+        DEBUG_MSG("in_cutoff_radius is %g, in_spatial_box_width[IDXC_Z] = %g\n", in_cutoff_radius, in_spatial_box_width[IDXC_Z]);
+        while(in_cutoff_radius <= in_spatial_box_width[IDXC_Z]/(long double)(idx_factor+1)){
             idx_factor += 1;
         }
         return idx_factor;
@@ -227,11 +228,10 @@ public:
     }
 
     long int last_cell_level_proc(const int dest_proc) const{
-        const real_number field_section_width_z = spatial_box_width[IDXC_Z]/real_number(field_grid_dim[IDXC_Z]);
-        const long int limite = static_cast<long int>((field_section_width_z*real_number(partition_interval_offset_per_proc[dest_proc+1])
-                                     - std::numeric_limits<real_number>::epsilon())/cutoff_radius);
-        if(static_cast<real_number>(limite)*cutoff_radius
-                == field_section_width_z*real_number(partition_interval_offset_per_proc[dest_proc+1])){
+        const long double field_section_width_z = spatial_box_width[IDXC_Z]/(long double)(field_grid_dim[IDXC_Z]);
+        const long int limite = static_cast<long int>((field_section_width_z*(long double)(partition_interval_offset_per_proc[dest_proc+1])
+                                     - std::numeric_limits<long double>::epsilon())/cutoff_radius);
+        if(limite == nb_cell_levels[IDXC_Z]){
             return limite-1;
         }
         return limite;
@@ -296,7 +296,9 @@ public:
         }
 
         const long int my_top_z_cell_level = last_cell_level_proc(my_rank);
+        assert(my_top_z_cell_level < nb_cell_levels[IDXC_Z]);
         const long int my_down_z_cell_level = first_cell_level_proc(my_rank);
+        assert(0 <= my_down_z_cell_level);
         const long int my_nb_cell_levels = 1+my_top_z_cell_level-my_down_z_cell_level;
 
         current_offset_particles_for_partition[0] = 0;
@@ -426,6 +428,18 @@ public:
         assert(mpiRequests.size() == 0);
         neigDescriptors.clear();
 
+        DEBUG_MSG("my_top_z_cell_level = %d\n", int(my_top_z_cell_level));
+        DEBUG_MSG("my_down_z_cell_level = %d\n", int(my_down_z_cell_level));
+        DEBUG_MSG("first_cell_level_proc(0) = %d\n", int(first_cell_level_proc(0)));
+        DEBUG_MSG("first_cell_level_proc(1) = %d\n", int(first_cell_level_proc(1)));
+        DEBUG_MSG("first_cell_level_proc(2) = %d\n", int(first_cell_level_proc(2)));
+        DEBUG_MSG("first_cell_level_proc(3) = %d\n", int(first_cell_level_proc(3)));
+        DEBUG_MSG("nb_cell_levels[IDXC_Z] = %d\n", int(nb_cell_levels[IDXC_Z]));
+        DEBUG_MSG("last_cell_level_proc(0) = %d\n", int(last_cell_level_proc(0)));
+        DEBUG_MSG("last_cell_level_proc(1) = %d\n", int(last_cell_level_proc(1)));
+        DEBUG_MSG("last_cell_level_proc(2) = %d\n", int(last_cell_level_proc(2)));
+        DEBUG_MSG("last_cell_level_proc(3) = %d\n", int(last_cell_level_proc(3)));
+
         // Find process with at least one neighbor
         {
             int dest_proc = (my_rank+1)%nb_processes_involved;
@@ -438,6 +452,7 @@ public:
                         && (my_top_z_cell_level-1+2)%nb_cell_levels[IDXC_Z] <= last_cell_level_proc(dest_proc)){
                     nb_levels_to_send += 1;
                 }
+                DEBUG_MSG("looking at dest_proc = %d ; nb_levels_to_send = %d\n", dest_proc, nb_levels_to_send);
 
                 NeighborDescriptor descriptor;
                 descriptor.destProc = dest_proc;
@@ -461,6 +476,7 @@ public:
                         && first_cell_level_proc(src_proc) <= (my_down_z_cell_level-1+2)%nb_cell_levels[IDXC_Z]){
                     nb_levels_to_recv += 1;
                 }
+                DEBUG_MSG("looking at src_proc = %d ; nb_levels_to_recv = %d\n", src_proc, nb_levels_to_recv);
 
                 NeighborDescriptor descriptor;
                 descriptor.destProc = src_proc;
@@ -475,6 +491,8 @@ public:
             }
         }
 
+        DEBUG_MSG_WAIT(current_com, "hello, found processes with neighbours.\n");
+
         //////////////////////////////////////////////////////////////////////
         /// Exchange the number of particles in each partition
         /// Could involve only here but I do not think it will be a problem