Fixes for parallel collision detection

doc/sphinx/advanced_methods.rst

@@ -68,6 +68,12 @@ Several modes are available for different types of binding.
   * `part_type_vs`: Particle type assigned to the virtual site created during the collision. 
   * `distance_glued_particle_to_vs`: Distance of the virtual site to the particle being bound to it (`small` particle).
 
+  Note: When the type of a particle is changed on collision, this makes the 
+  particle inert with regards to further collision. Should a particle  of
+  type `part_type_to_be_glued` collide with two particles in a single 
+  time step, no guarantees are made with regards to which partner is selected.
+  In particular, there is no guarantee that the choice is unbiased.
+
 
 
 

src/core/cells.cpp

@@ -399,3 +399,15 @@ void cells_update_ghosts() {
     /* Communication step: ghost information */
     ghost_communicator(&cell_structure.update_ghost_pos_comm);
 }
+
+Cell *find_current_cell(const Particle &p) {
+  auto c = cell_structure.position_to_cell(p.r.p.data());
+  if (c) {
+    return c;
+  } else if (!p.l.ghost) {
+    // Old pos must lie within the cell system
+    return cell_structure.position_to_cell(p.l.p_old.data());
+  } else {
+    return nullptr;
+  }
+}

src/core/cells.hpp

@@ -175,7 +175,7 @@ struct CellStructure {
       \param  pos Position of a particle.
       \return pointer to cell  where to put the particle.
   */
-  Cell *(*position_to_cell)(double pos[3]);
+  Cell *(*position_to_cell)(const double pos[3]);
 };
 
 /*@}*/
@@ -312,4 +312,13 @@ void local_sort_particles();
 
 /*@}*/
 
+/* @brief Finds the cell in which a particle is stored
+
+   Uses position_to_cell on p.r.p. If this is not on the node's domain,
+   uses position at last Verlet list rebuild (p.l.p_old).
+
+   @return pointer to the cell or nullptr if the particle is not on the node
+*/
+Cell *find_current_cell(const Particle &p);
+
 #endif

src/core/collision.cpp

@@ -28,6 +28,7 @@
 #include "interaction_data.hpp"
 #include "particle_data.hpp"
 #include "rotation.hpp"
+#include "utils/mpi/all_compare.hpp"
 #include "virtual_sites/VirtualSitesRelative.hpp"
 
 #ifdef COLLISION_DETECTION_DEBUG
@@ -55,7 +56,12 @@ Collision_parameters collision_params;
 /** @brief Return true if a bond between the centers of the colliding particles
  * needs to be placed. At this point, all modes need this */
 inline bool bind_centers() {
-  return collision_params.mode != COLLISION_MODE_OFF;
+  // Note that the glue to surface mode adds bonds between the centers
+  // but does so later in the process. This is needed to gaurantee that
+  // a particle can only be glued once, even if queued twice in a single
+  // time step
+  return collision_params.mode != COLLISION_MODE_OFF &&
+         collision_params.mode != COLLISION_MODE_GLUE_TO_SURF;
 }
 
 bool validate_collision_parameters() {
@@ -100,12 +106,6 @@ bool validate_collision_parameters() {
   }
 #endif
 
-  if ((collision_params.mode != COLLISION_MODE_OFF) && (n_nodes > 1)) {
-    runtimeErrorMsg() << "The collision detection schemes are currently not "
-                         "available in parallel simulations";
-    return false;
-  }
-
   // Check if bonded ia exist
   if ((collision_params.mode & COLLISION_MODE_BOND) &&
       (collision_params.bond_centers >= bonded_ia_params.size())) {
@@ -228,30 +228,30 @@ void queue_collision(const int part1, const int part2) {
 
 /** @brief Calculate position of vs for GLUE_TO_SURFACE mode
  *    Reutnrs id of particle to bind vs to */
-int glue_to_surface_calc_vs_pos(const Particle *const p1,
-                                const Particle *const p2, Vector3d &pos) {
-  int bind_vs_to_pid;
+const Particle &glue_to_surface_calc_vs_pos(const Particle &p1,
+                                            const Particle &p2, Vector3d &pos) {
   double vec21[3];
   double c;
-  get_mi_vector(vec21, p1->r.p, p2->r.p);
+  get_mi_vector(vec21, p1.r.p, p2.r.p);
   const double dist_betw_part = sqrt(sqrlen(vec21));
 
   // Find out, which is the particle to be glued.
-  if ((p1->p.type == collision_params.part_type_to_be_glued) &&
-      (p2->p.type == collision_params.part_type_to_attach_vs_to)) {
+  if ((p1.p.type == collision_params.part_type_to_be_glued) &&
+      (p2.p.type == collision_params.part_type_to_attach_vs_to)) {
     c = 1 - collision_params.dist_glued_part_to_vs / dist_betw_part;
-    bind_vs_to_pid = p2->p.identity;
-  } else if ((p2->p.type == collision_params.part_type_to_be_glued) &&
-             (p1->p.type == collision_params.part_type_to_attach_vs_to)) {
+  } else if ((p2.p.type == collision_params.part_type_to_be_glued) &&
+             (p1.p.type == collision_params.part_type_to_attach_vs_to)) {
     c = collision_params.dist_glued_part_to_vs / dist_betw_part;
-    bind_vs_to_pid = p1->p.identity;
   } else {
     throw std::runtime_error("This should never be thrown. Bug.");
   }
   for (int i = 0; i < 3; i++) {
-    pos[i] = p2->r.p[i] + vec21[i] * c;
+    pos[i] = p2.r.p[i] + vec21[i] * c;
   }
-  return bind_vs_to_pid;
+  if (p1.p.type == collision_params.part_type_to_attach_vs_to)
+    return p1;
+  else
+    return p2;
 }
 
 void bind_at_point_of_collision_calc_vs_pos(const Particle *const p1,
@@ -361,6 +361,7 @@ void coldet_do_three_particle_bond(Particle &p, Particle &p1, Particle &p2) {
   bondT[0] = bond_id;
   bondT[1] = p1.p.identity;
   bondT[2] = p2.p.identity;
+
   local_change_bond(p.p.identity, bondT, 0);
 }
 
@@ -392,16 +393,21 @@ void bind_at_poc_create_bond_between_vs(const int current_vs_pid,
     // Create bond between the virtual particles
     bondG[0] = collision_params.bond_vs;
     bondG[1] = current_vs_pid - 2;
-    local_change_bond(current_vs_pid - 1, bondG, 0);
+    // Only add bond if vs was created on this node
+    if (local_particles[current_vs_pid - 1])
+      local_change_bond(current_vs_pid - 1, bondG, 0);
     break;
   }
   case 2: {
     // Create 1st bond between the virtual particles
     bondG[0] = collision_params.bond_vs;
     bondG[1] = c.pp1;
     bondG[2] = c.pp2;
-    local_change_bond(current_vs_pid - 1, bondG, 0);
-    local_change_bond(current_vs_pid - 2, bondG, 0);
+    // Only add bond if vs was created on this node
+    if (local_particles[current_vs_pid - 1])
+      local_change_bond(current_vs_pid - 1, bondG, 0);
+    if (local_particles[current_vs_pid - 2])
+      local_change_bond(current_vs_pid - 2, bondG, 0);
     break;
   }
   }
@@ -516,15 +522,17 @@ void three_particle_binding_domain_decomposition(
   for (auto &c : gathered_queue) {
     // If we have both particles, at least as ghosts, Get the corresponding cell
     // indices
-    if ((local_particles[c.pp1] != NULL) && (local_particles[c.pp2] != NULL)) {
+    if ((local_particles[c.pp1]) && (local_particles[c.pp2])) {
       Particle &p1 = *local_particles[c.pp1];
       Particle &p2 = *local_particles[c.pp2];
-      auto cell1 = cell_structure.position_to_cell(p1.r.p.data());
-      auto cell2 = cell_structure.position_to_cell(p2.r.p.data());
+      auto cell1 = find_current_cell(p1);
+      auto cell2 = find_current_cell(p2);
 
-      three_particle_binding_do_search(cell1, p1, p2);
-      if (cell1 != cell2)
+      if (cell1)
+        three_particle_binding_do_search(cell1, p1, p2);
+      if (cell2 && cell1 != cell2)
         three_particle_binding_do_search(cell2, p1, p2);
+
     } // If local particles exist
   }   // Loop over total collisions
 }
@@ -540,6 +548,10 @@ void handle_collisions() {
     }
   }
 
+  // Note that the glue to surface mode adds bonds between the centers
+  // but does so later in the process. This is needed to gaurantee that
+  // a particle can only be glued once, even if queued twice in a single
+  // time step
   if (bind_centers()) {
     for (auto &c : local_collision_queue) {
       // put the bond to the non-ghost particle; at least one partner always is
@@ -576,21 +588,38 @@ void handle_collisions() {
     for (auto &c : gathered_queue) {
 
       // Get particle pointers
-      Particle *const p1 = local_particles[c.pp1];
-      Particle *const p2 = local_particles[c.pp2];
-
-      // If we have none of the two partic;es, only increase the counter for the
-      // next id to use
-      if (p1 == NULL and p2 == NULL) {
+      Particle *p1 = local_particles[c.pp1];
+      Particle *p2 = local_particles[c.pp2];
+
+      // Only nodes take part in particle creation and binding
+      // that see both particles
+
+      // If we cannot access both particles, both are ghosts,
+      // ore one is ghost and one is not accessible
+      // we only increase the counter for the ext id to use based on the
+      // number of particles created by other nodes
+      if (((!p1 or p1->l.ghost) and (!p2 or p2->l.ghost)) or !p1 or !p2) {
+        // Increase local counters
         if (collision_params.mode & COLLISION_MODE_VS) {
           added_particle(current_vs_pid);
           current_vs_pid++;
         }
         // For glue to surface, we have only one vs
         added_particle(current_vs_pid);
         current_vs_pid++;
-
-      } else { // We have at least one
+        if (collision_params.mode == COLLISION_MODE_GLUE_TO_SURF) {
+          if (p1)
+            if (p1->p.type == collision_params.part_type_to_be_glued) {
+              p1->p.type = collision_params.part_type_after_glueing;
+            }
+          if (p2)
+            if (p2->p.type == collision_params.part_type_to_be_glued) {
+              p2->p.type = collision_params.part_type_after_glueing;
+            }
+        } // mode glue to surface
+
+      } else { // We consider the pair because one particle
+               // is local to the node and the other is local or ghost
 
         // Calculate initial position for new vs, which is in the local node's
         // domain
@@ -611,27 +640,71 @@ void handle_collisions() {
           p1->p.rotation = ROTATION_X | ROTATION_Y | ROTATION_Z;
           p2->p.rotation = ROTATION_X | ROTATION_Y | ROTATION_Z;
 
-          // The vs placement is done by the node on which p1 is non-ghost
+          // Positions of the virtual sites
+          bind_at_point_of_collision_calc_vs_pos(p1, p2, pos1, pos2);
+
+          // place virtual sites on the node where the base particle is not a
+          // ghost
           if (!p1->l.ghost) {
-            bind_at_point_of_collision_calc_vs_pos(p1, p2, pos1, pos2);
             place_vs_and_relate_to_particle(current_vs_pid, pos1, c.pp1,
                                             initial_pos);
-            current_vs_pid++;
+            // Particle storage locations may have changed due to
+            // added particle
+            p1 = local_particles[c.pp1];
+            p2 = local_particles[c.pp2];
+          } else // update the books
+            added_particle(current_vs_pid);
+
+          // Increment counter
+          current_vs_pid++;
+
+          // Same for particle 2
+          if (!p2->l.ghost) {
             place_vs_and_relate_to_particle(current_vs_pid, pos2, c.pp2,
                                             initial_pos);
-            current_vs_pid++;
-            bind_at_poc_create_bond_between_vs(current_vs_pid, c);
-          } else { // Just update the books
+            // Particle storage locations may have changed due to
+            // added particle
+            p1 = local_particles[c.pp1];
+            p2 = local_particles[c.pp2];
+          } else // update the books
             added_particle(current_vs_pid);
-            current_vs_pid++;
 
-            added_particle(current_vs_pid);
-            current_vs_pid++;
-          }
-        }
+          // Increment counter
+          current_vs_pid++;
+
+          // Create bonds between the vs.
+
+          bind_at_poc_create_bond_between_vs(current_vs_pid, c);
+        } // mode VS
+
         if (collision_params.mode & COLLISION_MODE_GLUE_TO_SURF) {
+          // If particles are made inert by a type change on collision:
+          // We skip the pair if one of the particles has already reacted
+          // but we still increase the particle counters, as other nodes
+          // can not always know whether or not a vs is placed
+          if (collision_params.part_type_after_glueing !=
+              collision_params.part_type_to_be_glued) {
+            if ((p1->p.type == collision_params.part_type_after_glueing) ||
+                (p2->p.type == collision_params.part_type_after_glueing)) {
+
+              added_particle(current_vs_pid);
+              current_vs_pid++;
+              continue;
+            }
+          }
+
           Vector3d pos;
-          const int pid = glue_to_surface_calc_vs_pos(p1, p2, pos);
+          const Particle &attach_vs_to =
+              glue_to_surface_calc_vs_pos(*p1, *p2, pos);
+
+          // Add a bond between the centers of the colliding particles
+          // The bond is placed on the node that has p1
+          if (!p1->l.ghost) {
+            int bondG[2];
+            bondG[0] = collision_params.bond_centers;
+            bondG[1] = c.pp2;
+            local_change_bond(c.pp1, bondG, 0);
+          }
 
           // Change type of partilce being attached, to make it inert
           if (p1->p.type == collision_params.part_type_to_be_glued) {
@@ -642,27 +715,39 @@ void handle_collisions() {
           }
 
           // Vs placement happens on the node that has p1
-          if (!p1->l.ghost) {
-            place_vs_and_relate_to_particle(current_vs_pid, pos, pid,
-                                            initial_pos);
+          if (!attach_vs_to.l.ghost) {
+            place_vs_and_relate_to_particle(
+                current_vs_pid, pos, attach_vs_to.p.identity, initial_pos);
+            // Particle storage locations may have changed due to
+            // added particle
+            p1 = local_particles[c.pp1];
+            p2 = local_particles[c.pp2];
             current_vs_pid++;
           } else { // Just update the books
             added_particle(current_vs_pid);
             current_vs_pid++;
           }
           glue_to_surface_bind_part_to_vs(p1, p2, current_vs_pid, c);
         }
-      } // Loop over all collisions in the queue
+      } // we considered the pair
+    }   // Loop over all collisions in the queue
+#ifdef ADDITIONAL_CHECKS
+    if (!Utils::Mpi::all_compare(comm_cart, current_vs_pid)) {
+      throw std::runtime_error("Nodes disagree about current_vs_pid");
+    }
+    if (!Utils::Mpi::all_compare(comm_cart, max_seen_particle)) {
+      throw std::runtime_error("Nodes disagree about max_seen_particle");
+    }
+#endif
 
-      // If any node had a collision, all nodes need to do on_particle_change
-      // and resort
+    // If any node had a collision, all nodes need to do on_particle_change
+    // and resort
 
-      if (gathered_queue.size() > 0) {
-        on_particle_change();
-        announce_resort_particles();
-        cells_update_ghosts();
-      }
-    }  // total_collision>0
+    if (gathered_queue.size() > 0) {
+      on_particle_change();
+      announce_resort_particles();
+      cells_update_ghosts();
+    }
   }    // are we in one of the vs_based methods
 #endif // defined VIRTUAL_SITES_RELATIVE
 

src/core/domain_decomposition.cpp

@@ -31,7 +31,7 @@
 /** Returns pointer to the cell which corresponds to the position if
     the position is in the nodes spatial domain otherwise a nullptr
     pointer. */
-Cell *dd_save_position_to_cell(double pos[3]);
+Cell *dd_save_position_to_cell(const double pos[3]);
 
 /************************************************/
 /** \name Defines */
@@ -526,7 +526,7 @@ void dd_init_cell_interactions() {
 /** Returns pointer to the cell which corresponds to the position if
     the position is in the nodes spatial domain otherwise a nullptr
     pointer. */
-Cell *dd_save_position_to_cell(double pos[3]) {
+Cell *dd_save_position_to_cell(const double pos[3]) {
   int i, cpos[3];
 
   for (i = 0; i < 3; i++) {