Added the ghost task for black holes.

src/black_holes/Default/black_holes_iact.h

@@ -31,12 +31,10 @@
  * @param a Current scale factor.
  * @param H Current Hubble parameter.
  */
-__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_bh_density(const float r2, const float *dx,
-                                 const float hi, const float hj,
-                                 struct bpart *restrict bi,
-                                 const struct part *restrict pj, const float a,
-                                 const float H) {
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_bh_density(
+    const float r2, const float *dx, const float hi, const float hj,
+    struct bpart *restrict bi, const struct part *restrict pj, const float a,
+    const float H) {
 
   float wi, wi_dx;
 
@@ -76,11 +74,10 @@ runner_iact_nonsym_bh_density(const float r2, const float *dx,
  * @param H Current Hubble parameter.
  */
 __attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_bh_feedback(const float r2, const float *dx,
-                                  const float hi, const float hj,
-                                  struct bpart *restrict bi,
-                                  struct part *restrict pj, const float a,
-                                  const float H) {
+runner_iact_nonsym_bh_feedback(const float r2, const float *dx, const float hi,
+                               const float hj, struct bpart *restrict bi,
+                               struct part *restrict pj, const float a,
+                               const float H) {
 #ifdef DEBUG_INTERACTIONS_BH
   /* Update ngb counters */
   if (si->num_ngb_force < MAX_NUM_OF_NEIGHBOURS_BH)

src/runner.c

@@ -500,6 +500,314 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
   if (timer) TIMER_TOC(timer_do_stars_ghost);
 }
 
+/**
+ * @brief Intermediate task after the density to check that the smoothing
+ * lengths are correct.
+ *
+ * @param r The runner thread.
+ * @param c The cell.
+ * @param timer Are we timing this ?
+ */
+void runner_do_black_holes_ghost(struct runner *r, struct cell *c, int timer) {
+
+  struct bpart *restrict bparts = c->black_holes.parts;
+  const struct engine *e = r->e;
+  const struct cosmology *cosmo = e->cosmology;
+  const float black_holes_h_max = e->hydro_properties->h_max;
+  const float black_holes_h_min = e->hydro_properties->h_min;
+  const float eps = e->hydro_properties->h_tolerance;
+  const float black_holes_eta_dim =
+      pow_dimension(e->hydro_properties->eta_neighbours);
+  const int max_smoothing_iter = e->hydro_properties->max_smoothing_iterations;
+  int redo = 0, bcount = 0;
+
+  /* Running value of the maximal smoothing length */
+  double h_max = c->black_holes.h_max;
+
+  TIMER_TIC;
+
+  /* Anything to do here? */
+  if (c->black_holes.count == 0) return;
+  if (!cell_is_active_black_holes(c, e)) return;
+
+  /* Recurse? */
+  if (c->split) {
+    for (int k = 0; k < 8; k++) {
+      if (c->progeny[k] != NULL) {
+        runner_do_black_holes_ghost(r, c->progeny[k], 0);
+
+        /* Update h_max */
+        h_max = max(h_max, c->progeny[k]->black_holes.h_max);
+      }
+    }
+  } else {
+
+    /* Init the list of active particles that have to be updated. */
+    int *sid = NULL;
+    float *h_0 = NULL;
+    float *left = NULL;
+    float *right = NULL;
+    if ((sid = (int *)malloc(sizeof(int) * c->black_holes.count)) == NULL)
+      error("Can't allocate memory for sid.");
+    if ((h_0 = (float *)malloc(sizeof(float) * c->black_holes.count)) == NULL)
+      error("Can't allocate memory for h_0.");
+    if ((left = (float *)malloc(sizeof(float) * c->black_holes.count)) == NULL)
+      error("Can't allocate memory for left.");
+    if ((right = (float *)malloc(sizeof(float) * c->black_holes.count)) == NULL)
+      error("Can't allocate memory for right.");
+    for (int k = 0; k < c->black_holes.count; k++)
+      if (bpart_is_active(&bparts[k], e)) {
+        sid[bcount] = k;
+        h_0[bcount] = bparts[k].h;
+        left[bcount] = 0.f;
+        right[bcount] = black_holes_h_max;
+        ++bcount;
+      }
+
+    /* While there are particles that need to be updated... */
+    for (int num_reruns = 0; bcount > 0 && num_reruns < max_smoothing_iter;
+         num_reruns++) {
+
+      /* Reset the redo-count. */
+      redo = 0;
+
+      /* Loop over the remaining active parts in this cell. */
+      for (int i = 0; i < bcount; i++) {
+
+        /* Get a direct pointer on the part. */
+        struct bpart *bp = &bparts[sid[i]];
+
+#ifdef SWIFT_DEBUG_CHECKS
+        /* Is this part within the timestep? */
+        if (!bpart_is_active(bp, e))
+          error("Ghost applied to inactive particle");
+#endif
+
+        /* Get some useful values */
+        const float h_init = h_0[i];
+        const float h_old = bp->h;
+        const float h_old_dim = pow_dimension(h_old);
+        const float h_old_dim_minus_one = pow_dimension_minus_one(h_old);
+
+        float h_new;
+        int has_no_neighbours = 0;
+
+        if (bp->density.wcount == 0.f) { /* No neighbours case */
+
+          /* Flag that there were no neighbours */
+          has_no_neighbours = 1;
+
+          /* Double h and try again */
+          h_new = 2.f * h_old;
+
+        } else {
+
+          /* Finish the density calculation */
+          black_holes_end_density(bp, cosmo);
+
+          /* Compute one step of the Newton-Raphson scheme */
+          const float n_sum = bp->density.wcount * h_old_dim;
+          const float n_target = black_holes_eta_dim;
+          const float f = n_sum - n_target;
+          const float f_prime =
+              bp->density.wcount_dh * h_old_dim +
+              hydro_dimension * bp->density.wcount * h_old_dim_minus_one;
+
+          /* Improve the bisection bounds */
+          if (n_sum < n_target)
+            left[i] = max(left[i], h_old);
+          else if (n_sum > n_target)
+            right[i] = min(right[i], h_old);
+
+#ifdef SWIFT_DEBUG_CHECKS
+          /* Check the validity of the left and right bounds */
+          if (left[i] > right[i])
+            error("Invalid left (%e) and right (%e)", left[i], right[i]);
+#endif
+
+          /* Skip if h is already h_max and we don't have enough neighbours */
+          /* Same if we are below h_min */
+          if (((bp->h >= black_holes_h_max) && (f < 0.f)) ||
+              ((bp->h <= black_holes_h_min) && (f > 0.f))) {
+
+            black_holes_reset_feedback(bp);
+
+            /* Ok, we are done with this particle */
+            continue;
+          }
+
+          /* Normal case: Use Newton-Raphson to get a better value of h */
+
+          /* Avoid floating point exception from f_prime = 0 */
+          h_new = h_old - f / (f_prime + FLT_MIN);
+
+          /* Be verbose about the particles that struggle to converge */
+          if (num_reruns > max_smoothing_iter - 10) {
+
+            message(
+                "Smoothing length convergence problem: iter=%d p->id=%lld "
+                "h_init=%12.8e h_old=%12.8e h_new=%12.8e f=%f f_prime=%f "
+                "n_sum=%12.8e n_target=%12.8e left=%12.8e right=%12.8e",
+                num_reruns, bp->id, h_init, h_old, h_new, f, f_prime, n_sum,
+                n_target, left[i], right[i]);
+          }
+
+          /* Safety check: truncate to the range [ h_old/2 , 2h_old ]. */
+          h_new = min(h_new, 2.f * h_old);
+          h_new = max(h_new, 0.5f * h_old);
+
+          /* Verify that we are actually progrssing towards the answer */
+          h_new = max(h_new, left[i]);
+          h_new = min(h_new, right[i]);
+        }
+
+        /* Check whether the particle has an inappropriate smoothing length */
+        if (fabsf(h_new - h_old) > eps * h_old) {
+
+          /* Ok, correct then */
+
+          /* Case where we have been oscillating around the solution */
+          if ((h_new == left[i] && h_old == right[i]) ||
+              (h_old == left[i] && h_new == right[i])) {
+
+            /* Bissect the remaining interval */
+            bp->h = pow_inv_dimension(
+                0.5f * (pow_dimension(left[i]) + pow_dimension(right[i])));
+
+          } else {
+
+            /* Normal case */
+            bp->h = h_new;
+          }
+
+          /* If below the absolute maximum, try again */
+          if (bp->h < black_holes_h_max && bp->h > black_holes_h_min) {
+
+            /* Flag for another round of fun */
+            sid[redo] = sid[i];
+            h_0[redo] = h_0[i];
+            left[redo] = left[i];
+            right[redo] = right[i];
+            redo += 1;
+
+            /* Re-initialise everything */
+            black_holes_init_bpart(bp);
+
+            /* Off we go ! */
+            continue;
+
+          } else if (bp->h <= black_holes_h_min) {
+
+            /* Ok, this particle is a lost cause... */
+            bp->h = black_holes_h_min;
+
+          } else if (bp->h >= black_holes_h_max) {
+
+            /* Ok, this particle is a lost cause... */
+            bp->h = black_holes_h_max;
+
+            /* Do some damage control if no neighbours at all were found */
+            if (has_no_neighbours) {
+              black_holes_bpart_has_no_neighbours(bp, cosmo);
+            }
+
+          } else {
+            error(
+                "Fundamental problem with the smoothing length iteration "
+                "logic.");
+          }
+        }
+
+        /* We now have a particle whose smoothing length has converged */
+
+        /* Check if h_max has increased */
+        h_max = max(h_max, bp->h);
+
+        black_holes_reset_feedback(bp);
+      }
+
+      /* We now need to treat the particles whose smoothing length had not
+       * converged again */
+
+      /* Re-set the counter for the next loop (potentially). */
+      bcount = redo;
+      if (bcount > 0) {
+
+        /* Climb up the cell hierarchy. */
+        for (struct cell *finger = c; finger != NULL; finger = finger->parent) {
+
+          /* Run through this cell's density interactions. */
+          for (struct link *l = finger->black_holes.density; l != NULL;
+               l = l->next) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+            if (l->t->ti_run < r->e->ti_current)
+              error("Density task should have been run.");
+#endif
+
+            /* Self-interaction? */
+            if (l->t->type == task_type_self)
+              runner_doself_subset_branch_bh_density(r, finger, bparts, sid,
+                                                     bcount);
+
+            /* Otherwise, pair interaction? */
+            else if (l->t->type == task_type_pair) {
+
+              /* Left or right? */
+              if (l->t->ci == finger)
+                runner_dopair_subset_branch_bh_density(r, finger, bparts, sid,
+                                                       bcount, l->t->cj);
+              else
+                runner_dopair_subset_branch_bh_density(r, finger, bparts, sid,
+                                                       bcount, l->t->ci);
+            }
+
+            /* Otherwise, sub-self interaction? */
+            else if (l->t->type == task_type_sub_self)
+              runner_dosub_subset_bh_density(r, finger, bparts, sid, bcount,
+                                             NULL, 1);
+
+            /* Otherwise, sub-pair interaction? */
+            else if (l->t->type == task_type_sub_pair) {
+
+              /* Left or right? */
+              if (l->t->ci == finger)
+                runner_dosub_subset_bh_density(r, finger, bparts, sid, bcount,
+                                               l->t->cj, 1);
+              else
+                runner_dosub_subset_bh_density(r, finger, bparts, sid, bcount,
+                                               l->t->ci, 1);
+            }
+          }
+        }
+      }
+    }
+
+    if (bcount) {
+      error("Smoothing length failed to converge on %i particles.", bcount);
+    }
+
+    /* Be clean */
+    free(left);
+    free(right);
+    free(sid);
+    free(h_0);
+  }
+
+  /* Update h_max */
+  c->black_holes.h_max = h_max;
+
+  /* The ghost may not always be at the top level.
+   * Therefore we need to update h_max between the super- and top-levels */
+  if (c->black_holes.ghost) {
+    for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) {
+      atomic_max_d(&tmp->black_holes.h_max, h_max);
+    }
+  }
+
+  if (timer) TIMER_TOC(timer_do_black_holes_ghost);
+}
+
 /**
  * @brief Calculate gravity acceleration from external potential
  *
@@ -3656,6 +3964,9 @@ void *runner_main(void *data) {
         case task_type_stars_ghost:
           runner_do_stars_ghost(r, ci, 1);
           break;
+        case task_type_bh_ghost:
+          runner_do_black_holes_ghost(r, ci, 1);
+          break;
         case task_type_drift_part:
           runner_do_drift_part(r, ci, 1);
           break;

src/runner_doiact_black_holes.h

@@ -31,8 +31,7 @@
 #define _DO_SYM_PAIR1_BH(f) PASTE(runner_do_sym_pair_bh, f)
 #define DO_SYM_PAIR1_BH _DO_SYM_PAIR1_BH(FUNCTION)
 
-#define _DO_NONSYM_PAIR1_BH_NAIVE(f) \
-  PASTE(runner_do_nonsym_pair_bh_naive, f)
+#define _DO_NONSYM_PAIR1_BH_NAIVE(f) PASTE(runner_do_nonsym_pair_bh_naive, f)
 #define DO_NONSYM_PAIR1_BH_NAIVE _DO_NONSYM_PAIR1_BH_NAIVE(FUNCTION)
 
 #define _DOPAIR1_BH_NAIVE(f) PASTE(runner_dopair_bh_naive, f)
@@ -41,19 +40,16 @@
 #define _DOPAIR1_SUBSET_BH(f) PASTE(runner_dopair_subset_bh, f)
 #define DOPAIR1_SUBSET_BH _DOPAIR1_SUBSET_BH(FUNCTION)
 
-#define _DOPAIR1_SUBSET_BH_NAIVE(f) \
-  PASTE(runner_dopair_subset_bh_naive, f)
+#define _DOPAIR1_SUBSET_BH_NAIVE(f) PASTE(runner_dopair_subset_bh_naive, f)
 #define DOPAIR1_SUBSET_BH_NAIVE _DOPAIR1_SUBSET_BH_NAIVE(FUNCTION)
 
 #define _DOSELF1_SUBSET_BH(f) PASTE(runner_doself_subset_bh, f)
 #define DOSELF1_SUBSET_BH _DOSELF1_SUBSET_BH(FUNCTION)
 
-#define _DOSELF1_SUBSET_BRANCH_BH(f) \
-  PASTE(runner_doself_subset_branch_bh, f)
+#define _DOSELF1_SUBSET_BRANCH_BH(f) PASTE(runner_doself_subset_branch_bh, f)
 #define DOSELF1_SUBSET_BRANCH_BH _DOSELF1_SUBSET_BRANCH_BH(FUNCTION)
 
-#define _DOPAIR1_SUBSET_BRANCH_BH(f) \
-  PASTE(runner_dopair_subset_branch_bh, f)
+#define _DOPAIR1_SUBSET_BRANCH_BH(f) PASTE(runner_dopair_subset_branch_bh, f)
 #define DOPAIR1_SUBSET_BRANCH_BH _DOPAIR1_SUBSET_BRANCH_BH(FUNCTION)
 
 #define _DOSUB_SUBSET_BH(f) PASTE(runner_dosub_subset_bh, f)
@@ -102,8 +98,8 @@ void DOSELF1_BH(struct runner *r, struct cell *c, int timer) {
   TIMER_TIC;
 
   const struct engine *e = r->e;
-  //const int with_cosmology = e->policy & engine_policy_cosmology;
-  //const integertime_t ti_current = e->ti_current;
+  // const int with_cosmology = e->policy & engine_policy_cosmology;
+  // const integertime_t ti_current = e->ti_current;
   const struct cosmology *cosmo = e->cosmology;
 
   /* Anything to do here? */
@@ -118,7 +114,7 @@ void DOSELF1_BH(struct runner *r, struct cell *c, int timer) {
   const int count = c->hydro.count;
   struct bpart *restrict bparts = c->black_holes.parts;
   struct part *restrict parts = c->hydro.parts;
-  //struct xpart *restrict xparts = c->hydro.xparts;
+  // struct xpart *restrict xparts = c->hydro.xparts;
 
   /* Loop over the bparts in ci. */
   for (int sid = 0; sid < bcount; sid++) {
@@ -140,7 +136,7 @@ void DOSELF1_BH(struct runner *r, struct cell *c, int timer) {
 
       /* Get a pointer to the jth particle. */
       struct part *restrict pj = &parts[pjd];
-      //struct xpart *restrict xpj = &xparts[pjd];
+      // struct xpart *restrict xpj = &xparts[pjd];
       const float hj = pj->h;
 
       /* Early abort? */
@@ -176,7 +172,7 @@ void DOSELF1_BH(struct runner *r, struct cell *c, int timer) {
  * @param cj The second #cell
  */
 void DO_NONSYM_PAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
-                                 struct cell *restrict cj) {
+                              struct cell *restrict cj) {
 
 #ifdef SWIFT_DEBUG_CHECKS
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
@@ -187,8 +183,8 @@ void DO_NONSYM_PAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
 #endif
 
   const struct engine *e = r->e;
-  //const int with_cosmology = e->policy & engine_policy_cosmology;
-  //const integertime_t ti_current = e->ti_current;
+  // const int with_cosmology = e->policy & engine_policy_cosmology;
+  // const integertime_t ti_current = e->ti_current;
   const struct cosmology *cosmo = e->cosmology;
 
   /* Anything to do here? */
@@ -203,7 +199,7 @@ void DO_NONSYM_PAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
   const int count_j = cj->hydro.count;
   struct bpart *restrict bparts_i = ci->black_holes.parts;
   struct part *restrict parts_j = cj->hydro.parts;
-  //struct xpart *restrict xparts_j = cj->hydro.xparts;
+  // struct xpart *restrict xparts_j = cj->hydro.xparts;
 
   /* Get the relative distance between the pairs, wrapping. */
   double shift[3] = {0.0, 0.0, 0.0};
@@ -234,7 +230,7 @@ void DO_NONSYM_PAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
 
       /* Get a pointer to the jth particle. */
       struct part *restrict pj = &parts_j[pjd];
-      //struct xpart *restrict xpj = &xparts_j[pjd];
+      // struct xpart *restrict xpj = &xparts_j[pjd];
       const float hj = pj->h;
 
       /* Skip inhibited particles. */
@@ -261,7 +257,7 @@ void DO_NONSYM_PAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
 }
 
 void DOPAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
-                         struct cell *restrict cj, int timer) {
+                      struct cell *restrict cj, int timer) {
 
   TIMER_TIC;
 
@@ -289,23 +285,23 @@ void DOPAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
  *
  * @param r The #runner.
  * @param ci The first #cell.
- * @param sparts_i The #part to interact with @c cj.
+ * @param bparts_i The #part to interact with @c cj.
  * @param ind The list of indices of particles in @c ci to interact with.
- * @param scount The number of particles in @c ind.
+ * @param bcount The number of particles in @c ind.
  * @param cj The second #cell.
  * @param shift The shift vector to apply to the particles in ci.
  */
 void DOPAIR1_SUBSET_BH_NAIVE(struct runner *r, struct cell *restrict ci,
-                                struct bpart *restrict bparts_i,
-                                int *restrict ind, int bcount,
-                                struct cell *restrict cj, const double *shift) {
+                             struct bpart *restrict bparts_i, int *restrict ind,
+                             const int bcount, struct cell *restrict cj,
+                             const double *shift) {
 
 #ifdef SWIFT_DEBUG_CHECKS
   if (ci->nodeID != engine_rank) error("Should be run on a different node");
 #endif
 
   const struct engine *e = r->e;
-  //const integertime_t ti_current = e->ti_current;
+  // const integertime_t ti_current = e->ti_current;
   const struct cosmology *cosmo = e->cosmology;
 
   /* Cosmological terms */
@@ -314,7 +310,7 @@ void DOPAIR1_SUBSET_BH_NAIVE(struct runner *r, struct cell *restrict ci,
 
   const int count_j = cj->hydro.count;
   struct part *restrict parts_j = cj->hydro.parts;
-  //struct xpart *restrict xparts_j = cj->hydro.xparts;
+  // struct xpart *restrict xparts_j = cj->hydro.xparts;
 
   /* Early abort? */
   if (count_j == 0) return;
@@ -341,7 +337,7 @@ void DOPAIR1_SUBSET_BH_NAIVE(struct runner *r, struct cell *restrict ci,
 
       /* Get a pointer to the jth particle. */
       struct part *restrict pj = &parts_j[pjd];
-      //struct xpart *restrict xpj = &xparts_j[pjd];
+      // struct xpart *restrict xpj = &xparts_j[pjd];
 
       /* Skip inhibited particles */
       if (part_is_inhibited(pj, e)) continue;
@@ -375,20 +371,20 @@ void DOPAIR1_SUBSET_BH_NAIVE(struct runner *r, struct cell *restrict ci,
  *
  * @param r The #runner.
  * @param ci The first #cell.
- * @param sparts The #spart to interact.
+ * @param bparts The #spart to interact.
  * @param ind The list of indices of particles in @c ci to interact with.
- * @param scount The number of particles in @c ind.
+ * @param bcount The number of particles in @c ind.
  */
 void DOSELF1_SUBSET_BH(struct runner *r, struct cell *restrict ci,
-                          struct bpart *restrict bparts, int *restrict ind,
-                          int bcount) {
+                       struct bpart *restrict bparts, int *restrict ind,
+                       const int bcount) {
 
 #ifdef SWIFT_DEBUG_CHECKS
   if (ci->nodeID != engine_rank) error("Should be run on a different node");
 #endif
 
   const struct engine *e = r->e;
-  //const integertime_t ti_current = e->ti_current;
+  // const integertime_t ti_current = e->ti_current;
   const struct cosmology *cosmo = e->cosmology;
 
   /* Cosmological terms */
@@ -397,7 +393,7 @@ void DOSELF1_SUBSET_BH(struct runner *r, struct cell *restrict ci,
 
   const int count_i = ci->hydro.count;
   struct part *restrict parts_j = ci->hydro.parts;
-  //struct xpart *restrict xparts_j = ci->hydro.xparts;
+  // struct xpart *restrict xparts_j = ci->hydro.xparts;
 
   /* Early abort? */
   if (count_i == 0) return;
@@ -423,7 +419,7 @@ void DOSELF1_SUBSET_BH(struct runner *r, struct cell *restrict ci,
 
       /* Get a pointer to the jth particle. */
       struct part *restrict pj = &parts_j[pjd];
-      //struct xpart *restrict xpj = &xparts_j[pjd];
+      // struct xpart *restrict xpj = &xparts_j[pjd];
 
       /* Early abort? */
       if (part_is_inhibited(pj, e)) continue;
@@ -455,13 +451,13 @@ void DOSELF1_SUBSET_BH(struct runner *r, struct cell *restrict ci,
  *
  * @param r The #runner.
  * @param ci The first #cell.
- * @param sparts The #spart to interact.
+ * @param bparts The #spart to interact.
  * @param ind The list of indices of particles in @c ci to interact with.
- * @param scount The number of particles in @c ind.
+ * @param bcount The number of particles in @c ind.
  */
 void DOSELF1_SUBSET_BRANCH_BH(struct runner *r, struct cell *restrict ci,
-                                 struct bpart *restrict bparts,
-                                 int *restrict ind, int bcount) {
+                              struct bpart *restrict bparts, int *restrict ind,
+                              const int bcount) {
 
   DOSELF1_SUBSET_BH(r, ci, bparts, ind, bcount);
 }
@@ -473,15 +469,15 @@ void DOSELF1_SUBSET_BRANCH_BH(struct runner *r, struct cell *restrict ci,
  *
  * @param r The #runner.
  * @param ci The first #cell.
- * @param sparts_i The #spart to interact with @c cj.
+ * @param bparts_i The #spart to interact with @c cj.
  * @param ind The list of indices of particles in @c ci to interact with.
- * @param scount The number of particles in @c ind.
+ * @param bcount The number of particles in @c ind.
  * @param cj The second #cell.
  */
 void DOPAIR1_SUBSET_BRANCH_BH(struct runner *r, struct cell *restrict ci,
-                                 struct bpart *restrict bparts_i,
-                                 int *restrict ind, int bcount,
-                                 struct cell *restrict cj) {
+                              struct bpart *restrict bparts_i,
+                              int *restrict ind, int const bcount,
+                              struct cell *restrict cj) {
 
   const struct engine *e = r->e;
 
@@ -501,7 +497,8 @@ void DOPAIR1_SUBSET_BRANCH_BH(struct runner *r, struct cell *restrict ci,
 }
 
 void DOSUB_SUBSET_BH(struct runner *r, struct cell *ci, struct bpart *bparts,
-                        int *ind, int bcount, struct cell *cj, int gettimer) {
+                     int *ind, const int bcount, struct cell *cj,
+                     int gettimer) {
 
   const struct engine *e = r->e;
   struct space *s = e->s;
@@ -518,7 +515,8 @@ void DOSUB_SUBSET_BH(struct runner *r, struct cell *ci, struct bpart *bparts,
       if (ci->progeny[k] != NULL) {
         if (&bparts[ind[0]] >= &ci->progeny[k]->black_holes.parts[0] &&
             &bparts[ind[0]] <
-                &ci->progeny[k]->black_holes.parts[ci->progeny[k]->black_holes.count]) {
+                &ci->progeny[k]
+                     ->black_holes.parts[ci->progeny[k]->black_holes.count]) {
           sub = ci->progeny[k];
           break;
         }
@@ -562,10 +560,10 @@ void DOSUB_SUBSET_BH(struct runner *r, struct cell *ci, struct bpart *bparts,
         const int pjd = csp->pairs[k].pjd;
         if (ci->progeny[pid] == sub && cj->progeny[pjd] != NULL)
           DOSUB_SUBSET_BH(r, ci->progeny[pid], bparts, ind, bcount,
-                             cj->progeny[pjd], 0);
+                          cj->progeny[pjd], 0);
         if (ci->progeny[pid] != NULL && cj->progeny[pjd] == sub)
           DOSUB_SUBSET_BH(r, cj->progeny[pjd], bparts, ind, bcount,
-                             ci->progeny[pid], 0);
+                          ci->progeny[pid], 0);
       }
     }