diff --git a/src/black_holes/EAGLE/black_holes_iact.h b/src/black_holes/EAGLE/black_holes_iact.h
index 3fbcb327fbf60038716335c7ff9fe166ad4d6b5a..0bfe2a51e8327dcb33e4d0f390cad530633094fd 100644
--- a/src/black_holes/EAGLE/black_holes_iact.h
+++ b/src/black_holes/EAGLE/black_holes_iact.h
@@ -230,6 +230,8 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
* smoothing length */
if (v2_pec < G * M / (kernel_gamma * h)) {
+ // MATTHIEU: Also add distance check: factor * Plummer softening
+
/* This particle is swallowed by the BH with the largest ID of all the
* candidates wanting to swallow it */
if (bj->merger_data.swallow_id < bi->id) {
diff --git a/src/cell.c b/src/cell.c
index 795467e3d9679f785bbd3b95417efbb69b9e7b0b..763f5117f17eec3927446e6a2dddab4b9b4ad4ec 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1988,7 +1988,8 @@ void cell_clean_links(struct cell *c, void *data) {
c->stars.feedback = NULL;
c->black_holes.density = NULL;
c->black_holes.swallow = NULL;
- c->black_holes.do_swallow = NULL;
+ c->black_holes.do_gas_swallow = NULL;
+ c->black_holes.do_bh_swallow = NULL;
c->black_holes.feedback = NULL;
}
@@ -4001,7 +4002,7 @@ int cell_unskip_black_holes_tasks(struct cell *c, struct scheduler *s) {
}
/* Un-skip the swallow tasks involved with this cell. */
- for (struct link *l = c->black_holes.do_swallow; l != NULL; l = l->next) {
+ for (struct link *l = c->black_holes.do_gas_swallow; l != NULL; l = l->next) {
struct task *t = l->t;
struct cell *ci = t->ci;
struct cell *cj = t->cj;
@@ -4049,23 +4050,14 @@ int cell_unskip_black_holes_tasks(struct cell *c, struct scheduler *s) {
/* Unskip all the other task types. */
if (c->nodeID == nodeID && cell_is_active_black_holes(c, e)) {
- /* for (struct link *l = c->black_holes.swallow; l != NULL; l = l->next)
- */
- /* scheduler_activate(s, l->t); */
- /* for (struct link *l = c->black_holes.do_swallow; l != NULL; l =
- * l->next)
- */
- /* scheduler_activate(s, l->t); */
- /* for (struct link *l = c->black_holes.feedback; l != NULL; l = l->next)
- */
- /* scheduler_activate(s, l->t); */
-
if (c->black_holes.density_ghost != NULL)
scheduler_activate(s, c->black_holes.density_ghost);
if (c->black_holes.swallow_ghost[0] != NULL)
scheduler_activate(s, c->black_holes.swallow_ghost[0]);
if (c->black_holes.swallow_ghost[1] != NULL)
scheduler_activate(s, c->black_holes.swallow_ghost[1]);
+ if (c->black_holes.swallow_ghost[2] != NULL)
+ scheduler_activate(s, c->black_holes.swallow_ghost[2]);
if (c->black_holes.black_holes_in != NULL)
scheduler_activate(s, c->black_holes.black_holes_in);
if (c->black_holes.black_holes_out != NULL)
diff --git a/src/cell.h b/src/cell.h
index 6fcc0afb79dc3b76c70acdcbde87c6d2264aba27..df5ef068501f54b7448582e13a5990ad115245c2 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -651,7 +651,7 @@ struct cell {
struct task *density_ghost;
/*! The star ghost task itself */
- struct task *swallow_ghost[2];
+ struct task *swallow_ghost[3];
/*! Linked list of the tasks computing this cell's BH density. */
struct link *density;
@@ -661,7 +661,10 @@ struct cell {
struct link *swallow;
/*! Linked list of the tasks processing the particles to swallow */
- struct link *do_swallow;
+ struct link *do_gas_swallow;
+
+ /*! Linked list of the tasks processing the particles to swallow */
+ struct link *do_bh_swallow;
/*! Linked list of the tasks computing this cell's BH feedback. */
struct link *feedback;
diff --git a/src/engine.c b/src/engine.c
index a42ead2d878f102dadaa925510f41296e48dde69..eeeca1992b7fce06b55952d7b2602210e911d67e 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -3327,11 +3327,13 @@ void engine_skip_force_and_kick(struct engine *e) {
t->type == task_type_extra_ghost ||
t->type == task_type_bh_swallow_ghost1 ||
t->type == task_type_bh_swallow_ghost2 ||
+ t->type == task_type_bh_swallow_ghost3 ||
t->subtype == task_subtype_gradient ||
t->subtype == task_subtype_stars_feedback ||
t->subtype == task_subtype_bh_feedback ||
t->subtype == task_subtype_bh_swallow ||
- t->subtype == task_subtype_do_swallow ||
+ t->subtype == task_subtype_do_gas_swallow ||
+ t->subtype == task_subtype_do_bh_swallow ||
t->subtype == task_subtype_bpart_rho ||
t->subtype == task_subtype_part_swallow ||
t->subtype == task_subtype_bpart_swallow ||
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index af119c9370d83eace384690b9b58aaa95e49867e..128dcb7e6cb6257086bf1b630a157c311ae73677 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -672,11 +672,10 @@ void engine_addtasks_recv_black_holes(struct engine *e, struct cell *c,
for (struct link *l = c->black_holes.swallow; l != NULL; l = l->next) {
scheduler_addunlock(s, t_rho, l->t);
- // scheduler_addunlock(s, l->t, t_swallow);
scheduler_addunlock(s, l->t, t_gas_swallow);
}
- for (struct link *l = c->black_holes.do_swallow; l != NULL; l = l->next) {
- // scheduler_addunlock(s, t_swallow, l->t);
+ for (struct link *l = c->black_holes.do_gas_swallow; l != NULL;
+ l = l->next) {
scheduler_addunlock(s, t_gas_swallow, l->t);
scheduler_addunlock(s, l->t, t_feedback);
}
@@ -1093,8 +1092,12 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c) {
c->black_holes.density_ghost = scheduler_addtask(
s, task_type_bh_density_ghost, task_subtype_none, 0, 0, c, NULL);
- c->black_holes.swallow_ghost[1] = scheduler_addtask(
- s, task_type_bh_swallow_ghost2, task_subtype_none, 0, 0, c, NULL);
+ c->black_holes.swallow_ghost[1] =
+ scheduler_addtask(s, task_type_bh_swallow_ghost2, task_subtype_none,
+ 0, /* implicit =*/1, c, NULL);
+
+ c->black_holes.swallow_ghost[2] = scheduler_addtask(
+ s, task_type_bh_swallow_ghost3, task_subtype_none, 0, 0, c, NULL);
scheduler_addunlock(s, c->super->kick2, c->black_holes.black_holes_in);
scheduler_addunlock(s, c->black_holes.black_holes_out,
@@ -1718,7 +1721,8 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
struct task *t_star_feedback = NULL;
struct task *t_bh_density = NULL;
struct task *t_bh_swallow = NULL;
- struct task *t_do_swallow = NULL;
+ struct task *t_do_gas_swallow = NULL;
+ struct task *t_do_bh_swallow = NULL;
struct task *t_bh_feedback = NULL;
for (int ind = 0; ind < num_elements; ind++) {
@@ -1774,8 +1778,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
sched, task_type_self, task_subtype_bh_density, flags, 0, ci, NULL);
t_bh_swallow = scheduler_addtask(
sched, task_type_self, task_subtype_bh_swallow, flags, 0, ci, NULL);
- t_do_swallow = scheduler_addtask(
- sched, task_type_self, task_subtype_do_swallow, flags, 0, ci, NULL);
+ t_do_gas_swallow =
+ scheduler_addtask(sched, task_type_self,
+ task_subtype_do_gas_swallow, flags, 0, ci, NULL);
+ t_do_bh_swallow =
+ scheduler_addtask(sched, task_type_self, task_subtype_do_bh_swallow,
+ flags, 0, ci, NULL);
t_bh_feedback =
scheduler_addtask(sched, task_type_self, task_subtype_bh_feedback,
flags, 0, ci, NULL);
@@ -1793,7 +1801,8 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
if (with_black_holes && bcount_i > 0) {
engine_addlink(e, &ci->black_holes.density, t_bh_density);
engine_addlink(e, &ci->black_holes.swallow, t_bh_swallow);
- engine_addlink(e, &ci->black_holes.do_swallow, t_do_swallow);
+ engine_addlink(e, &ci->black_holes.do_gas_swallow, t_do_gas_swallow);
+ engine_addlink(e, &ci->black_holes.do_bh_swallow, t_do_bh_swallow);
engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
}
@@ -1851,13 +1860,20 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
scheduler_addunlock(sched, t_bh_swallow,
ci->hydro.super->black_holes.swallow_ghost[0]);
- scheduler_addunlock(
- sched, ci->hydro.super->black_holes.swallow_ghost[0], t_do_swallow);
- scheduler_addunlock(sched, t_do_swallow,
+ scheduler_addunlock(sched,
+ ci->hydro.super->black_holes.swallow_ghost[0],
+ t_do_gas_swallow);
+ scheduler_addunlock(sched, t_do_gas_swallow,
ci->hydro.super->black_holes.swallow_ghost[1]);
scheduler_addunlock(sched,
ci->hydro.super->black_holes.swallow_ghost[1],
+ t_do_bh_swallow);
+ scheduler_addunlock(sched, t_do_bh_swallow,
+ ci->hydro.super->black_holes.swallow_ghost[2]);
+
+ scheduler_addunlock(sched,
+ ci->hydro.super->black_holes.swallow_ghost[2],
t_bh_feedback);
scheduler_addunlock(sched, t_bh_feedback,
ci->hydro.super->black_holes.black_holes_out);
@@ -1916,8 +1932,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
sched, task_type_pair, task_subtype_bh_density, flags, 0, ci, cj);
t_bh_swallow = scheduler_addtask(
sched, task_type_pair, task_subtype_bh_swallow, flags, 0, ci, cj);
- t_do_swallow = scheduler_addtask(
- sched, task_type_pair, task_subtype_do_swallow, flags, 0, ci, cj);
+ t_do_gas_swallow =
+ scheduler_addtask(sched, task_type_pair,
+ task_subtype_do_gas_swallow, flags, 0, ci, cj);
+ t_do_bh_swallow =
+ scheduler_addtask(sched, task_type_pair, task_subtype_do_bh_swallow,
+ flags, 0, ci, cj);
t_bh_feedback = scheduler_addtask(
sched, task_type_pair, task_subtype_bh_feedback, flags, 0, ci, cj);
}
@@ -1939,8 +1959,10 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
engine_addlink(e, &cj->black_holes.density, t_bh_density);
engine_addlink(e, &ci->black_holes.swallow, t_bh_swallow);
engine_addlink(e, &cj->black_holes.swallow, t_bh_swallow);
- engine_addlink(e, &ci->black_holes.do_swallow, t_do_swallow);
- engine_addlink(e, &cj->black_holes.do_swallow, t_do_swallow);
+ engine_addlink(e, &ci->black_holes.do_gas_swallow, t_do_gas_swallow);
+ engine_addlink(e, &cj->black_holes.do_gas_swallow, t_do_gas_swallow);
+ engine_addlink(e, &ci->black_holes.do_bh_swallow, t_do_bh_swallow);
+ engine_addlink(e, &cj->black_holes.do_bh_swallow, t_do_bh_swallow);
engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
engine_addlink(e, &cj->black_holes.feedback, t_bh_feedback);
}
@@ -2030,12 +2052,18 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
scheduler_addunlock(sched,
ci->hydro.super->black_holes.swallow_ghost[0],
- t_do_swallow);
- scheduler_addunlock(sched, t_do_swallow,
+ t_do_gas_swallow);
+ scheduler_addunlock(sched, t_do_gas_swallow,
ci->hydro.super->black_holes.swallow_ghost[1]);
scheduler_addunlock(sched,
ci->hydro.super->black_holes.swallow_ghost[1],
+ t_do_bh_swallow);
+ scheduler_addunlock(sched, t_do_bh_swallow,
+ ci->hydro.super->black_holes.swallow_ghost[2]);
+
+ scheduler_addunlock(sched,
+ ci->hydro.super->black_holes.swallow_ghost[2],
t_bh_feedback);
scheduler_addunlock(sched, t_bh_feedback,
ci->hydro.super->black_holes.black_holes_out);
@@ -2102,12 +2130,18 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
scheduler_addunlock(sched,
cj->hydro.super->black_holes.swallow_ghost[0],
- t_do_swallow);
- scheduler_addunlock(sched, t_do_swallow,
+ t_do_gas_swallow);
+ scheduler_addunlock(sched, t_do_gas_swallow,
cj->hydro.super->black_holes.swallow_ghost[1]);
scheduler_addunlock(sched,
cj->hydro.super->black_holes.swallow_ghost[1],
+ t_do_bh_swallow);
+ scheduler_addunlock(sched, t_do_bh_swallow,
+ cj->hydro.super->black_holes.swallow_ghost[2]);
+
+ scheduler_addunlock(sched,
+ cj->hydro.super->black_holes.swallow_ghost[2],
t_bh_feedback);
scheduler_addunlock(sched, t_bh_feedback,
cj->hydro.super->black_holes.black_holes_out);
@@ -2172,9 +2206,13 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
scheduler_addtask(sched, task_type_sub_self,
task_subtype_bh_swallow, flags, 0, ci, NULL);
- t_do_swallow =
+ t_do_gas_swallow =
+ scheduler_addtask(sched, task_type_sub_self,
+ task_subtype_do_gas_swallow, flags, 0, ci, NULL);
+
+ t_do_bh_swallow =
scheduler_addtask(sched, task_type_sub_self,
- task_subtype_do_swallow, flags, 0, ci, NULL);
+ task_subtype_do_bh_swallow, flags, 0, ci, NULL);
t_bh_feedback =
scheduler_addtask(sched, task_type_sub_self,
@@ -2193,7 +2231,8 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
if (with_black_holes && bcount_i > 0) {
engine_addlink(e, &ci->black_holes.density, t_bh_density);
engine_addlink(e, &ci->black_holes.swallow, t_bh_swallow);
- engine_addlink(e, &ci->black_holes.do_swallow, t_do_swallow);
+ engine_addlink(e, &ci->black_holes.do_gas_swallow, t_do_gas_swallow);
+ engine_addlink(e, &ci->black_holes.do_bh_swallow, t_do_bh_swallow);
engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
}
@@ -2257,13 +2296,20 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
scheduler_addunlock(sched, t_bh_swallow,
ci->hydro.super->black_holes.swallow_ghost[0]);
- scheduler_addunlock(
- sched, ci->hydro.super->black_holes.swallow_ghost[0], t_do_swallow);
- scheduler_addunlock(sched, t_do_swallow,
+ scheduler_addunlock(sched,
+ ci->hydro.super->black_holes.swallow_ghost[0],
+ t_do_gas_swallow);
+ scheduler_addunlock(sched, t_do_gas_swallow,
ci->hydro.super->black_holes.swallow_ghost[1]);
scheduler_addunlock(sched,
ci->hydro.super->black_holes.swallow_ghost[1],
+ t_do_bh_swallow);
+ scheduler_addunlock(sched, t_do_bh_swallow,
+ ci->hydro.super->black_holes.swallow_ghost[2]);
+
+ scheduler_addunlock(sched,
+ ci->hydro.super->black_holes.swallow_ghost[2],
t_bh_feedback);
scheduler_addunlock(sched, t_bh_feedback,
ci->hydro.super->black_holes.black_holes_out);
@@ -2326,9 +2372,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
t_bh_swallow =
scheduler_addtask(sched, task_type_sub_pair,
task_subtype_bh_swallow, flags, 0, ci, cj);
- t_do_swallow =
+ t_do_gas_swallow =
scheduler_addtask(sched, task_type_sub_pair,
- task_subtype_do_swallow, flags, 0, ci, cj);
+ task_subtype_do_gas_swallow, flags, 0, ci, cj);
+ t_do_bh_swallow =
+ scheduler_addtask(sched, task_type_sub_pair,
+ task_subtype_do_bh_swallow, flags, 0, ci, cj);
t_bh_feedback =
scheduler_addtask(sched, task_type_sub_pair,
task_subtype_bh_feedback, flags, 0, ci, cj);
@@ -2351,8 +2400,10 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
engine_addlink(e, &cj->black_holes.density, t_bh_density);
engine_addlink(e, &ci->black_holes.swallow, t_bh_swallow);
engine_addlink(e, &cj->black_holes.swallow, t_bh_swallow);
- engine_addlink(e, &ci->black_holes.do_swallow, t_do_swallow);
- engine_addlink(e, &cj->black_holes.do_swallow, t_do_swallow);
+ engine_addlink(e, &ci->black_holes.do_gas_swallow, t_do_gas_swallow);
+ engine_addlink(e, &cj->black_holes.do_gas_swallow, t_do_gas_swallow);
+ engine_addlink(e, &ci->black_holes.do_bh_swallow, t_do_bh_swallow);
+ engine_addlink(e, &cj->black_holes.do_bh_swallow, t_do_bh_swallow);
engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
engine_addlink(e, &cj->black_holes.feedback, t_bh_feedback);
}
@@ -2441,12 +2492,18 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
scheduler_addunlock(sched,
ci->hydro.super->black_holes.swallow_ghost[0],
- t_do_swallow);
- scheduler_addunlock(sched, t_do_swallow,
+ t_do_gas_swallow);
+ scheduler_addunlock(sched, t_do_gas_swallow,
ci->hydro.super->black_holes.swallow_ghost[1]);
scheduler_addunlock(sched,
ci->hydro.super->black_holes.swallow_ghost[1],
+ t_do_bh_swallow);
+ scheduler_addunlock(sched, t_do_bh_swallow,
+ ci->hydro.super->black_holes.swallow_ghost[2]);
+
+ scheduler_addunlock(sched,
+ ci->hydro.super->black_holes.swallow_ghost[2],
t_bh_feedback);
scheduler_addunlock(sched, t_bh_feedback,
ci->hydro.super->black_holes.black_holes_out);
@@ -2515,12 +2572,18 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
scheduler_addunlock(sched,
cj->hydro.super->black_holes.swallow_ghost[0],
- t_do_swallow);
- scheduler_addunlock(sched, t_do_swallow,
+ t_do_gas_swallow);
+ scheduler_addunlock(sched, t_do_gas_swallow,
cj->hydro.super->black_holes.swallow_ghost[1]);
scheduler_addunlock(sched,
cj->hydro.super->black_holes.swallow_ghost[1],
+ t_do_bh_swallow);
+ scheduler_addunlock(sched, t_do_bh_swallow,
+ cj->hydro.super->black_holes.swallow_ghost[2]);
+
+ scheduler_addunlock(sched,
+ cj->hydro.super->black_holes.swallow_ghost[2],
t_bh_feedback);
scheduler_addunlock(sched, t_bh_feedback,
cj->hydro.super->black_holes.black_holes_out);
diff --git a/src/engine_marktasks.c b/src/engine_marktasks.c
index 73cb820d6b221940b22e2b7ef37ad613f51d8ff8..a776684e8651e4bbbdd5f5ea965631d131636769 100644
--- a/src/engine_marktasks.c
+++ b/src/engine_marktasks.c
@@ -207,14 +207,26 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
}
else if (t_type == task_type_self &&
- t_subtype == task_subtype_do_swallow) {
+ t_subtype == task_subtype_do_gas_swallow) {
if (ci_active_black_holes) {
scheduler_activate(s, t);
}
}
else if (t_type == task_type_sub_self &&
- t_subtype == task_subtype_do_swallow) {
+ t_subtype == task_subtype_do_gas_swallow) {
+ if (ci_active_black_holes) scheduler_activate(s, t);
+ }
+
+ else if (t_type == task_type_self &&
+ t_subtype == task_subtype_do_bh_swallow) {
+ if (ci_active_black_holes) {
+ scheduler_activate(s, t);
+ }
+ }
+
+ else if (t_type == task_type_sub_self &&
+ t_subtype == task_subtype_do_bh_swallow) {
if (ci_active_black_holes) scheduler_activate(s, t);
}
@@ -404,7 +416,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
/* Black_Holes density */
else if ((t_subtype == task_subtype_bh_density ||
t_subtype == task_subtype_bh_swallow ||
- t_subtype == task_subtype_do_swallow ||
+ t_subtype == task_subtype_do_gas_swallow ||
t_subtype == task_subtype_bh_feedback) &&
(ci_active_black_holes || cj_active_black_holes) &&
(ci_nodeID == nodeID || cj_nodeID == nodeID)) {
@@ -892,7 +904,8 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
/* Black hole ghost tasks ? */
else if (t_type == task_type_bh_density_ghost ||
t_type == task_type_bh_swallow_ghost1 ||
- t_type == task_type_bh_swallow_ghost2) {
+ t_type == task_type_bh_swallow_ghost2 ||
+ t_type == task_type_bh_swallow_ghost3) {
if (cell_is_active_black_holes(t->ci, e)) scheduler_activate(s, t);
}
diff --git a/src/runner.c b/src/runner.c
index 1ab2db48b650ceccf27d6dafec920ed5ecce15ca..83ad8af827326d27a5d44bac467f1c9ec31e80d0 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -3198,20 +3198,23 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
} else { /* gpart is inactive */
/* Count the number of inhibited particles */
- if (gpart_is_inhibited(gp, e)) g_inhibited++;
+ if (gpart_is_inhibited(gp, e)) {
+ g_inhibited++;
+ } else {
- const integertime_t ti_end =
- get_integer_time_end(ti_current, gp->time_bin);
+ const integertime_t ti_end =
+ get_integer_time_end(ti_current, gp->time_bin);
- /* What is the next sync-point ? */
- ti_gravity_end_min = min(ti_end, ti_gravity_end_min);
- ti_gravity_end_max = max(ti_end, ti_gravity_end_max);
+ /* What is the next sync-point ? */
+ ti_gravity_end_min = min(ti_end, ti_gravity_end_min);
+ ti_gravity_end_max = max(ti_end, ti_gravity_end_max);
- const integertime_t ti_beg =
- get_integer_time_begin(ti_current + 1, gp->time_bin);
+ const integertime_t ti_beg =
+ get_integer_time_begin(ti_current + 1, gp->time_bin);
- /* What is the next starting point for this cell ? */
- ti_gravity_beg_max = max(ti_beg, ti_gravity_beg_max);
+ /* What is the next starting point for this cell ? */
+ ti_gravity_beg_max = max(ti_beg, ti_gravity_beg_max);
+ }
}
}
}
@@ -3319,22 +3322,26 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
} else {
/* Count the number of inhibited particles */
- if (bpart_is_inhibited(bp, e)) ++b_inhibited;
+ if (bpart_is_inhibited(bp, e)) {
+ message("BH inhibited!");
+ ++b_inhibited;
+ } else {
- const integertime_t ti_end =
- get_integer_time_end(ti_current, bp->time_bin);
+ const integertime_t ti_end =
+ get_integer_time_end(ti_current, bp->time_bin);
- const integertime_t ti_beg =
- get_integer_time_begin(ti_current + 1, bp->time_bin);
+ const integertime_t ti_beg =
+ get_integer_time_begin(ti_current + 1, bp->time_bin);
- ti_black_holes_end_min = min(ti_end, ti_black_holes_end_min);
- ti_black_holes_end_max = max(ti_end, ti_black_holes_end_max);
- ti_gravity_end_min = min(ti_end, ti_gravity_end_min);
- ti_gravity_end_max = max(ti_end, ti_gravity_end_max);
+ ti_black_holes_end_min = min(ti_end, ti_black_holes_end_min);
+ ti_black_holes_end_max = max(ti_end, ti_black_holes_end_max);
+ ti_gravity_end_min = min(ti_end, ti_gravity_end_min);
+ ti_gravity_end_max = max(ti_end, ti_gravity_end_max);
- /* What is the next starting point for this cell ? */
- ti_black_holes_beg_max = max(ti_beg, ti_black_holes_beg_max);
- ti_gravity_beg_max = max(ti_beg, ti_gravity_beg_max);
+ /* What is the next starting point for this cell ? */
+ ti_black_holes_beg_max = max(ti_beg, ti_black_holes_beg_max);
+ ti_gravity_beg_max = max(ti_beg, ti_gravity_beg_max);
+ }
}
}
@@ -3731,7 +3738,7 @@ void runner_do_end_grav_force(struct runner *r, struct cell *c, int timer) {
* @param c The #cell.
* @param timer Are we timing this?
*/
-void runner_do_swallow(struct runner *r, struct cell *c, int timer) {
+void runner_do_gas_swallow(struct runner *r, struct cell *c, int timer) {
struct engine *e = r->e;
struct space *s = e->s;
@@ -3759,7 +3766,7 @@ void runner_do_swallow(struct runner *r, struct cell *c, int timer) {
if (c->progeny[k] != NULL) {
struct cell *restrict cp = c->progeny[k];
- runner_do_swallow(r, cp, 0);
+ runner_do_gas_swallow(r, cp, 0);
}
}
} else {
@@ -3883,7 +3890,7 @@ void runner_do_swallow(struct runner *r, struct cell *c, int timer) {
* @param c The #cell.
* @param timer Are we timing this?
*/
-void runner_do_swallow_self(struct runner *r, struct cell *c, int timer) {
+void runner_do_gas_swallow_self(struct runner *r, struct cell *c, int timer) {
#ifdef SWIFT_DEBUG_CHECKS
if (c->nodeID != r->e->nodeID) error("Running self task on foreign node");
@@ -3891,7 +3898,7 @@ void runner_do_swallow_self(struct runner *r, struct cell *c, int timer) {
error("Running self task on inactive cell");
#endif
- runner_do_swallow(r, c, timer);
+ runner_do_gas_swallow(r, c, timer);
}
/**
@@ -3902,8 +3909,8 @@ void runner_do_swallow_self(struct runner *r, struct cell *c, int timer) {
* @param cj Second #cell.
* @param timer Are we timing this?
*/
-void runner_do_swallow_pair(struct runner *r, struct cell *ci, struct cell *cj,
- int timer) {
+void runner_do_gas_swallow_pair(struct runner *r, struct cell *ci,
+ struct cell *cj, int timer) {
const struct engine *e = r->e;
@@ -3916,8 +3923,213 @@ void runner_do_swallow_pair(struct runner *r, struct cell *ci, struct cell *cj,
/* Run the swallowing loop only in the cell that is the neighbour of the
* active BH */
- if (cell_is_active_black_holes(cj, e)) runner_do_swallow(r, ci, timer);
- if (cell_is_active_black_holes(ci, e)) runner_do_swallow(r, cj, timer);
+ if (cell_is_active_black_holes(cj, e)) runner_do_gas_swallow(r, ci, timer);
+ if (cell_is_active_black_holes(ci, e)) runner_do_gas_swallow(r, cj, timer);
+}
+
+/**
+ * @brief Process all the BH particles in a cell that have been flagged for
+ * swallowing by a black hole.
+ *
+ * This is done by recursing down to the leaf-level and skipping the sub-cells
+ * that have not been drifted as they would not have any particles with
+ * swallowing flag. We then loop over the particles with a flag and look into
+ * the space-wide list of black holes for the particle with the corresponding
+ * ID. If found, the BH swallows the BH particle and the BH particle is
+ * removed. If the cell is local, we may be looking for a foreign BH, in which
+ * case, we do not update the BH (that will be done on its node) but just remove
+ * the BH particle.
+ *
+ * @param r The thread #runner.
+ * @param c The #cell.
+ * @param timer Are we timing this?
+ */
+void runner_do_bh_swallow(struct runner *r, struct cell *c, int timer) {
+
+ struct engine *e = r->e;
+ struct space *s = e->s;
+ struct bpart *bparts = s->bparts;
+ const size_t nr_bpart = s->nr_bparts;
+#ifdef WITH_MPI
+ struct bpart *bparts_foreign = s->bparts_foreign;
+ const size_t nr_bparts_foreign = s->nr_bparts_foreign;
+#endif
+
+ struct bpart *cell_bparts = c->black_holes.parts;
+
+ /* Early abort?
+ * (We only want cells for which we drifted the gas as these are
+ * the only ones that could have gas particles that have been flagged
+ * for swallowing) */
+ if (c->black_holes.count == 0 ||
+ c->black_holes.ti_old_part != e->ti_current) {
+ return;
+ }
+
+ /* Loop over the progeny ? */
+ if (c->split) {
+ for (int k = 0; k < 8; k++) {
+ if (c->progeny[k] != NULL) {
+ struct cell *restrict cp = c->progeny[k];
+
+ runner_do_bh_swallow(r, cp, 0);
+ }
+ }
+ } else {
+
+ /* Loop over all the gas particles in the cell
+ * Note that the cell (and hence the parts) may be local or foreign. */
+ const size_t nr_cell_bparts = c->black_holes.count;
+ for (size_t k = 0; k < nr_cell_bparts; k++) {
+
+ /* Get a handle on the part. */
+ struct bpart *const cell_bp = &cell_bparts[k];
+
+ /* Ignore inhibited particles (they have already been removed!) */
+ if (bpart_is_inhibited(cell_bp, e)) continue;
+
+ /* Get the ID of the black holes that will swallow this part */
+ const long long swallow_id =
+ black_holes_get_bpart_swallow_id(&cell_bp->merger_data);
+
+ /* Has this particle been flagged for swallowing? */
+ if (swallow_id >= 0) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (cell_bp->ti_drift != e->ti_current)
+ error("Trying to swallow an un-drifted particle.");
+#endif
+
+ /* ID of the BH swallowing this particle */
+ const long long BH_id = swallow_id;
+
+ /* Have we found this particle's BH already? */
+ int found = 0;
+
+ /* Let's look for the hungry black hole in the local list */
+ for (size_t i = 0; i < nr_bpart; ++i) {
+
+ /* Get a handle on the bpart. */
+ struct bpart *bp = &bparts[i];
+
+ if (bp->id == BH_id) {
+
+ /* Lock the space as we are going to work directly on the bpart list
+ */
+ lock_lock(&s->lock);
+
+ /* Swallow the gas particle (i.e. update the BH properties) */
+ // black_holes_swallow_part(bp, p, xp, e->cosmology);
+
+ /* Release the space as we are done updating the bpart */
+ if (lock_unlock(&s->lock) != 0)
+ error("Failed to unlock the space.");
+
+ message("BH %lld swallowing particle %lld", bp->id, cell_bp->id);
+
+ /* If the gas particle is local, remove it */
+ if (c->nodeID == e->nodeID) {
+
+ message("BH %lld removing particle %lld", bp->id, cell_bp->id);
+
+ /* Finally, remove the gas particle from the system */
+ struct gpart *gp = cell_bp->gpart;
+ cell_remove_bpart(e, c, cell_bp);
+ cell_remove_gpart(e, c, gp);
+ }
+
+ /* In any case, prevent the particle from being re-swallowed */
+ black_holes_mark_bpart_as_merged(&cell_bp->merger_data);
+
+ found = 1;
+ break;
+ }
+
+ } /* Loop over local BHs */
+
+#ifdef WITH_MPI
+
+ /* We could also be in the case of a local gas particle being
+ * swallowed by a foreign BH. In this case, we won't update the
+ * BH but just remove the particle from the local list. */
+ if (c->nodeID == e->nodeID && !found) {
+
+ /* Let's look for the foreign hungry black hole */
+ for (size_t i = 0; i < nr_bparts_foreign; ++i) {
+
+ /* Get a handle on the bpart. */
+ struct bpart *bp = &bparts_foreign[i];
+
+ if (bp->id == BH_id) {
+
+ message("BH %lld removing particle %lld (foreign BH case)",
+ bp->id, cell_bp->id);
+
+ /* Finally, remove the gas particle from the system */
+ // struct gpart *gp = p->gpart;
+ // cell_remove_part(e, c, p, xp);
+ // cell_remove_gpart(e, c, gp);
+
+ found = 1;
+ break;
+ }
+ } /* Loop over foreign BHs */
+ } /* Is the cell local? */
+#endif
+
+ /* If we have a local particle, we must have found the BH in one
+ * of our list of black holes. */
+ if (c->nodeID == e->nodeID && !found) {
+ error("BH particle %lld could not find BH %lld to be swallowed",
+ cell_bp->id, swallow_id);
+ }
+ } /* Part was flagged for swallowing */
+ } /* Loop over the parts */
+ } /* Cell is not split */
+}
+
+/**
+ * @brief Processing of bh particles to swallow - self task case.
+ *
+ * @param r The thread #runner.
+ * @param c The #cell.
+ * @param timer Are we timing this?
+ */
+void runner_do_bh_swallow_self(struct runner *r, struct cell *c, int timer) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->nodeID != r->e->nodeID) error("Running self task on foreign node");
+ if (!cell_is_active_black_holes(c, r->e))
+ error("Running self task on inactive cell");
+#endif
+
+ runner_do_bh_swallow(r, c, timer);
+}
+
+/**
+ * @brief Processing of bh particles to swallow - pair task case.
+ *
+ * @param r The thread #runner.
+ * @param ci First #cell.
+ * @param cj Second #cell.
+ * @param timer Are we timing this?
+ */
+void runner_do_bh_swallow_pair(struct runner *r, struct cell *ci,
+ struct cell *cj, int timer) {
+
+ const struct engine *e = r->e;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ci->nodeID != e->nodeID && cj->nodeID != e->nodeID)
+ error("Running pair task on foreign node");
+ if (!cell_is_active_black_holes(ci, e) && !cell_is_active_black_holes(cj, e))
+ error("Running pair task with two inactive cells");
+#endif
+
+ /* Run the swallowing loop only in the cell that is the neighbour of the
+ * active BH */
+ if (cell_is_active_black_holes(cj, e)) runner_do_bh_swallow(r, ci, timer);
+ if (cell_is_active_black_holes(ci, e)) runner_do_bh_swallow(r, cj, timer);
}
/**
@@ -4340,12 +4552,15 @@ void *runner_main(void *data) {
runner_doself_branch_bh_density(r, ci);
else if (t->subtype == task_subtype_bh_swallow)
runner_doself_branch_bh_swallow(r, ci);
- else if (t->subtype == task_subtype_do_swallow)
- runner_do_swallow_self(r, ci, 1);
+ else if (t->subtype == task_subtype_do_gas_swallow)
+ runner_do_gas_swallow_self(r, ci, 1);
+ else if (t->subtype == task_subtype_do_bh_swallow)
+ runner_do_bh_swallow_self(r, ci, 1);
else if (t->subtype == task_subtype_bh_feedback)
runner_doself_branch_bh_feedback(r, ci);
else
- error("Unknown/invalid task subtype (%d).", t->subtype);
+ error("Unknown/invalid task subtype (%s).",
+ subtaskID_names[t->subtype]);
break;
case task_type_pair:
@@ -4369,8 +4584,8 @@ void *runner_main(void *data) {
runner_dopair_branch_bh_density(r, ci, cj);
else if (t->subtype == task_subtype_bh_swallow)
runner_dopair_branch_bh_swallow(r, ci, cj);
- else if (t->subtype == task_subtype_do_swallow) {
- runner_do_swallow_pair(r, ci, cj, 1);
+ else if (t->subtype == task_subtype_do_gas_swallow) {
+ runner_do_gas_swallow_pair(r, ci, cj, 1);
} else if (t->subtype == task_subtype_bh_feedback)
runner_dopair_branch_bh_feedback(r, ci, cj);
else
@@ -4396,8 +4611,8 @@ void *runner_main(void *data) {
runner_dosub_self_bh_density(r, ci, 1);
else if (t->subtype == task_subtype_bh_swallow)
runner_dosub_self_bh_swallow(r, ci, 1);
- else if (t->subtype == task_subtype_do_swallow)
- runner_do_swallow_self(r, ci, 1);
+ else if (t->subtype == task_subtype_do_gas_swallow)
+ runner_do_gas_swallow_self(r, ci, 1);
else if (t->subtype == task_subtype_bh_feedback)
runner_dosub_self_bh_feedback(r, ci, 1);
else
@@ -4423,8 +4638,8 @@ void *runner_main(void *data) {
runner_dosub_pair_bh_density(r, ci, cj, 1);
else if (t->subtype == task_subtype_bh_swallow)
runner_dosub_pair_bh_swallow(r, ci, cj, 1);
- else if (t->subtype == task_subtype_do_swallow) {
- runner_do_swallow_pair(r, ci, cj, 1);
+ else if (t->subtype == task_subtype_do_gas_swallow) {
+ runner_do_gas_swallow_pair(r, ci, cj, 1);
} else if (t->subtype == task_subtype_bh_feedback)
runner_dosub_pair_bh_feedback(r, ci, cj, 1);
else
@@ -4464,7 +4679,7 @@ void *runner_main(void *data) {
case task_type_bh_density_ghost:
runner_do_black_holes_density_ghost(r, ci, 1);
break;
- case task_type_bh_swallow_ghost2:
+ case task_type_bh_swallow_ghost3:
runner_do_black_holes_swallow_ghost(r, ci, 1);
break;
case task_type_drift_part:
diff --git a/src/space.c b/src/space.c
index d6dbbc22adda0d913f2d8ffcff46d6e5f99a6ad0..f622fda5f6e47205aeed2ce1f108848a81d2e00c 100644
--- a/src/space.c
+++ b/src/space.c
@@ -234,9 +234,11 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->black_holes.density_ghost = NULL;
c->black_holes.swallow_ghost[0] = NULL;
c->black_holes.swallow_ghost[1] = NULL;
+ c->black_holes.swallow_ghost[2] = NULL;
c->black_holes.density = NULL;
c->black_holes.swallow = NULL;
- c->black_holes.do_swallow = NULL;
+ c->black_holes.do_gas_swallow = NULL;
+ c->black_holes.do_bh_swallow = NULL;
c->black_holes.feedback = NULL;
c->kick1 = NULL;
c->kick2 = NULL;
diff --git a/src/task.c b/src/task.c
index 9d7167e0b82cf8459676f5216aabb1c0f6ce7fb8..20b5447b097dc9ea8c89e7140b6a6cc9b36ff64a 100644
--- a/src/task.c
+++ b/src/task.c
@@ -94,18 +94,19 @@ const char *taskID_names[task_type_count] = {"none",
"bh_density_ghost",
"bh_swallow_ghost1",
"bh_swallow_ghost2",
+ "bh_swallow_ghost3",
"fof_self",
"fof_pair"};
/* Sub-task type names. */
const char *subtaskID_names[task_subtype_count] = {
- "none", "density", "gradient", "force",
- "limiter", "grav", "external_grav", "tend_part",
- "tend_gpart", "tend_spart", "tend_bpart", "xv",
- "rho", "part_swallow", "gpart", "multipole",
- "spart", "stars_density", "stars_feedback", "sf_count",
- "bpart_rho", "bpart_swallow", "bpart_feedback", "bh_density",
- "bh_swallow", "do_swallow", "bh_feedback"};
+ "none", "density", "gradient", "force",
+ "limiter", "grav", "external_grav", "tend_part",
+ "tend_gpart", "tend_spart", "tend_bpart", "xv",
+ "rho", "part_swallow", "gpart", "multipole",
+ "spart", "stars_density", "stars_feedback", "sf_count",
+ "bpart_rho", "bpart_swallow", "bpart_feedback", "bh_density",
+ "bh_swallow", "do_gas_swallow", "do_bh_swallow", "bh_feedback"};
#ifdef WITH_MPI
/* MPI communicators for the subtypes. */
@@ -177,7 +178,7 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
case task_type_drift_bpart:
case task_type_bh_density_ghost:
- case task_type_bh_swallow_ghost2:
+ case task_type_bh_swallow_ghost3:
return task_action_bpart;
break;
@@ -202,10 +203,14 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
case task_subtype_bh_density:
case task_subtype_bh_feedback:
case task_subtype_bh_swallow:
- case task_subtype_do_swallow:
+ case task_subtype_do_gas_swallow:
return task_action_all;
break;
+ case task_subtype_do_bh_swallow:
+ return task_action_bpart;
+ break;
+
case task_subtype_grav:
case task_subtype_external_grav:
return task_action_gpart;
@@ -446,9 +451,11 @@ void task_unlock(struct task *t) {
} else if ((subtype == task_subtype_bh_density) ||
(subtype == task_subtype_bh_feedback) ||
(subtype == task_subtype_bh_swallow) ||
- (subtype == task_subtype_do_swallow)) {
+ (subtype == task_subtype_do_gas_swallow)) {
cell_bunlocktree(ci);
cell_unlocktree(ci);
+ } else if (subtype == task_subtype_do_bh_swallow) {
+ cell_bunlocktree(ci);
} else {
cell_unlocktree(ci);
}
@@ -470,11 +477,14 @@ void task_unlock(struct task *t) {
} else if ((subtype == task_subtype_bh_density) ||
(subtype == task_subtype_bh_feedback) ||
(subtype == task_subtype_bh_swallow) ||
- (subtype == task_subtype_do_swallow)) {
+ (subtype == task_subtype_do_gas_swallow)) {
cell_bunlocktree(ci);
cell_bunlocktree(cj);
cell_unlocktree(ci);
cell_unlocktree(cj);
+ } else if (subtype == task_subtype_do_bh_swallow) {
+ cell_bunlocktree(ci);
+ cell_bunlocktree(cj);
} else {
cell_unlocktree(ci);
cell_unlocktree(cj);
@@ -599,7 +609,7 @@ int task_lock(struct task *t) {
} else if ((subtype == task_subtype_bh_density) ||
(subtype == task_subtype_bh_feedback) ||
(subtype == task_subtype_bh_swallow) ||
- (subtype == task_subtype_do_swallow)) {
+ (subtype == task_subtype_do_gas_swallow)) {
if (ci->black_holes.hold) return 0;
if (ci->hydro.hold) return 0;
if (cell_blocktree(ci) != 0) return 0;
@@ -607,6 +617,9 @@ int task_lock(struct task *t) {
cell_bunlocktree(ci);
return 0;
}
+ } else if (subtype == task_subtype_do_bh_swallow) {
+ if (ci->black_holes.hold) return 0;
+ if (cell_blocktree(ci) != 0) return 0;
} else { /* subtype == hydro */
if (ci->hydro.hold) return 0;
if (cell_locktree(ci) != 0) return 0;
@@ -656,7 +669,7 @@ int task_lock(struct task *t) {
} else if ((subtype == task_subtype_bh_density) ||
(subtype == task_subtype_bh_feedback) ||
(subtype == task_subtype_bh_swallow) ||
- (subtype == task_subtype_do_swallow)) {
+ (subtype == task_subtype_do_gas_swallow)) {
/* Lock the BHs and the gas particles in both cells */
if (ci->black_holes.hold || cj->black_holes.hold) return 0;
if (ci->hydro.hold || cj->hydro.hold) return 0;
@@ -676,6 +689,13 @@ int task_lock(struct task *t) {
cell_unlocktree(ci);
return 0;
}
+ } else if (subtype == task_subtype_do_bh_swallow) {
+ if (ci->black_holes.hold || cj->black_holes.hold) return 0;
+ if (cell_blocktree(ci) != 0) return 0;
+ if (cell_blocktree(cj) != 0) {
+ cell_bunlocktree(ci);
+ return 0;
+ }
} else { /* subtype == hydro */
/* Lock the parts in both cells */
if (ci->hydro.hold || cj->hydro.hold) return 0;
@@ -799,8 +819,11 @@ void task_get_group_name(int type, int subtype, char *cluster) {
case task_subtype_bh_swallow:
strcpy(cluster, "BHSwallow");
break;
- case task_subtype_do_swallow:
- strcpy(cluster, "DoSwallow");
+ case task_subtype_do_gas_swallow:
+ strcpy(cluster, "DoGasSwallow");
+ break;
+ case task_subtype_do_bh_swallow:
+ strcpy(cluster, "DoBHSwallow");
break;
case task_subtype_bh_feedback:
strcpy(cluster, "BHFeedback");
diff --git a/src/task.h b/src/task.h
index 995f39357f810e068b3fe4e4304bfab8f9734080..a28e31c37bccb661c5a1ee496183c3652a428de4 100644
--- a/src/task.h
+++ b/src/task.h
@@ -87,8 +87,9 @@ enum task_types {
task_type_bh_in, /* Implicit */
task_type_bh_out, /* Implicit */
task_type_bh_density_ghost,
- task_type_bh_swallow_ghost1,
+ task_type_bh_swallow_ghost1, /* Implicit */
task_type_bh_swallow_ghost2,
+ task_type_bh_swallow_ghost3, /* Implicit */
task_type_fof_self,
task_type_fof_pair,
task_type_count
@@ -123,7 +124,8 @@ enum task_subtypes {
task_subtype_bpart_feedback,
task_subtype_bh_density,
task_subtype_bh_swallow,
- task_subtype_do_swallow,
+ task_subtype_do_gas_swallow,
+ task_subtype_do_bh_swallow,
task_subtype_bh_feedback,
task_subtype_count
} __attribute__((packed));