Merge branch 'no_spart_memmove' into 'master'

Replaced the memmove in cell_add_spart() with a faster alternative

See merge request swift/swiftsim!2321

Author: MatthieuSchaller

src/cell_convert_part.c

@@ -29,44 +29,91 @@
 #include "active.h"
 #include "engine.h"
 #include "hydro.h"
+#include "memswap.h"
 #include "sink_properties.h"
 
 /**
  * @brief Recursively update the pointer and counter for #spart after the
  * addition of a new particle.
  *
+ * Recall that the particles are in depth-first ordering. This means all
+ * the leaves' content is next to each other in memory.
+ * To insert a empty spot in the desired leaf (the one at the bottom of
+ * the progeny list), we can:
+ * (1) take the first particle in a leaf and move it to the
+ * start of the next leaf.
+ * (2a) If the leaf is the one where we want the gap, increase the particle
+ * counter.
+ * (2b) If not, shift the particle array pointer by one spot to account
+ * for the gap which was created on the main progeny branch.
+ * (3) Since we know the delta a spart is moved, we can update its #gpart
+ * connection at the same time.
+ * (4) Since all the cells beyond the insertion point have had their
+ * particle order altered, we need to void the sort arrays.
+ * (5) In cells that are not leaves, don't touch the particles, but move the
+ * counter or pointer.
+ *
  * @param c The cell we are working on.
  * @param progeny_list The list of the progeny index at each level for the
  * leaf-cell where the particle was added.
- * @param main_branch Are we in a cell directly above the leaf where the new
- * particle was added?
+ * @param main_branch Are we in a cell directly above the leaf where we
+ * want to create a gap?
+ * @param spart_to_insert The #spart to insert at the start of the next
+ * leaf encountered in the hierarchy.
  */
 void cell_recursively_shift_sparts(struct cell *c,
                                    const int progeny_list[space_cell_maxdepth],
-                                   const int main_branch) {
+                                   const int main_branch,
+                                   struct spart *const spart_to_insert) {
+
+  /* Recurse to the lower levels */
   if (c->split) {
-    /* No need to recurse in progenies located before the insestion point */
+
+    /* No need to recurse in progenies located before the insertion point */
     const int first_progeny = main_branch ? progeny_list[(int)c->depth] : 0;
 
     for (int k = first_progeny; k < 8; ++k) {
-      if (c->progeny[k] != NULL)
+      if (c->progeny[k] != NULL) {
+
+        /* Is the progeny on the main branch? */
+        const int progeny_on_main_branch = main_branch && (k == first_progeny);
+
         cell_recursively_shift_sparts(c->progeny[k], progeny_list,
-                                      main_branch && (k == first_progeny));
+                                      progeny_on_main_branch, spart_to_insert);
+      }
     }
   }
 
-  /* When directly above the leaf with the new particle: increase the particle
-   * count */
-  /* When after the leaf with the new particle: shift by one position */
+  /* Indicate that the cell is not sorted and cancel the pointer sorting
+   * arrays. */
+  c->stars.sorted = 0;
+  cell_free_stars_sorts(c);
+
+  /* If we are in a leaf, we can operate on the particle arrays */
+  if (!c->split) {
+
+    /* If we are in a leaf, then we want to:
+     * - Swap the first particle in the array with the one in our copy buffer
+     * - Update the gpart link of the new particle in the buffer since we know
+     * where it is going (it is moving by exactly the cell's content).
+     */
+    if (c->stars.count > 0) {
+
+      memswap(c->stars.parts, spart_to_insert, sizeof(struct spart));
+      spart_to_insert->gpart->id_or_neg_offset -= c->stars.count;
+
+    } else {
+      /* Nothing to do */
+    }
+  }
+
+  /* Now we can increase the cell counter or shift the pointer*/
   if (main_branch) {
-    c->stars.count++;
 
-    /* Indicate that the cell is not sorted and cancel the pointer sorting
-     * arrays. */
-    c->stars.sorted = 0;
-    cell_free_stars_sorts(c);
+    c->stars.count++;
 
   } else {
+
     c->stars.parts++;
   }
 }
@@ -154,6 +201,7 @@ void cell_recursively_shift_gparts(struct cell *c,
  * time bin.
  */
 struct spart *cell_add_spart(struct engine *e, struct cell *const c) {
+
   /* Perform some basic consitency checks */
   if (c->nodeID != engine_rank) error("Adding spart on a foreign node");
   if (c->stars.ti_old_part != e->ti_current) error("Undrifted cell!");
@@ -205,39 +253,36 @@ struct spart *cell_add_spart(struct engine *e, struct cell *const c) {
     return NULL;
   }
 
+#ifdef SWIFT_DEBUG_CHECKS
   /* Number of particles to shift in order to get a free space. */
   const size_t n_copy = &top->stars.parts[top->stars.count] - c->stars.parts;
-
-#ifdef SWIFT_DEBUG_CHECKS
   if (c->stars.parts + n_copy > top->stars.parts + top->stars.count)
     error("Copying beyond the allowed range");
 #endif
 
-  if (n_copy > 0) {
-    // MATTHIEU: This can be improved. We don't need to copy everything, just
-    // need to swap a few particles.
-    memmove(&c->stars.parts[1], &c->stars.parts[0],
-            n_copy * sizeof(struct spart));
-
-    /* Update the spart->gpart links (shift by 1) */
-    for (size_t i = 0; i < n_copy; ++i) {
-
+  /* Recursively shift all the stars to get a free spot at the start of the
+   * current cell */
+  struct spart temp;
 #ifdef SWIFT_DEBUG_CHECKS
-      if (c->stars.parts[i + 1].gpart == NULL) {
-        error("Incorrectly linked spart!");
-      }
+  bzero(&temp, sizeof(struct spart));
+  temp.id = 666;
 #endif
-      c->stars.parts[i + 1].gpart->id_or_neg_offset--;
-    }
-  }
+  cell_recursively_shift_sparts(top, progeny, /* main_branch=*/1, &temp);
 
-  /* Recursively shift all the stars to get a free spot at the start of the
-   * current cell*/
-  cell_recursively_shift_sparts(top, progeny, /* main_branch=*/1);
+  /* Finish all the operations by putting the last particle
+   * we have in hands in the first spot beyond the top-level cell's
+   * current range.
+   * Note that top->stars.count has been updated in the recursion */
+  struct spart *final_spart = top->stars.parts + (top->stars.count - 1);
+#ifdef SWIFT_DEBUG_CHECKS
+  /* Verify that we are indeed overwriting an 'extra' spart */
+  if (final_spart->id != -42)
+    error("Incorrect shift beyond the original top-level range");
+#endif
+  memcpy(final_spart, &temp, sizeof(struct spart));
 
   /* Make sure the gravity will be recomputed for this particle in the next
-   * step
-   */
+   * step */
   struct cell *top2 = c;
   while (top2->parent != NULL) {
     top2->stars.ti_old_part = e->ti_current;
@@ -251,6 +296,9 @@ struct spart *cell_add_spart(struct engine *e, struct cell *const c) {
 
   /* We now have an empty spart as the first particle in that cell */
   struct spart *sp = &c->stars.parts[0];
+#ifdef SWIFT_DEBUG_CHECKS
+  if (sp->id != 666) error("Did not get the right particle!");
+#endif
   bzero(sp, sizeof(struct spart));
 
   /* Give it a decent position */