diff --git a/src/black_holes/Default/black_holes_iact.h b/src/black_holes/Default/black_holes_iact.h
index 16d5599d82c6df6227be0489a487378e8d7b711f..382906487f12ced85221f7634a92ce1fa293a620 100644
--- a/src/black_holes/Default/black_holes_iact.h
+++ b/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)
diff --git a/src/runner.c b/src/runner.c
index 77f8679eb394a6f1fed90b2228fe59dcd3bad3c1..1542a420f917b3cb659ca7088f16bc15036e33c2 100644
--- a/src/runner.c
+++ b/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;
diff --git a/src/runner_doiact_black_holes.h b/src/runner_doiact_black_holes.h
index f4c939399bc28969201e30c368a4566f14a257fa..47604f2a1ea21d7e89f72a0ef4d441c45254b96c 100644
--- a/src/runner_doiact_black_holes.h
+++ b/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);
}
}
@@ -666,7 +664,7 @@ void DOPAIR1_BRANCH_BH(struct runner *r, struct cell *ci, struct cell *cj) {
* redundant computations to find the sid on-the-fly.
*/
void DOSUB_PAIR1_BH(struct runner *r, struct cell *ci, struct cell *cj,
- int gettimer) {
+ int gettimer) {
TIMER_TIC;
@@ -721,7 +719,7 @@ void DOSUB_PAIR1_BH(struct runner *r, struct cell *ci, struct cell *cj,
if (!cell_are_part_drifted(cj, e))
error("Interacting undrifted cells (parts).");
}
-
+
if (do_cj) {
/* Make sure both cells are drifted to the current timestep. */
diff --git a/src/task.c b/src/task.c
index cf1e2bab688dba29443267b33647c444b0cfd11e..c0a168277a0fcb532d07dade5eb3126af850d972 100644
--- a/src/task.c
+++ b/src/task.c
@@ -88,7 +88,8 @@ const char *taskID_names[task_type_count] = {"none",
"stars_ghost_in",
"stars_ghost",
"stars_ghost_out",
- "stars_sort"};
+ "stars_sort",
+ "bh_ghost"};
/* Sub-task type names. */
const char *subtaskID_names[task_subtype_count] = {
diff --git a/src/task.h b/src/task.h
index 3252ead471e1145e8fff70eb03a686dad16d50c3..3e1e3778240866389714620f02fe5e3742d70ef6 100644
--- a/src/task.h
+++ b/src/task.h
@@ -84,6 +84,7 @@ enum task_types {
task_type_stars_ghost,
task_type_stars_ghost_out, /* Implicit */
task_type_stars_sort,
+ task_type_bh_ghost,
task_type_count
} __attribute__((packed));