Use the max h of the particles in the subset rather than in the cell to decide...

Use the max h of the particles in the subset rather than in the cell to decide on  the recursion strategy in the hydro ghost

src/cell.h

@@ -713,6 +713,16 @@ cell_can_recurse_in_pair_hydro_task(const struct cell *c) {
                        c->hydro.dx_max_part_old) < 0.5f * c->dmin);
 }
 
+__attribute__((always_inline)) INLINE static int
+cell_can_recurse_in_subset_pair_hydro_task(const struct cell *c,
+                                           const float h_max_subset) {
+
+  /* If so, is the cut-off radius plus the max distance the parts have moved */
+  /* smaller than the sub-cell sizes ? */
+  return (kernel_gamma * h_max_subset + c->hydro.dx_max_part_old) <
+         0.5f * c->dmin;
+}
+
 /**
  * @brief Can a sub-pair hydro task recurse to a lower level based
  * on the status of the particles in the cell.
@@ -756,6 +766,20 @@ cell_can_recurse_in_self_hydro_task(const struct cell *c) {
   return c->split && (kernel_gamma * c->hydro.h_max_old < 0.5f * c->dmin);
 }
 
+/**
+ * @brief Can a hydro task recurse to a lower level based
+ * on the status of the particles in the cell.
+ *
+ * @param c The #cell.
+ */
+__attribute__((always_inline)) INLINE static int
+cell_can_recurse_in_self_subset_hydro_task(const struct cell *c,
+                                           const float h_max_subset) {
+
+  /* Is the cell not smaller than the smoothing length? */
+  return kernel_gamma * h_max_subset < 0.5f * c->dmin;
+}
+
 /**
  * @brief Can a sub-self hydro task recurse to a lower level based
  * on the status of the particles in the cell.

src/runner_doiact_functions_hydro.h

@@ -2838,7 +2838,7 @@ struct cell *FIND_SUB(const struct cell *const c,
 
 void DOSUB_PAIR_SUBSET(struct runner *r, struct cell *ci, struct part *parts,
                        const int *ind, const int count, struct cell *cj,
-                       const int gettimer) {
+                       const float h_max_subset, const int gettimer) {
 
   const struct engine *e = r->e;
   struct space *s = e->s;
@@ -2850,8 +2850,8 @@ void DOSUB_PAIR_SUBSET(struct runner *r, struct cell *ci, struct part *parts,
   if (!cell_is_active_hydro(ci, e)) return;
 
   /* Recurse? */
-  if (cell_can_recurse_in_pair_hydro_task(ci) &&
-      cell_can_recurse_in_pair_hydro_task(cj)) {
+  if (ci->split && cj->split &&
+      cell_can_recurse_in_subset_pair_hydro_task(ci, h_max_subset)) {
 
     /* Find in which sub-cell of ci the particles are */
     struct cell *const sub = FIND_SUB(ci, parts, ind);
@@ -2866,29 +2866,33 @@ void DOSUB_PAIR_SUBSET(struct runner *r, struct cell *ci, struct part *parts,
       const int pjd = csp->pairs[k].pjd;
       if (ci->progeny[pid] == sub && cj->progeny[pjd] != NULL)
         DOSUB_PAIR_SUBSET(r, ci->progeny[pid], parts, ind, count,
-                          cj->progeny[pjd],
+                          cj->progeny[pjd], h_max_subset,
                           /*gettimer=*/0);
       if (ci->progeny[pid] != NULL && cj->progeny[pjd] == sub)
         DOSUB_PAIR_SUBSET(r, cj->progeny[pjd], parts, ind, count,
-                          ci->progeny[pid],
+                          ci->progeny[pid], h_max_subset,
                           /*gettimer=*/0);
     }
 
-  }
-  /* Otherwise, compute the pair directly. */
-  else if (cell_is_active_hydro(ci, e)) {
+  } else if (cell_is_active_hydro(ci, e)) {
+
+    /* Otherwise, compute the pair directly. */
 
     /* Do any of the cells need to be drifted first? */
     if (!cell_are_part_drifted(cj, e)) error("Cell should be drifted!");
 
     DOPAIR_SUBSET_BRANCH(r, ci, parts, ind, count, cj);
+
+  } else {
+    error("Recursion logic");
   }
 
   if (gettimer) TIMER_TOC(timer_dosub_subset);
 }
 
 void DOSUB_SELF_SUBSET(struct runner *r, struct cell *ci, struct part *parts,
-                       const int *ind, const int count, const int gettimer) {
+                       const int *ind, const int count,
+                       const float h_max_subset, const int gettimer) {
 
   const struct engine *e = r->e;
 
@@ -2897,17 +2901,18 @@ void DOSUB_SELF_SUBSET(struct runner *r, struct cell *ci, struct part *parts,
   if (!cell_is_active_hydro(ci, e)) return;
 
   /* Recurse? */
-  if (ci->split && cell_can_recurse_in_self_hydro_task(ci)) {
+  if (ci->split &&
+      cell_can_recurse_in_self_subset_hydro_task(ci, h_max_subset)) {
 
     /* Find in which sub-cell of ci the particles are */
     struct cell *const sub = FIND_SUB(ci, parts, ind);
 
     /* Loop over all progeny. */
-    DOSUB_SELF_SUBSET(r, sub, parts, ind, count, /*gettimer=*/0);
+    DOSUB_SELF_SUBSET(r, sub, parts, ind, count, h_max_subset, /*gettimer=*/0);
     for (int j = 0; j < 8; j++)
       if (ci->progeny[j] != sub && ci->progeny[j] != NULL)
         DOSUB_PAIR_SUBSET(r, sub, parts, ind, count, ci->progeny[j],
-                          /*gettimer=*/0);
+                          h_max_subset, /*gettimer=*/0);
   }
 
   /* Otherwise, compute self-interaction. */

src/runner_doiact_hydro.h

@@ -161,7 +161,8 @@ void DOPAIR_SUBSET_BRANCH(struct runner *r, const struct cell *restrict ci,
 
 void DOSUB_PAIR_SUBSET(struct runner *r, struct cell *ci, struct part *parts,
                        const int *ind, const int count, struct cell *cj,
-                       const int gettimer);
+                       const float h_max_subset, const int gettimer);
 
 void DOSUB_SELF_SUBSET(struct runner *r, struct cell *ci, struct part *parts,
-                       const int *ind, const int count, const int gettimer);
+                       const int *ind, const int count,
+                       const float h_max_subset, const int gettimer);

src/runner_ghost.c

@@ -1457,7 +1457,8 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
 
             /* Otherwise, sub-self interaction? */
             else if (l->t->type == task_type_sub_self)
-              runner_dosub_self_subset_density(r, finger, parts, pid, count, 1);
+              runner_dosub_self_subset_density(r, finger, parts, pid, count,
+                                               h_max_subset, 1);
 
             /* Otherwise, sub-pair interaction? */
             else if (l->t->type == task_type_sub_pair) {
@@ -1465,10 +1466,10 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
               /* Left or right? */
               if (l->t->ci == finger)
                 runner_dosub_pair_subset_density(r, finger, parts, pid, count,
-                                                 l->t->cj, 1);
+                                                 l->t->cj, h_max_subset, 1);
               else
                 runner_dosub_pair_subset_density(r, finger, parts, pid, count,
-                                                 l->t->ci, 1);
+                                                 l->t->ci, h_max_subset, 1);
             }
           }
         }