Merge branch 'feature/collisions' into develop

CMakeLists.txt

@@ -364,6 +364,7 @@ set(hpp_for_lib
     ${PROJECT_SOURCE_DIR}/cpp/particles/lock_free_bool_array.hpp
     ${PROJECT_SOURCE_DIR}/cpp/particles/p2p_computer_empty.hpp
     ${PROJECT_SOURCE_DIR}/cpp/particles/p2p_computer.hpp
+    ${PROJECT_SOURCE_DIR}/cpp/particles/p2p_ghost_collision_base.hpp
     ${PROJECT_SOURCE_DIR}/cpp/particles/p2p_ghost_collisions.hpp
     ${PROJECT_SOURCE_DIR}/cpp/particles/p2p_cylinder_collisions.hpp
     ${PROJECT_SOURCE_DIR}/cpp/particles/p2p_distr_mpi.hpp

README.rst

@@ -257,7 +257,7 @@ Afterwards, please run variations of the following command:
 
     .. code:: bash
 
-        python ${TURTLE REPOSITORY}/tests/DNS/test_scaling.py D \
+        python ${TURTLE_REPOSITORY}/tests/DNS/test_scaling.py D \
             -n 128 \
             --nprocesses 4 \
             --ncores 1 \

TurTLE/DNS.py

@@ -445,6 +445,7 @@ class DNS(_code):
                 'NSVE_Stokes_particles',
                 'NSVEparticles',
                 'static_field',
+                'static_field_with_ghost_collisions',
                 'kraichnan_field']:
             assert (self.parameters['niter_todo'] % self.parameters['niter_part'] == 0)
             assert (self.parameters['niter_out']  % self.parameters['niter_part'] == 0)
@@ -488,6 +489,8 @@ class DNS(_code):
                                                  4),
                                      dtype = np.int64)
             ofile['checkpoint'] = int(0)
+            if self.dns_type in ['static_field_with_ghost_collisions']:
+                ofile.create_group('statistics/collisions')
         if (self.dns_type in ['NSVE', 'NSVE_no_output']):
             return None
         return None
@@ -646,6 +649,10 @@ class DNS(_code):
                 'static_field',
                 help = 'static field with basic fluid tracers')
 
+        parser_static_field_with_ghost_collisions = subparsers.add_parser(
+                'static_field_with_ghost_collisions',
+                help = 'static field with basic fluid tracers and ghost collisions')
+
         parser_kraichnan_field = subparsers.add_parser(
                 'kraichnan_field',
                 help = 'Kraichnan field with basic fluid tracers')
@@ -673,6 +680,7 @@ class DNS(_code):
                 'NSVE_Stokes_particles',
                 'NSVEp_extra',
                 'static_field',
+                'static_field_with_ghost_collisions',
                 'kraichnan_field']:
             eval('self.simulation_parser_arguments({0})'.format('parser_' + pp))
             eval('self.job_parser_arguments({0})'.format('parser_' + pp))
@@ -749,6 +757,7 @@ class DNS(_code):
                 'NSVEparticles_no_output',
                 'NSVEp_extra_sampling',
                 'static_field',
+                'static_field_with_ghost_collisions',
                 'kraichnan_field']:
             for k in self.NSVEp_extra_parameters.keys():
                 self.parameters[k] = self.NSVEp_extra_parameters[k]
@@ -823,6 +832,7 @@ class DNS(_code):
             if self.dns_type in [
                     'kraichnan_field',
                     'static_field',
+                    'static_field_with_ghost_collisions',
                     'NSVEparticles',
                     'NSVEcomplex_particles',
                     'NSVE_Stokes_particles',
@@ -1138,6 +1148,7 @@ class DNS(_code):
         if self.dns_type in [
                 'kraichnan_field',
                 'static_field',
+                'static_field_with_ghost_collisions',
                 'NSVEparticles',
                 'NSVEcomplex_particles',
                 'NSVE_Stokes_particles',

cpp/base.hpp

@@ -43,6 +43,11 @@ inline int MOD(int a, int n)
     return ((a%n) + n) % n;
 }
 
+template <typename T> int sgn(T val)
+{
+        return (T(0) < val) - (val < T(0));
+}
+
 /////////////////////////////////////////////////////////////
 /////////////////////////////////////////////////////////////
 

cpp/full_code/static_field.cpp

@@ -64,19 +64,24 @@ int static_field<rnumber>::initialize(void)
     this->vorticity = new field<rnumber, FFTW, THREE>(
             nx, ny, nz,
             this->comm,
-	    fftw_planner_string_to_flag[this->fftw_plan_rigor]);
+            fftw_planner_string_to_flag[this->fftw_plan_rigor]);
     this->vorticity->real_space_representation = false;
 
     this->velocity = new field<rnumber, FFTW, THREE>(
             nx, ny, nz,
             this->comm,
-	    fftw_planner_string_to_flag[this->fftw_plan_rigor]);
+            fftw_planner_string_to_flag[this->fftw_plan_rigor]);
     this->velocity->real_space_representation = false;
-    //read vorticity field
+    //read static vorticity field from iteration 0
+    std::string checkpoint_fname = (
+            std::string(this->simname) +
+            std::string("_checkpoint_") +
+            std::to_string(0) +
+            std::string(".h5"));
     this->vorticity->io(
-	    this->get_current_fname(),
+	    checkpoint_fname,
 	    "vorticity",
-	    this->iteration,
+	    0,
 	    true);
     this->kk = new kspace<FFTW, SMOOTH>(
             this->vorticity->clayout, this->dkx, this->dky, this->dkz);
@@ -96,7 +101,7 @@ int static_field<rnumber>::initialize(void)
                 this->kk,                     // (kspace object, contains dkx, dky, dkz)
                 tracers0_integration_steps, // to check coherency between parameters and hdf input file (nb rhs)
                 (long long int)nparticles,  // to check coherency between parameters and hdf input file
-	            this->get_current_fname(),    // particles input filename
+                this->get_current_fname(),    // particles input filename
                 std::string("/tracers0/state/") + std::to_string(this->iteration), // dataset name for initial input
                 std::string("/tracers0/rhs/")  + std::to_string(this->iteration),  // dataset name for initial input
                 tracers0_neighbours,        // parameter (interpolation no neighbours)
@@ -144,6 +149,7 @@ int static_field<rnumber>::write_checkpoint(void)
             this->ps->getLocalNbParticles(),
             this->iteration);
     this->particles_output_writer_mpi->close_file();
+    this->write_iteration();
     return EXIT_SUCCESS;
 }
 

cpp/hdf5_tools.cpp

@@ -3,6 +3,63 @@
 #include <cfloat>
 #include <climits>
 
+template <> hid_t hdf5_tools::hdf5_type_id<char>()
+{
+    return H5T_NATIVE_CHAR;
+}
+template <> hid_t hdf5_tools::hdf5_type_id<signed char>()
+{
+    return H5T_NATIVE_SCHAR;
+}
+template <> hid_t hdf5_tools::hdf5_type_id<unsigned char>()
+{
+    return H5T_NATIVE_UCHAR;
+}
+template <> hid_t hdf5_tools::hdf5_type_id<short>()
+{
+    return H5T_NATIVE_SHORT;
+}
+template <> hid_t hdf5_tools::hdf5_type_id<unsigned short>()
+{
+    return H5T_NATIVE_USHORT;
+}
+template <> hid_t hdf5_tools::hdf5_type_id<int>()
+{
+    return H5T_NATIVE_INT;
+}
+template <> hid_t hdf5_tools::hdf5_type_id<unsigned>()
+{
+    return H5T_NATIVE_UINT;
+}
+template <> hid_t hdf5_tools::hdf5_type_id<long>()
+{
+    return H5T_NATIVE_LONG;
+}
+template <> hid_t hdf5_tools::hdf5_type_id<unsigned long>()
+{
+    return H5T_NATIVE_ULONG;
+}
+template <> hid_t hdf5_tools::hdf5_type_id<long long>()
+{
+    return H5T_NATIVE_LLONG;
+}
+template <> hid_t hdf5_tools::hdf5_type_id<unsigned long long>()
+{
+    return H5T_NATIVE_ULLONG;
+}
+template <> hid_t hdf5_tools::hdf5_type_id<float>()
+{
+    return H5T_NATIVE_FLOAT;
+}
+template <> hid_t hdf5_tools::hdf5_type_id<double>()
+{
+    return H5T_NATIVE_DOUBLE;
+}
+template <> hid_t hdf5_tools::hdf5_type_id<long double>()
+{
+    return H5T_NATIVE_LDOUBLE;
+}
+
 int hdf5_tools::require_size_single_dataset(hid_t dset, int tsize)
 {
     TIMEZONE("hdf5_tools::require_size_single_dataset");
@@ -178,6 +235,36 @@ number hdf5_tools::read_value(
     return result;
 }
 
+template <typename number>
+int hdf5_tools::write_value_with_single_rank(
+        const hid_t group,
+        const std::string dset_name,
+        const number value)
+{
+    hid_t dset, mem_dtype;
+    if (typeid(number) == typeid(int))
+        mem_dtype = H5Tcopy(H5T_NATIVE_INT);
+    else if (typeid(number) == typeid(double))
+        mem_dtype = H5Tcopy(H5T_NATIVE_DOUBLE);
+    if (H5Lexists(group, dset_name.c_str(), H5P_DEFAULT))
+    {
+        dset = H5Dopen(group, dset_name.c_str(), H5P_DEFAULT);
+    }
+    else
+    {
+        hid_t fspace;
+        hsize_t count[1];
+        count[0] = 1;
+        fspace = H5Screate_simple(1, count, NULL);
+        dset = H5Dcreate(group, dset_name.c_str(), mem_dtype, fspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+        H5Sclose(fspace);
+    }
+    H5Dwrite(dset, mem_dtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &value);
+    H5Dclose(dset);
+    H5Tclose(mem_dtype);
+    return EXIT_SUCCESS;
+}
+
 template <typename dtype>
 std::vector<dtype> hdf5_tools::read_vector_with_single_rank(
         const int myrank,
@@ -244,6 +331,44 @@ std::string hdf5_tools::read_string(
     }
 }
 
+template <class partsize_t>
+int hdf5_tools::write_particle_ID_pairs_with_single_rank(
+        const std::vector<partsize_t> v,
+        const hid_t group,
+        const std::string dset_name)
+{
+    TIMEZONE("hdf5_tools::write_particle_ID_pairs_with_single_rank");
+    // the vector contains pair information, so its size must be a multiple of 2
+    assert((v.size() % 2) == 0);
+    // file space creation
+    hid_t fspace;
+    hsize_t dims[2];
+    dims[0] = v.size()/2;
+    dims[1] = 2;
+    fspace = H5Screate_simple(2, dims, NULL);
+    // create dataset
+    hsize_t dset_id = H5Dcreate(
+            group,
+            dset_name.c_str(),
+            hdf5_tools::hdf5_type_id<partsize_t>(),
+            fspace,
+            H5P_DEFAULT,
+            H5P_DEFAULT,
+            H5P_DEFAULT);
+    // write data
+    H5Dwrite(
+            dset_id,
+            hdf5_tools::hdf5_type_id<partsize_t>(),
+            H5S_ALL,
+            H5S_ALL,
+            H5P_DEFAULT,
+            &v.front());
+    // clean up
+    H5Dclose(dset_id);
+    H5Sclose(fspace);
+    return EXIT_SUCCESS;
+}
+
 template
 std::vector<int> hdf5_tools::read_vector<int>(
         const hid_t,
@@ -280,3 +405,21 @@ double hdf5_tools::read_value<double>(
         const hid_t,
         const std::string);
 
+template
+int hdf5_tools::write_value_with_single_rank<int>(
+        const hid_t group,
+        const std::string dset_name,
+        const int value);
+
+template
+int hdf5_tools::write_value_with_single_rank<double>(
+        const hid_t group,
+        const std::string dset_name,
+        const double value);
+
+template
+int hdf5_tools::write_particle_ID_pairs_with_single_rank(
+        const std::vector<long long> v,
+        const hid_t group,
+        const std::string dset_name);
+

cpp/hdf5_tools.hpp

@@ -33,6 +33,9 @@
 
 namespace hdf5_tools
 {
+    // see https://support.hdfgroup.org/HDF5/doc/H5.user/Datatypes.html
+    template <typename data_type> hid_t hdf5_type_id();
+
     int grow_single_dataset(
             hid_t dset,
             int tincrement);
@@ -84,6 +87,18 @@ namespace hdf5_tools
     number read_value(
             const hid_t group,
             const std::string dset_name);
+
+    template <typename number>
+    int write_value_with_single_rank(
+            const hid_t group,
+            const std::string dset_name,
+            const number value);
+
+    template <class partsize_t>
+    int write_particle_ID_pairs_with_single_rank(
+            const std::vector<partsize_t> v,
+            const hid_t group,
+            const std::string dset_name);
 }
 
 #endif//HDF5_TOOLS_HPP

cpp/particles/abstract_particles_system.hpp

@@ -27,6 +27,7 @@
 #define ABSTRACT_PARTICLES_SYSTEM_HPP
 
 #include <memory>
+#include <vector>
 
 //- Not generic to enable sampling begin
 #include "field.hpp"
@@ -39,6 +40,8 @@ class abstract_particles_system {
 public:
     virtual ~abstract_particles_system(){}
 
+    virtual void delete_particles_from_indexes(const std::vector<partsize_t>& inIndexToDelete) = 0;
+
     virtual void compute() = 0;
 
     virtual void compute_p2p() = 0;

cpp/particles/p2p_computer.hpp

@@ -70,7 +70,9 @@ public:
     }
 
     template <int size_particle_positions, int size_particle_rhs>
-    void compute_interaction(const real_number pos_part1[], real_number rhs_part1[],
+    void compute_interaction(const partsize_t /*idx_part1*/,
+                             const real_number pos_part1[], real_number rhs_part1[],
+                             const partsize_t /*idx_part2*/,
                              const real_number pos_part2[], real_number rhs_part2[],
                              const real_number dist_pow2, const real_number cutoff,
                              const real_number /*xshift_coef*/, const real_number /*yshift_coef*/, const real_number /*zshift_coef*/) const{

cpp/particles/p2p_computer_empty.hpp

@@ -40,7 +40,9 @@ public:
     }
 
     template <int size_particle_positions, int size_particle_rhs>
-    void compute_interaction(const real_number /*pos_part1*/[], real_number /*rhs_part1*/[],
+    void compute_interaction(const partsize_t /*idx_part1*/,
+                             const real_number /*pos_part1*/[], real_number /*rhs_part1*/[],
+                             const partsize_t /*idx_part2*/,
                              const real_number /*pos_part2*/[], real_number /*rhs_part2*/[],
                              const real_number /*dist_pow2*/,  const real_number /*cutoff*/,
                              const real_number /*xshift_coef*/, const real_number /*yshift_coef*/, const real_number /*zshift_coef*/) const{

cpp/particles/p2p_cylinder_collisions.hpp

@@ -24,54 +24,138 @@
 #define P2P_CYLINDER_GHOST_COLLISIONS_HPP
 
 #include <cstring>
+#include <math.h>       /* for sqrt, abs */
+
+#include <set>
+#include <utility>
+#include <vector>
+#include "particles/p2p_ghost_collision_base.hpp"
+
 
 template <class real_number, class partsize_t>
-class p2p_cylinder_ghost_collisions{
-    long int collision_counter;
+class p2p_cylinder_ghost_collisions: public p2p_ghost_collision_base<real_number, partsize_t>
+{
+    private:
+        // Following doubles are needed for the collision computation
+        double width;
+        double length;
 
 public:
-    p2p_cylinder_ghost_collisions() : collision_counter(0){}
+    p2p_cylinder_ghost_collisions() : width(1.),length(1.){}
 
     // Copy constructor use a counter set to zero
-    p2p_cylinder_ghost_collisions(const p2p_cylinder_ghost_collisions&) : collision_counter(0){}
-
-    template <int size_particle_rhs>
-    void init_result_array(real_number /*rhs*/[], const partsize_t /*nbParticles*/) const{
-    }
-
-    template <int size_particle_rhs>
-    void reduce_particles_rhs(real_number /*rhs_dst*/[], const real_number /*rhs_src*/[], const partsize_t /*nbParticles*/) const{
-    }
+    p2p_cylinder_ghost_collisions(const p2p_cylinder_ghost_collisions&){}
 
     template <int size_particle_positions, int size_particle_rhs>
-    void compute_interaction(const real_number /*pos_part1*/[], real_number /*rhs_part1*/[],
-                             const real_number /*pos_part2*/[], real_number /*rhs_part2*/[],
+    void compute_interaction(const partsize_t idx_part1,
+                             const real_number pos_part1[], real_number /*rhs_part1*/[],
+                             const partsize_t idx_part2,
+                             const real_number pos_part2[], real_number /*rhs_part2*/[],
                              const real_number /*dist_pow2*/,  const real_number /*cutoff*/,
                              const real_number /*xshift_coef*/, const real_number /*yshift_coef*/, const real_number /*zshift_coef*/){
-        bool cylinders_overlap = false;
+        double x, y, z, pq, xp, xq, t, s, x_dist, y_dist, z_dist ,min_distance, pi, pi2;
 
-        /* TODO: check if cylinders overlap or not, set `cylinders_overlap` accordingly */
+        pi2 = atan(1.0)*8.0;
+        pi = atan(1.0)*4.0;
 
-        if(cylinders_overlap)
-            collision_counter += 1;
+        /* Relative position vector of the two particles (x,y,z)^T */
+        /* complicated usage of fmod, fabs and sgn because C++ fmod does not do what it should. */
+        x = std::fmod(pos_part2[0],pi2)-std::fmod(pos_part1[0],pi2);
+        y = std::fmod(pos_part2[1],pi2)-std::fmod(pos_part1[1],pi2);
+        z = std::fmod(pos_part2[2],pi2)-std::fmod(pos_part1[2],pi2);
+
+        x = ( std::fmod( std::fabs(x) +pi ,pi2) - pi ) * sgn(x) ;
+        y = ( std::fmod( std::fabs(y) +pi ,pi2) - pi ) * sgn(y) ;
+        z = ( std::fmod( std::fabs(z) +pi ,pi2) - pi ) * sgn(z) ;
+
+        if( sqrt(x*x+y*y+z*z) <= length )
+        {
+
+        /* p and q are the orientation vectors of the first and second particles. */
+        /* pq, xp, xq are scalar products of these vectors with x, relative position */
+        pq = pos_part1[3] * pos_part2[3] + pos_part1[4] * pos_part2[4] + pos_part1[5] * pos_part2[5];
+        xp = x * pos_part1[3] + y * pos_part1[4] + z * pos_part1[5];
+        xq = x * pos_part2[3] + y * pos_part2[4] + z * pos_part2[5];
+        /* t and s parametrize the two rods. Find min distance: */
+        assert(this->length > 0);
+        t = 2.0/(this->length*(pq*pq-1.0))*(-xp+pq*xq);
+        s = 2.0/(this->length*(pq*pq-1.0))*(-pq*xp+xq);
+        /* Test if -1<s<1 and -1<t<1 */
+        if( abs(t)<=1.0 and abs(s)<=1.0 ) {
+            /* Get minimal distance in case of both t and s in {-1,1}. Else: check edges */
+            x_dist = this->length*0.5*t*pos_part1[3]-this->length*0.5*s*pos_part2[3]-x;
+            y_dist = this->length*0.5*t*pos_part1[4]-this->length*0.5*s*pos_part2[4]-y;
+            z_dist = this->length*0.5*t*pos_part1[5]-this->length*0.5*s*pos_part2[5]-z;
+            min_distance = sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist);
+        }
+        else
+        {
+            min_distance = this->length;
+            /* t fixed at 1, find min along s */
+            t = 1.0;
+            s = t*pq-2.0/this->length*xq;
+            if( abs(s)>1.0 ) { s = s / abs(s) ;}
+            x_dist = this->length*0.5*t*pos_part1[3]-this->length*0.5*s*pos_part2[3]-x;
+            y_dist = this->length*0.5*t*pos_part1[4]-this->length*0.5*s*pos_part2[4]-y;
+            z_dist = this->length*0.5*t*pos_part1[5]-this->length*0.5*s*pos_part2[5]-z;
+            min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance );
+            /* t fixed at -1, find min along s */
+            t = -1.0;
+            s = t*pq-2.0/this->length*xq;
+            if( abs(s)>1.0 ) { s = s / abs(s) ;}
+            x_dist = this->length*0.5*t*pos_part1[3]-this->length*0.5*s*pos_part2[3]-x;
+            y_dist = this->length*0.5*t*pos_part1[4]-this->length*0.5*s*pos_part2[4]-y;
+            z_dist = this->length*0.5*t*pos_part1[5]-this->length*0.5*s*pos_part2[5]-z;
+            min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance );
+            /* s fixed at 1, find min along t */
+            s = 1.0;
+            t = s*pq+2.0/this->length*xp;
+            if( abs(t)>1.0 ) { t = t / abs(t) ;}
+            x_dist = this->length*0.5*t*pos_part1[3]-this->length*0.5*s*pos_part2[3]-x;
+            y_dist = this->length*0.5*t*pos_part1[4]-this->length*0.5*s*pos_part2[4]-y;
+            z_dist = this->length*0.5*t*pos_part1[5]-this->length*0.5*s*pos_part2[5]-z;
+            min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance );
+            /* s fixed at -1, find min along t */
+            s = -1.0;
+            t = s*pq+2.0/this->length*xp;
+            if( abs(t)>1.0 ) { t = t / abs(t) ;}
+            x_dist = this->length*0.5*t*pos_part1[3]-this->length*0.5*s*pos_part2[3]-x;
+            y_dist = this->length*0.5*t*pos_part1[4]-this->length*0.5*s*pos_part2[4]-y;
+            z_dist = this->length*0.5*t*pos_part1[5]-this->length*0.5*s*pos_part2[5]-z;
+            min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance );
+        }
+        /* If cylinders overlap count it as a collision */
+        if( min_distance<=width ){
+            std::pair <partsize_t, partsize_t> single_collision_pair(idx_part1, idx_part2);
+            this->collision_pairs_local.insert(single_collision_pair);
+            //DEBUG_MSG("inside compute interaction idx_part1 = %ld and idx_part2 = %ld\n", idx_part1, idx_part2);
+            assert(idx_part1!=idx_part2);
+        }
+
+		}
     }
 
-    void merge(const p2p_cylinder_ghost_collisions& other){
-        collision_counter += other.collision_counter;
+    void set_width(const double WIDTH)
+    {
+        this->width = WIDTH;
     }
 
-    constexpr static bool isEnable() {
-        return true;
+    double get_width()
+    {
+        return this->width;
     }
 
-    long int get_collision_counter() const{
-        return collision_counter;
+    void set_length(const double LENGTH)
+    {
+        this->length = LENGTH;
     }
 
-    void reset_collision_counter(){
-        collision_counter = 0;
+    double get_length()
+    {
+        return this->length;
     }
 };
 
 
 #endif // P2P_CYLINDER_GHOST_COLLISIONS_HPP
+

cpp/particles/p2p_ghost_collision_base.hpp

+/******************************************************************************
+*                                                                             *
+*  Copyright 2019 Max Planck Institute for Dynamics and Self-Organization     *
+*                                                                             *
+*  This file is part of bfps.                                                 *
+*                                                                             *
+*  bfps is free software: you can redistribute it and/or modify               *
+*  it under the terms of the GNU General Public License as published          *
+*  by the Free Software Foundation, either version 3 of the License,          *
+*  or (at your option) any later version.                                     *
+*                                                                             *
+*  bfps is distributed in the hope that it will be useful,                    *
+*  but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+*  GNU General Public License for more details.                               *
+*                                                                             *
+*  You should have received a copy of the GNU General Public License          *
+*  along with bfps.  If not, see <http://www.gnu.org/licenses/>               *
+*                                                                             *
+* Contact: Cristian.Lalescu@ds.mpg.de                                         *
+*                                                                             *
+******************************************************************************/
+#ifndef P2P_GHOST_COLLISION_BASE_HPP
+#define P2P_GHOST_COLLISION_BASE_HPP
+
+#include <cstring>
+#include <set>
+#include <utility>
+#include <vector>
+template < class partsize_t>
+std::vector<partsize_t> pairs2vec(std::set <std::pair <partsize_t,partsize_t>> ID_pairs){
+  std::vector<partsize_t> v(2*ID_pairs.size());
+  int i = 0;
+  for(auto p: ID_pairs)
+  {
+    v[2*i] = p.first;
+    v[2*i+1] = p.second;
+    i++;
+  }
+  return v;
+}
+
+template < class partsize_t>
+std::set <std::pair <partsize_t,partsize_t>> vec2pairs(std::vector<partsize_t> v){
+  std::set <std::pair <partsize_t,partsize_t>> ID_pairs;
+  assert(v.size()%2 == 0);
+  for(int i=0; i < int(v.size())/2; i++)
+    {
+      //DEBUG_MSG((std::to_string(v[2*i])+" "+std::to_string(v[2*i+1])+"\n").c_str());
+      std::pair <partsize_t, partsize_t> single_collision_pair(v[2*i],v[2*i+1]);
+      ID_pairs.insert(single_collision_pair);
+    }
+  return ID_pairs;
+}
+
+template <class partsize_t>
+void print_pair_vec(std::vector<partsize_t> vec)