Merge branch 'DOPAIR2_improvements' into 'master'

Reverse the loop to find the exit condition in DOPAIR2() to go from the centre to the outside

Follows @jwillis' suggestion in #362.

Also implements an early exit in the outer loop. If no particle can be in range any more (based on hi_max and hj_max), no need to keep checking the next particles.

See merge request !422

src/hydro/DebugInteractions/hydro.h

@@ -24,6 +24,7 @@
  * @brief Empty SPH implementation used solely to test the SELF/PAIR routines.
  */
 
+#include "equation_of_state.h"
 #include "hydro_properties.h"
 #include "hydro_space.h"
 #include "part.h"

src/runner_doiact.h

@@ -940,9 +940,17 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
   const int count_j = cj->count;
   struct part *restrict parts_i = ci->parts;
   struct part *restrict parts_j = cj->parts;
+
+  /* Maximal displacement since last rebuild */
+  const double dx_max = (ci->dx_max_sort + cj->dx_max_sort);
+
+  /* Position on the axis of the particles closest to the interface */
   const double di_max = sort_i[count_i - 1].d;
   const double dj_min = sort_j[0].d;
-  const double dx_max = (ci->dx_max_sort + cj->dx_max_sort);
+
+  /* Upper bound on maximal distance in the other cell */
+  const double max_i = max(hi_max, hj_max) * kernel_gamma + dx_max - rshift;
+  const double max_j = max(hi_max, hj_max) * kernel_gamma + dx_max;
 
   if (cell_is_active(ci, e)) {
 
@@ -950,6 +958,9 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
      * We start from the centre and move outwards. */
     for (int pid = count_i - 1; pid >= 0; pid--) {
 
+      /* Did we go beyond the upper bound ? */
+      if (sort_i[pid].d + max_i < dj_min) break;
+
       /* Get a hold of the ith part in ci. */
       struct part *pi = &parts_i[sort_i[pid].i];
       const float hi = pi->h;
@@ -963,11 +974,8 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
       if (di < dj_min) continue;
 
       /* Where do we have to stop when looping over cell j? */
-      int last_pj = count_j - 1;
-      while (sort_j[last_pj].d > di && last_pj > 0) last_pj--;
-
-      last_pj += 1;
-      if (last_pj >= count_j) last_pj = count_j - 1;
+      int last_pj = 0;
+      while (last_pj < count_j - 1 && sort_j[last_pj].d < di) last_pj++;
 
       /* Get some additional information about pi */
       const float hig2 = hi * hi * kernel_gamma2;
@@ -990,6 +998,15 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
         float dx[3] = {pix - pjx, piy - pjy, piz - pjz};
         const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
 
+#ifdef SWIFT_DEBUG_CHECKS
+        /* Check that particles have been drifted to the current time */
+        if (pi->ti_drift != e->ti_current)
+          error("Particle pi not drifted to current time");
+        if (pj->ti_drift != e->ti_current)
+          error("Particle pj not drifted to current time");
+#endif
+
+        /* Are these particles actually neighbours? */
         if (r2 < hig2 || r2 < hjg2) {
 
           IACT_NONSYM(r2, dx, hi, hj, pi, pj);
@@ -1004,6 +1021,9 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
      * We start from the centre and move outwards. */
     for (int pjd = 0; pjd < count_j; pjd++) {
 
+      /* Did we go beyond the upper bound ? */
+      if (sort_j[pjd].d - max_j > di_max - rshift) break;
+
       /* Get a hold of the jth part in cj. */
       struct part *pj = &parts_j[sort_j[pjd].i];
       const float hj = pj->h;
@@ -1016,12 +1036,8 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
       if (dj > di_max - rshift) continue;
 
       /* Where do we have to stop when looping over cell i? */
-      int first_pi = 0;
-      while (sort_i[first_pi].d - rshift < dj && first_pi < count_i - 1)
-        first_pi++;
-
-      first_pi -= 1;
-      if (first_pi < 0) first_pi = 0;
+      int first_pi = count_i;
+      while (first_pi > 0 && sort_i[first_pi].d - rshift > dj) first_pi--;
 
       /* Get some additional information about pj */
       const float hjg2 = hj * hj * kernel_gamma2;
@@ -1044,6 +1060,15 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
         float dx[3] = {pjx - pix, pjy - piy, pjz - piz};
         const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
 
+#ifdef SWIFT_DEBUG_CHECKS
+        /* Check that particles have been drifted to the current time */
+        if (pi->ti_drift != e->ti_current)
+          error("Particle pi not drifted to current time");
+        if (pj->ti_drift != e->ti_current)
+          error("Particle pj not drifted to current time");
+#endif
+
+        /* Are these particles actually neighbours? */
         if (r2 < hjg2 || r2 < hig2) {
 
           IACT_NONSYM(r2, dx, hj, hi, pj, pi);