diff --git a/AUTHORS b/AUTHORS
index ff77660e38846a8dfec6bbbfc8cdd62e13d67114..21268f93b663ceaa1bc9326555ac98c30ca3c3fe 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -30,4 +30,5 @@ Orestis Karapiperis karapiperis@lorentz.leidenuniv.nl
Stan Verhoeve s06verhoeve@gmail.com
Nikyta Shchutskyi shchutskyi@lorentz.leidenuniv.nl
Will Roper w.roper@sussex.ac.uk
-Darwin Roduit darwin.roduit@alumni.epfl.ch
+Darwin Roduit darwin.roduit@alumni.epfl.ch
+Jonathan Davies j.j.davies@ljmu.ac.uk
diff --git a/doc/RTD/source/NewOption/sink_adding_new_scheme.rst b/doc/RTD/source/NewOption/sink_adding_new_scheme.rst
index b25829a18ae63f677df961982d607ee7b52c8cd5..023446e1fdab844ade130d53eedac08c19170cd1 100644
--- a/doc/RTD/source/NewOption/sink_adding_new_scheme.rst
+++ b/doc/RTD/source/NewOption/sink_adding_new_scheme.rst
@@ -22,6 +22,19 @@ Before forming a sink, the code loops over all gas and sink particles to gather
Then, to decide if we can turn a gas particle into a sink particle, the function ``sink_is_forming()`` is called. Before forming a sink particle, there is a call to ``sink_should_convert_to_sink()``. This function determines whether the gas particle must transform into a sink. Both functions return either 0 or 1.
+
+Gas-sink density interactions
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The first interaction task to be run for the sinks is the density task. This task updates the smoothing length for the sink particle, unless a fixed cutoff radius is being used (coming soon). It can also calculate the contributions made by neigbouring gas particles to the density, sound speed, velocity etc. at the location of the sink. Code for these interactions should be added to ``sink_iact.h/runner_iact_nonsym_sinks_gas_density()``.
+
+Once the contributions of all neigbouring gas particles have been calculated, the density calculation is completed by the sink density ghost task. You can set what this task does with the functions ``sink_end_density()`` and ``sink_prepare_swallow()`` in ``sink.h``.
+
+The ``sink_end_density()`` function completes the calculation of the smoothing length (coming soon), and this is where you can finish density-based calculations by e.g. dividing mass-weighted contributions to the velocity field by the total density in the kernel. For examples of this, see the equivalent task for the black hole particles.
+
+The ``sink_prepare_swallow()`` task is where you can calculate density-based quantities that you might need to use in swallowing interactions later. For example, a Bondi-Hoyle accretion prescription should calculate an accretion rate and target mass to be accreted here.
+
+
Gas and sink flagging: finding whom to eat
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
diff --git a/src/cell.c b/src/cell.c
index a610135c87f479f5cc4a7eff52dc230d2286b47e..4f58054830f31a43dd96f8539262af0fc0b5ceb6 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -588,6 +588,7 @@ void cell_clean_links(struct cell *c, void *data) {
c->stars.prepare2 = NULL;
c->stars.feedback = NULL;
c->sinks.swallow = NULL;
+ c->sinks.density = NULL;
c->sinks.do_sink_swallow = NULL;
c->sinks.do_gas_swallow = NULL;
c->black_holes.density = NULL;
diff --git a/src/cell_sinks.h b/src/cell_sinks.h
index ac1d9d592030e83df9c8c48b8a0772e2981c201f..6bd2f9a9afed9847001a9a106b6f2dac9a1b433d 100644
--- a/src/cell_sinks.h
+++ b/src/cell_sinks.h
@@ -57,6 +57,12 @@ struct cell_sinks {
/*! Implicit tasks marking the entry of the sink block of tasks */
struct task *sink_in;
+ /*! The sink ghost task itself */
+ struct task *density_ghost;
+
+ /*! Linked list of the tasks computing this cell's sink density. */
+ struct link *density;
+
/*! Implicit tasks marking the end of sink swallow */
struct task *sink_ghost1;
diff --git a/src/cell_unskip.c b/src/cell_unskip.c
index 75c51f3f1347dcd385066f1c22a47b12d0837fe0..3ce7d791a712e5dcacbc7794af08e81513fa5fdd 100644
--- a/src/cell_unskip.c
+++ b/src/cell_unskip.c
@@ -2934,11 +2934,12 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
cell_activate_drift_sink(c, s);
}
- /* Un-skip the swallow tasks involved with this cell. */
- for (struct link *l = c->sinks.swallow; l != NULL; l = l->next) {
+ /* Un-skip the density tasks involved with this cell. */
+ for (struct link *l = c->sinks.density; l != NULL; l = l->next) {
struct task *t = l->t;
struct cell *ci = t->ci;
struct cell *cj = t->cj;
+
#ifdef WITH_MPI
const int ci_nodeID = ci->nodeID;
const int cj_nodeID = (cj != NULL) ? cj->nodeID : -1;
@@ -2949,7 +2950,6 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
const int ci_active =
cell_is_active_sinks(ci, e) || cell_is_active_hydro(ci, e);
-
const int cj_active = (cj != NULL) && (cell_is_active_sinks(cj, e) ||
cell_is_active_hydro(cj, e));
@@ -2959,21 +2959,17 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
scheduler_activate(s, t);
- /* Activate the drifts */
+ /* Activate drifts for self tasks */
if (t->type == task_type_self) {
cell_activate_drift_part(ci, s);
cell_activate_drift_sink(ci, s);
- if (with_timestep_sync) cell_activate_sync_part(ci, s);
}
- /* Activate the drifts */
+ /* Activate drifts for pair tasks */
else if (t->type == task_type_pair) {
-
- /* Activate the drift tasks. */
if (ci_nodeID == nodeID) cell_activate_drift_sink(ci, s);
if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
if (ci_nodeID == nodeID) cell_activate_sink_formation_tasks(ci->top, s);
-
if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
if (cj_nodeID == nodeID) cell_activate_drift_sink(cj, s);
if (cj_nodeID == nodeID) cell_activate_sink_formation_tasks(cj->top, s);
@@ -2997,7 +2993,6 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
* a pair task as to not miss any dependencies */
if (ci_nodeID == nodeID)
scheduler_activate(s, ci->hydro.super->sinks.sink_in);
-
if (cj_nodeID == nodeID)
scheduler_activate(s, cj->hydro.super->sinks.sink_in);
@@ -3012,6 +3007,32 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
}
}
+ /* Un-skip the swallow tasks involved with this cell. */
+ for (struct link *l = c->sinks.swallow; l != NULL; l = l->next) {
+ struct task *t = l->t;
+ struct cell *ci = t->ci;
+ struct cell *cj = t->cj;
+#ifdef WITH_MPI
+ const int ci_nodeID = ci->nodeID;
+ const int cj_nodeID = (cj != NULL) ? cj->nodeID : -1;
+#else
+ const int ci_nodeID = nodeID;
+ const int cj_nodeID = nodeID;
+#endif
+
+ const int ci_active =
+ cell_is_active_sinks(ci, e) || cell_is_active_hydro(ci, e);
+ const int cj_active = (cj != NULL) && (cell_is_active_sinks(cj, e) ||
+ cell_is_active_hydro(cj, e));
+
+ /* Only activate tasks that involve a local active cell. */
+ if ((ci_active || cj_active) &&
+ (ci_nodeID == nodeID || cj_nodeID == nodeID)) {
+
+ scheduler_activate(s, t);
+ }
+ }
+
/* Un-skip the do_sink_swallow tasks involved with this cell. */
for (struct link *l = c->sinks.do_sink_swallow; l != NULL; l = l->next) {
struct task *t = l->t;
@@ -3027,7 +3048,6 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
const int ci_active =
cell_is_active_sinks(ci, e) || cell_is_active_hydro(ci, e);
-
const int cj_active = (cj != NULL) && (cell_is_active_sinks(cj, e) ||
cell_is_active_hydro(cj, e));
@@ -3053,7 +3073,6 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
const int ci_active =
cell_is_active_sinks(ci, e) || cell_is_active_hydro(ci, e);
-
const int cj_active = (cj != NULL) && (cell_is_active_sinks(cj, e) ||
cell_is_active_hydro(cj, e));
@@ -3082,6 +3101,8 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
cell_activate_sink_formation_tasks(c->top, s);
cell_activate_super_sink_drifts(c->top, s);
}
+ if (c->sinks.density_ghost != NULL)
+ scheduler_activate(s, c->sinks.density_ghost);
if (c->sinks.sink_ghost1 != NULL)
scheduler_activate(s, c->sinks.sink_ghost1);
if (c->sinks.sink_ghost2 != NULL)
@@ -3100,7 +3121,6 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
if (c->csds != NULL) scheduler_activate(s, c->csds);
#endif
}
-
return rebuild;
}
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index 227cdfec9905b904c7f0209df79b4b2e1f77ca3f..5592e384ec9d4fc2a1609ad5744ea1b46397f8ad 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -1660,6 +1660,9 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c,
scheduler_addtask(s, task_type_sink_in, task_subtype_none, 0,
/* implicit = */ 1, c, NULL);
+ c->sinks.density_ghost = scheduler_addtask(
+ s, task_type_sink_density_ghost, task_subtype_none, 0, 0, c, NULL);
+
c->sinks.sink_ghost1 =
scheduler_addtask(s, task_type_sink_ghost1, task_subtype_none, 0,
/* implicit = */ 1, c, NULL);
@@ -2483,6 +2486,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
struct task *t_do_gas_swallow = NULL;
struct task *t_do_bh_swallow = NULL;
struct task *t_bh_feedback = NULL;
+ struct task *t_sink_density = NULL;
struct task *t_sink_swallow = NULL;
struct task *t_rt_gradient = NULL;
struct task *t_rt_transport = NULL;
@@ -2554,6 +2558,9 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
/* The sink tasks */
if (with_sink) {
+ t_sink_density =
+ scheduler_addtask(sched, task_type_self, task_subtype_sink_density,
+ flags, 0, ci, NULL);
t_sink_swallow =
scheduler_addtask(sched, task_type_self, task_subtype_sink_swallow,
flags, 0, ci, NULL);
@@ -2605,6 +2612,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
#endif
}
if (with_sink) {
+ engine_addlink(e, &ci->sinks.density, t_sink_density);
engine_addlink(e, &ci->sinks.swallow, t_sink_swallow);
engine_addlink(e, &ci->sinks.do_sink_swallow, t_sink_do_sink_swallow);
engine_addlink(e, &ci->sinks.do_gas_swallow, t_sink_do_gas_swallow);
@@ -2684,15 +2692,20 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
/* The sink's tasks. */
if (with_sink) {
- /* Do the sink_swallow */
- scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
- t_sink_swallow);
+ /* Sink density */
scheduler_addunlock(sched, ci->hydro.super->sinks.drift,
- t_sink_swallow);
- scheduler_addunlock(sched, ci->hydro.super->hydro.sorts,
- t_sink_swallow);
+ t_sink_density);
+ scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
+ t_sink_density);
scheduler_addunlock(sched, ci->hydro.super->sinks.sink_in,
+ t_sink_density);
+ scheduler_addunlock(sched, t_sink_density,
+ ci->hydro.super->sinks.density_ghost);
+
+ /* Do the sink_swallow */
+ scheduler_addunlock(sched, ci->hydro.super->sinks.density_ghost,
t_sink_swallow);
+
scheduler_addunlock(sched, t_sink_swallow,
ci->hydro.super->sinks.sink_ghost1);
@@ -2832,6 +2845,8 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
/* The sink tasks */
if (with_sink) {
+ t_sink_density = scheduler_addtask(
+ sched, task_type_pair, task_subtype_sink_density, flags, 0, ci, cj);
t_sink_swallow = scheduler_addtask(
sched, task_type_pair, task_subtype_sink_swallow, flags, 0, ci, cj);
t_sink_do_sink_swallow = scheduler_addtask(
@@ -2910,13 +2925,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
#endif
}
if (with_sink) {
- /* Formation */
+ engine_addlink(e, &ci->sinks.density, t_sink_density);
+ engine_addlink(e, &cj->sinks.density, t_sink_density);
engine_addlink(e, &ci->sinks.swallow, t_sink_swallow);
engine_addlink(e, &cj->sinks.swallow, t_sink_swallow);
- /* Merger */
engine_addlink(e, &ci->sinks.do_sink_swallow, t_sink_do_sink_swallow);
engine_addlink(e, &cj->sinks.do_sink_swallow, t_sink_do_sink_swallow);
- /* Accretion */
engine_addlink(e, &ci->sinks.do_gas_swallow, t_sink_do_gas_swallow);
engine_addlink(e, &cj->sinks.do_gas_swallow, t_sink_do_gas_swallow);
}
@@ -3039,14 +3053,18 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
if (with_sink) {
- /* Do the sink_swallow */
- scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
- t_sink_swallow);
+ /* Sink density */
scheduler_addunlock(sched, ci->hydro.super->sinks.drift,
- t_sink_swallow);
- scheduler_addunlock(sched, ci->hydro.super->hydro.sorts,
- t_sink_swallow);
+ t_sink_density);
+ scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
+ t_sink_density);
scheduler_addunlock(sched, ci->hydro.super->sinks.sink_in,
+ t_sink_density);
+ scheduler_addunlock(sched, t_sink_density,
+ ci->hydro.super->sinks.density_ghost);
+
+ /* Do the sink_swallow */
+ scheduler_addunlock(sched, ci->hydro.super->sinks.density_ghost,
t_sink_swallow);
scheduler_addunlock(sched, t_sink_swallow,
ci->hydro.super->sinks.sink_ghost1);
@@ -3194,14 +3212,18 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
if (with_sink) {
- /* Do the sink_swallow */
- scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
- t_sink_swallow);
+ /* Sink density */
scheduler_addunlock(sched, cj->hydro.super->sinks.drift,
- t_sink_swallow);
- scheduler_addunlock(sched, cj->hydro.super->hydro.sorts,
- t_sink_swallow);
+ t_sink_density);
+ scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
+ t_sink_density);
scheduler_addunlock(sched, cj->hydro.super->sinks.sink_in,
+ t_sink_density);
+ scheduler_addunlock(sched, t_sink_density,
+ cj->hydro.super->sinks.density_ghost);
+
+ /* Do the sink_swallow */
+ scheduler_addunlock(sched, cj->hydro.super->sinks.density_ghost,
t_sink_swallow);
scheduler_addunlock(sched, t_sink_swallow,
cj->hydro.super->sinks.sink_ghost1);
@@ -3355,6 +3377,9 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
/* The sink tasks */
if (with_sink) {
+ t_sink_density =
+ scheduler_addtask(sched, task_type_sub_self,
+ task_subtype_sink_density, flags, 0, ci, NULL);
t_sink_swallow =
scheduler_addtask(sched, task_type_sub_self,
task_subtype_sink_swallow, flags, 0, ci, NULL);
@@ -3411,6 +3436,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
#endif
}
if (with_sink) {
+ engine_addlink(e, &ci->sinks.density, t_sink_density);
engine_addlink(e, &ci->sinks.swallow, t_sink_swallow);
engine_addlink(e, &ci->sinks.do_sink_swallow, t_sink_do_sink_swallow);
engine_addlink(e, &ci->sinks.do_gas_swallow, t_sink_do_gas_swallow);
@@ -3495,14 +3521,18 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
if (with_sink) {
- /* Do the sink_swallow */
- scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
- t_sink_swallow);
+ /* Sink density */
scheduler_addunlock(sched, ci->hydro.super->sinks.drift,
- t_sink_swallow);
- scheduler_addunlock(sched, ci->hydro.super->hydro.sorts,
- t_sink_swallow);
+ t_sink_density);
+ scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
+ t_sink_density);
scheduler_addunlock(sched, ci->hydro.super->sinks.sink_in,
+ t_sink_density);
+ scheduler_addunlock(sched, t_sink_density,
+ ci->hydro.super->sinks.density_ghost);
+
+ /* Do the sink_swallow */
+ scheduler_addunlock(sched, ci->hydro.super->sinks.density_ghost,
t_sink_swallow);
scheduler_addunlock(sched, t_sink_swallow,
ci->hydro.super->sinks.sink_ghost1);
@@ -3648,6 +3678,9 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
/* The sink tasks */
if (with_sink) {
+ t_sink_density =
+ scheduler_addtask(sched, task_type_sub_pair,
+ task_subtype_sink_density, flags, 0, ci, cj);
t_sink_swallow =
scheduler_addtask(sched, task_type_sub_pair,
task_subtype_sink_swallow, flags, 0, ci, cj);
@@ -3732,15 +3765,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
#endif
}
if (with_sink) {
- /* Formation */
+ engine_addlink(e, &ci->sinks.density, t_sink_density);
+ engine_addlink(e, &cj->sinks.density, t_sink_density);
engine_addlink(e, &ci->sinks.swallow, t_sink_swallow);
engine_addlink(e, &cj->sinks.swallow, t_sink_swallow);
-
- /* Merger */
engine_addlink(e, &ci->sinks.do_sink_swallow, t_sink_do_sink_swallow);
engine_addlink(e, &cj->sinks.do_sink_swallow, t_sink_do_sink_swallow);
-
- /* Accretion */
engine_addlink(e, &ci->sinks.do_gas_swallow, t_sink_do_gas_swallow);
engine_addlink(e, &cj->sinks.do_gas_swallow, t_sink_do_gas_swallow);
}
@@ -3862,14 +3892,18 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
if (with_sink) {
- /* Do the sink_swallow */
- scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
- t_sink_swallow);
+ /* Sink density */
scheduler_addunlock(sched, ci->hydro.super->sinks.drift,
- t_sink_swallow);
- scheduler_addunlock(sched, ci->hydro.super->hydro.sorts,
- t_sink_swallow);
+ t_sink_density);
+ scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
+ t_sink_density);
scheduler_addunlock(sched, ci->hydro.super->sinks.sink_in,
+ t_sink_density);
+ scheduler_addunlock(sched, t_sink_density,
+ ci->hydro.super->sinks.density_ghost);
+
+ /* Do the sink_swallow */
+ scheduler_addunlock(sched, ci->hydro.super->sinks.density_ghost,
t_sink_swallow);
scheduler_addunlock(sched, t_sink_swallow,
ci->hydro.super->sinks.sink_ghost1);
@@ -4014,16 +4048,21 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
scheduler_addunlock(sched, t_star_feedback,
cj->hydro.super->stars.stars_out);
}
+
if (with_sink) {
- /* Do the sink_swallow */
- scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
- t_sink_swallow);
+ /* Sink density */
scheduler_addunlock(sched, cj->hydro.super->sinks.drift,
- t_sink_swallow);
- scheduler_addunlock(sched, cj->hydro.super->hydro.sorts,
- t_sink_swallow);
+ t_sink_density);
+ scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
+ t_sink_density);
scheduler_addunlock(sched, cj->hydro.super->sinks.sink_in,
+ t_sink_density);
+ scheduler_addunlock(sched, t_sink_density,
+ cj->hydro.super->sinks.density_ghost);
+
+ /* Do the sink_swallow */
+ scheduler_addunlock(sched, cj->hydro.super->sinks.density_ghost,
t_sink_swallow);
scheduler_addunlock(sched, t_sink_swallow,
cj->hydro.super->sinks.sink_ghost1);
diff --git a/src/runner.h b/src/runner.h
index 6beae93b3f5e2b0bd9c13f8e809b886ae02a4e64..73682b306dfa54e907b6a27ac0c9ec599b85d4d5 100644
--- a/src/runner.h
+++ b/src/runner.h
@@ -98,6 +98,7 @@ void runner_do_black_holes_density_ghost(struct runner *r, struct cell *c,
int timer);
void runner_do_black_holes_swallow_ghost(struct runner *r, struct cell *c,
int timer);
+void runner_do_sinks_density_ghost(struct runner *r, struct cell *c, int timer);
void runner_do_init_grav(struct runner *r, struct cell *c, int timer);
void runner_do_hydro_sort(struct runner *r, struct cell *c, int flag,
int cleanup, int rt_requests_sort, int clock);
diff --git a/src/runner_doiact_functions_sinks.h b/src/runner_doiact_functions_sinks.h
index 959d79eaa025ba29edc8cf19db7e26adc7fad59e..5128f785fddb423a8a1e8fd1736a0aae5da964ad 100644
--- a/src/runner_doiact_functions_sinks.h
+++ b/src/runner_doiact_functions_sinks.h
@@ -92,7 +92,8 @@ void DOSELF1_SINKS(struct runner *r, struct cell *c, int timer) {
if (r2 < ri2) {
IACT_SINKS_GAS(r2, dx, ri, hj, si, pj, with_cosmology, cosmo,
- e->gravity_properties, e->sink_properties);
+ e->gravity_properties, e->sink_properties,
+ e->ti_current, e->time);
}
} /* loop over the parts in ci. */
} /* loop over the sinks in ci. */
@@ -148,7 +149,8 @@ void DOSELF1_SINKS(struct runner *r, struct cell *c, int timer) {
if (r2 < ri2 || r2 < rj2) {
IACT_SINKS_SINK(r2, dx, ri, rj, si, sj, with_cosmology, cosmo,
- e->gravity_properties, e->sink_properties);
+ e->gravity_properties, e->sink_properties,
+ e->ti_current, e->time);
}
} /* loop over the sinks in ci. */
} /* loop over the sinks in ci. */
@@ -241,7 +243,8 @@ void DO_NONSYM_PAIR1_SINKS_NAIVE(struct runner *r, struct cell *restrict ci,
if (r2 < ri2) {
IACT_SINKS_GAS(r2, dx, ri, hj, si, pj, with_cosmology, cosmo,
- e->gravity_properties, e->sink_properties);
+ e->gravity_properties, e->sink_properties,
+ e->ti_current, e->time);
}
} /* loop over the parts in cj. */
} /* loop over the sinks in ci. */
@@ -297,7 +300,8 @@ void DO_NONSYM_PAIR1_SINKS_NAIVE(struct runner *r, struct cell *restrict ci,
if (r2 < ri2 || r2 < rj2) {
IACT_SINKS_SINK(r2, dx, ri, rj, si, sj, with_cosmology, cosmo,
- e->gravity_properties, e->sink_properties);
+ e->gravity_properties, e->sink_properties,
+ e->ti_current, e->time);
}
} /* loop over the sinks in cj. */
} /* loop over the sinks in ci. */
diff --git a/src/runner_doiact_sinks.c b/src/runner_doiact_sinks.c
index 5b02344c52133c58a88db70983514f5323bd1dfe..7e65564731371856eba74336b50a3c65bb58fb64 100644
--- a/src/runner_doiact_sinks.c
+++ b/src/runner_doiact_sinks.c
@@ -31,6 +31,12 @@
#include "space_getsid.h"
#include "timers.h"
+/* Import the sink density loop functions. */
+#define FUNCTION density
+#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
+#include "runner_doiact_functions_sinks.h"
+#include "runner_doiact_undef.h"
+
/* Import the sink swallow loop functions. */
#define FUNCTION swallow
#define FUNCTION_TASK_LOOP TASK_LOOP_SWALLOW
diff --git a/src/runner_ghost.c b/src/runner_ghost.c
index 8d0bcc544fe181e08a3b2ca00d1db5bd4628181f..38aaf9fd523af8a7c7b72a7dc8d561600f5c71e1 100644
--- a/src/runner_ghost.c
+++ b/src/runner_ghost.c
@@ -34,6 +34,7 @@
#include "feedback.h"
#include "mhd.h"
#include "rt.h"
+#include "sink.h"
#include "space_getsid.h"
#include "star_formation.h"
#include "stars.h"
@@ -1707,3 +1708,88 @@ void runner_do_rt_ghost2(struct runner *r, struct cell *c, int timer) {
if (timer) TIMER_TOC(timer_do_rt_ghost2);
}
+
+/**
+ * @brief Intermediate task after the density to finish density calculation
+ * and calculate accretion rates for the particle swallowing step
+ *
+ * @param r The runner thread.
+ * @param c The cell.
+ * @param timer Are we timing this ?
+ */
+void runner_do_sinks_density_ghost(struct runner *r, struct cell *c,
+ int timer) {
+
+ struct sink *restrict sinks = c->sinks.parts;
+ const struct engine *e = r->e;
+ const struct cosmology *cosmo = e->cosmology;
+ const int with_cosmology = e->policy & engine_policy_cosmology;
+
+ TIMER_TIC;
+
+ /* Anything to do here? */
+ if (c->sinks.count == 0) return;
+ if (!cell_is_active_sinks(c, e)) return;
+
+ /* Recurse? */
+ if (c->split) {
+ for (int k = 0; k < 8; k++) {
+ if (c->progeny[k] != NULL) {
+ runner_do_sinks_density_ghost(r, c->progeny[k], 0);
+ }
+ }
+ } else {
+
+ /* Init the list of active particles that have to be updated. */
+ int *sid = NULL;
+ if ((sid = (int *)malloc(sizeof(int) * c->sinks.count)) == NULL)
+ error("Can't allocate memory for sid.");
+
+ int scount = 0;
+ for (int k = 0; k < c->sinks.count; k++)
+ if (sink_is_active(&sinks[k], e)) {
+ sid[scount] = k;
+ ++scount;
+ }
+
+ /* Loop over the remaining active parts in this cell. */
+ for (int i = 0; i < scount; i++) {
+
+ /* Get a direct pointer on the part. */
+ struct sink *sp = &sinks[sid[i]];
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Is this part within the timestep? */
+ if (!sink_is_active(sp, e)) error("Ghost applied to inactive particle");
+#endif
+
+ /* Finish the density calculation */
+ sink_end_density(sp, cosmo);
+
+ if (sp->num_ngbs == 0) {
+ sinks_sink_has_no_neighbours(sp, cosmo);
+ }
+
+ /* Get particle time-step */
+ double dt;
+ if (with_cosmology) {
+ const integertime_t ti_step = get_integer_timestep(sp->time_bin);
+ const integertime_t ti_begin =
+ get_integer_time_begin(e->ti_current - 1, sp->time_bin);
+
+ dt = cosmology_get_delta_time(e->cosmology, ti_begin,
+ ti_begin + ti_step);
+ } else {
+ dt = get_timestep(sp->time_bin, e->time_base);
+ }
+
+ /* Calculate the accretion rate and accreted mass this timestep, for use
+ * in swallow loop */
+ sink_prepare_swallow(sp, e->sink_properties, e->physical_constants,
+ e->cosmology, e->cooling_func, e->entropy_floor,
+ e->time, with_cosmology, dt, e->ti_current);
+ }
+ }
+
+ if (timer) TIMER_TOC(timer_do_sinks_ghost);
+}
diff --git a/src/runner_main.c b/src/runner_main.c
index 157d6943c47cbe4836e16e63918b23c17c8b5d52..0ada1a76171405963c44ac2c218d4d7cb32041c7 100644
--- a/src/runner_main.c
+++ b/src/runner_main.c
@@ -112,6 +112,12 @@
#include "runner_doiact_black_holes.h"
#include "runner_doiact_undef.h"
+/* Import the sink density loop functions. */
+#define FUNCTION density
+#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
+#include "runner_doiact_sinks.h"
+#include "runner_doiact_undef.h"
+
/* Import the sink swallow loop functions. */
#define FUNCTION swallow
#define FUNCTION_TASK_LOOP TASK_LOOP_SWALLOW
@@ -237,6 +243,8 @@ void *runner_main(void *data) {
runner_doself1_branch_rt_gradient(r, ci);
else if (t->subtype == task_subtype_rt_transport)
runner_doself2_branch_rt_transport(r, ci);
+ else if (t->subtype == task_subtype_sink_density)
+ runner_doself_branch_sinks_density(r, ci);
else if (t->subtype == task_subtype_sink_swallow)
runner_doself_branch_sinks_swallow(r, ci);
else if (t->subtype == task_subtype_sink_do_gas_swallow)
@@ -285,6 +293,8 @@ void *runner_main(void *data) {
runner_dopair1_branch_rt_gradient(r, ci, cj);
else if (t->subtype == task_subtype_rt_transport)
runner_dopair2_branch_rt_transport(r, ci, cj);
+ else if (t->subtype == task_subtype_sink_density)
+ runner_dopair_branch_sinks_density(r, ci, cj);
else if (t->subtype == task_subtype_sink_swallow)
runner_dopair_branch_sinks_swallow(r, ci, cj);
else if (t->subtype == task_subtype_sink_do_gas_swallow)
@@ -331,6 +341,8 @@ void *runner_main(void *data) {
runner_dosub_self1_rt_gradient(r, ci, 1);
else if (t->subtype == task_subtype_rt_transport)
runner_dosub_self2_rt_transport(r, ci, 1);
+ else if (t->subtype == task_subtype_sink_density)
+ runner_dosub_self_sinks_density(r, ci, 1);
else if (t->subtype == task_subtype_sink_swallow)
runner_dosub_self_sinks_swallow(r, ci, 1);
else if (t->subtype == task_subtype_sink_do_gas_swallow)
@@ -377,6 +389,8 @@ void *runner_main(void *data) {
runner_dosub_pair1_rt_gradient(r, ci, cj, 1);
else if (t->subtype == task_subtype_rt_transport)
runner_dosub_pair2_rt_transport(r, ci, cj, 1);
+ else if (t->subtype == task_subtype_sink_density)
+ runner_dosub_pair_sinks_density(r, ci, cj, 1);
else if (t->subtype == task_subtype_sink_swallow)
runner_dosub_pair_sinks_swallow(r, ci, cj, 1);
else if (t->subtype == task_subtype_sink_do_gas_swallow)
@@ -436,6 +450,9 @@ void *runner_main(void *data) {
case task_type_bh_swallow_ghost3:
runner_do_black_holes_swallow_ghost(r, ci, 1);
break;
+ case task_type_sink_density_ghost:
+ runner_do_sinks_density_ghost(r, ci, 1);
+ break;
case task_type_drift_part:
runner_do_drift_part(r, ci, 1);
break;
diff --git a/src/scheduler.c b/src/scheduler.c
index cfec9e728b428be8f71c99d59b4e825312ebe781..358377b6914759289240d3c3c9be0c70cd3571be 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -2004,7 +2004,8 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
t->subtype == task_subtype_stars_prep2 ||
t->subtype == task_subtype_stars_feedback)
cost = 1.f * wscale * scount_i * count_i;
- else if (t->subtype == task_subtype_sink_swallow ||
+ else if (t->subtype == task_subtype_sink_density ||
+ t->subtype == task_subtype_sink_swallow ||
t->subtype == task_subtype_sink_do_gas_swallow)
cost = 1.f * wscale * count_i * sink_count_i;
else if (t->subtype == task_subtype_sink_do_sink_swallow)
@@ -2050,7 +2051,8 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
cost = 2.f * wscale * (scount_i * count_j + scount_j * count_i) *
sid_scale[t->flags];
- } else if (t->subtype == task_subtype_sink_swallow ||
+ } else if (t->subtype == task_subtype_sink_density ||
+ t->subtype == task_subtype_sink_swallow ||
t->subtype == task_subtype_sink_do_gas_swallow) {
if (t->ci->nodeID != nodeID)
cost = 3.f * wscale * count_i * sink_count_j * sid_scale[t->flags];
@@ -2126,7 +2128,8 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
sid_scale[t->flags];
}
- } else if (t->subtype == task_subtype_sink_swallow ||
+ } else if (t->subtype == task_subtype_sink_density ||
+ t->subtype == task_subtype_sink_swallow ||
t->subtype == task_subtype_sink_do_gas_swallow) {
if (t->ci->nodeID != nodeID) {
cost =
@@ -2195,7 +2198,8 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
t->subtype == task_subtype_stars_prep2 ||
t->subtype == task_subtype_stars_feedback) {
cost = 1.f * (wscale * scount_i) * count_i;
- } else if (t->subtype == task_subtype_sink_swallow ||
+ } else if (t->subtype == task_subtype_sink_density ||
+ t->subtype == task_subtype_sink_swallow ||
t->subtype == task_subtype_sink_do_gas_swallow) {
cost = 1.f * (wscale * sink_count_i) * count_i;
} else if (t->subtype == task_subtype_sink_do_sink_swallow) {
@@ -2237,6 +2241,9 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
case task_type_bh_swallow_ghost2:
if (t->ci == t->ci->hydro.super) cost = wscale * bcount_i;
break;
+ case task_type_sink_density_ghost:
+ if (t->ci == t->ci->hydro.super) cost = wscale * sink_count_i;
+ break;
case task_type_drift_part:
cost = wscale * count_i;
break;
diff --git a/src/sink/Default/sink.h b/src/sink/Default/sink.h
index 666583c5ece22897b6cbeccf65fa6eac6d14071c..f284b9d0f916b02bd7df578274953c5eb120f706 100644
--- a/src/sink/Default/sink.h
+++ b/src/sink/Default/sink.h
@@ -67,7 +67,10 @@ __attribute__((always_inline)) INLINE static void sink_init_part(
* @param sp The #sink particle to act upon.
*/
__attribute__((always_inline)) INLINE static void sink_init_sink(
- struct sink* sp) {}
+ struct sink* sp) {
+
+ sp->num_ngbs = 0;
+}
/**
* @brief Predict additional particle fields forward in time when drifting
@@ -96,6 +99,49 @@ __attribute__((always_inline)) INLINE static void sink_reset_predicted_values(
__attribute__((always_inline)) INLINE static void sink_kick_extra(
struct sink* sp, float dt) {}
+/**
+ * @brief Finishes the calculation of density on sinks
+ *
+ * @param si The particle to act upon
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void sink_end_density(
+ struct sink* si, const struct cosmology* cosmo) {}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #sink has 0
+ * ngbs.
+ *
+ * @param sp The particle to act upon
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void sinks_sink_has_no_neighbours(
+ struct sink* restrict sp, const struct cosmology* cosmo) {}
+
+/**
+ * @brief Compute the accretion rate of the sink and any quantities
+ * required swallowing based on an accretion rate
+ *
+ * Adapted from black_holes_prepare_feedback
+ *
+ * @param si The sink particle.
+ * @param props The properties of the sink scheme.
+ * @param constants The physical constants (in internal units).
+ * @param cosmo The cosmological model.
+ * @param cooling Properties of the cooling model.
+ * @param floor_props Properties of the entropy floor.
+ * @param time Time since the start of the simulation (non-cosmo mode).
+ * @param with_cosmology Are we running with cosmology?
+ * @param dt The time-step size (in physical internal units).
+ * @param ti_begin Integer time value at the beginning of timestep
+ */
+__attribute__((always_inline)) INLINE static void sink_prepare_swallow(
+ struct sink* restrict si, const struct sink_props* props,
+ const struct phys_const* constants, const struct cosmology* cosmo,
+ const struct cooling_function_data* cooling,
+ const struct entropy_floor_properties* floor_props, const double time,
+ const int with_cosmology, const double dt, const integertime_t ti_begin) {}
+
/**
* @brief Calculate if the gas has the potential of becoming
* a sink.
diff --git a/src/sink/Default/sink_iact.h b/src/sink/Default/sink_iact.h
index 72373da7292fbfd33cb26017a42786d35e79346c..bae7b6129f15c8d248911caef72f36e901e2a790 100644
--- a/src/sink/Default/sink_iact.h
+++ b/src/sink/Default/sink_iact.h
@@ -51,12 +51,47 @@ __attribute__((always_inline)) INLINE static void runner_iact_sink(
* @param a Current scale factor.
* @param H Current Hubble parameter.
* @param cut_off_radius Sink cut off radius.
+ * @param ti_current Current integer time value (for random numbers).
+ * @param time current physical time in the simulation
*/
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_sink(
const float r2, const float dx[3], const float hi, const float hj,
struct part *restrict pi, const struct part *restrict pj, const float a,
const float H, const float cut_off_radius) {}
+/**
+ * @brief Density interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param ri Comoving cut off radius of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param si First particle (sink).
+ * @param pj Second particle (gas, not updated).
+ * @param with_cosmology Are we doing a cosmological run?
+ * @param cosmo The cosmological model.
+ * @param grav_props The properties of the gravity scheme (softening, G, ...).
+ * @param sink_props the sink properties to use.
+ * @param ti_current Current integer time value (for random numbers).
+ * @param time current physical time in the simulation
+ */
+__attribute__((always_inline)) INLINE static void
+runner_iact_nonsym_sinks_gas_density(
+ const float r2, const float dx[3], const float ri, const float hj,
+ struct sink *si, const struct part *pj, const int with_cosmology,
+ const struct cosmology *cosmo, const struct gravity_props *grav_props,
+ const struct sink_props *sink_props, const integertime_t ti_current,
+ const double time) {
+
+ /* Get r. */
+ const float r = sqrtf(r2);
+
+ if (r < ri) {
+ /* Contribution to the number of neighbours in cutoff radius */
+ si->num_ngbs++;
+ }
+}
+
/**
* @brief Compute sink-sink swallow interaction (non-symmetric).
*
@@ -70,6 +105,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_sink(
* @param cosmo The cosmological parameters and properties.
* @param grav_props The gravity scheme parameters and properties.
* @param sink_props the sink properties to use.
+ * @param ti_current Current integer time value (for random numbers).
+ * @param time current physical time in the simulation
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_sinks_sink_swallow(
@@ -77,7 +114,8 @@ runner_iact_nonsym_sinks_sink_swallow(
struct sink *restrict si, struct sink *restrict sj,
const int with_cosmology, const struct cosmology *cosmo,
const struct gravity_props *grav_props,
- const struct sink_props *sink_properties) {}
+ const struct sink_props *sink_properties, const integertime_t ti_current,
+ const double time) {}
/**
* @brief Compute sink-gas swallow interaction (non-symmetric).
@@ -92,16 +130,16 @@ runner_iact_nonsym_sinks_sink_swallow(
* @param cosmo The cosmological parameters and properties.
* @param grav_props The gravity scheme parameters and properties.
* @param sink_props the sink properties to use.
+ * @param ti_current Current integer time value (for random numbers).
+ * @param time current physical time in the simulation
*/
__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_sinks_gas_swallow(const float r2, const float dx[3],
- const float ri, const float hj,
- struct sink *restrict si,
- struct part *restrict pj,
- const int with_cosmology,
- const struct cosmology *cosmo,
- const struct gravity_props *grav_props,
- const struct sink_props *sink_properties) {
-}
+runner_iact_nonsym_sinks_gas_swallow(
+ const float r2, const float dx[3], const float ri, const float hj,
+ struct sink *restrict si, struct part *restrict pj,
+ const int with_cosmology, const struct cosmology *cosmo,
+ const struct gravity_props *grav_props,
+ const struct sink_props *sink_properties, const integertime_t ti_current,
+ const double time) {}
#endif
diff --git a/src/sink/Default/sink_part.h b/src/sink/Default/sink_part.h
index 7df2a99a7d8198ee2bba807fa7bf526f9ea57d08..25a7b2e765a0316d1d5f3b78d47fb99382ca1c69 100644
--- a/src/sink/Default/sink_part.h
+++ b/src/sink/Default/sink_part.h
@@ -54,6 +54,9 @@ struct sink {
/*! Particle time bin */
timebin_t time_bin;
+ /*! Integer number of neighbours */
+ int num_ngbs;
+
#ifdef SWIFT_DEBUG_CHECKS
/* Time of the last drift */
diff --git a/src/sink/GEAR/sink.h b/src/sink/GEAR/sink.h
index 977c580ed2b8f37752becb32511e4ad23e159497..907d3d38aeee6ab55b5409324cf8e6f4d90e7f0e 100644
--- a/src/sink/GEAR/sink.h
+++ b/src/sink/GEAR/sink.h
@@ -192,6 +192,8 @@ __attribute__((always_inline)) INLINE static void sink_init_sink(
/* Reset to the mass of the sink */
sp->mass_tot_before_star_spawning = sp->mass;
+ sp->num_ngbs = 0;
+
#ifdef DEBUG_INTERACTIONS_SINKS
for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_SINKS; ++i)
sp->ids_ngbs_accretion[i] = -1;
@@ -234,6 +236,49 @@ __attribute__((always_inline)) INLINE static void sink_reset_predicted_values(
__attribute__((always_inline)) INLINE static void sink_kick_extra(
struct sink* sp, float dt) {}
+/**
+ * @brief Finishes the calculation of density on sinks
+ *
+ * @param si The particle to act upon
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void sink_end_density(
+ struct sink* si, const struct cosmology* cosmo) {}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #sink has 0
+ * ngbs.
+ *
+ * @param sp The particle to act upon
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void sinks_sink_has_no_neighbours(
+ struct sink* restrict sp, const struct cosmology* cosmo) {}
+
+/**
+ * @brief Compute the accretion rate of the sink and any quantities
+ * required swallowing based on an accretion rate
+ *
+ * Adapted from black_holes_prepare_feedback
+ *
+ * @param si The sink particle.
+ * @param props The properties of the sink scheme.
+ * @param constants The physical constants (in internal units).
+ * @param cosmo The cosmological model.
+ * @param cooling Properties of the cooling model.
+ * @param floor_props Properties of the entropy floor.
+ * @param time Time since the start of the simulation (non-cosmo mode).
+ * @param with_cosmology Are we running with cosmology?
+ * @param dt The time-step size (in physical internal units).
+ * @param ti_begin Integer time value at the beginning of timestep
+ */
+__attribute__((always_inline)) INLINE static void sink_prepare_swallow(
+ struct sink* restrict si, const struct sink_props* props,
+ const struct phys_const* constants, const struct cosmology* cosmo,
+ const struct cooling_function_data* cooling,
+ const struct entropy_floor_properties* floor_props, const double time,
+ const int with_cosmology, const double dt, const integertime_t ti_begin) {}
+
/**
* @brief Compute the rotational energy of the neighbouring gas particles.
*
diff --git a/src/sink/GEAR/sink_iact.h b/src/sink/GEAR/sink_iact.h
index 1c859cb696983f8dec35fcbdcee1d0785e0705ab..fc75ba8842bb7ceaae6d789b277a76b21206c29f 100644
--- a/src/sink/GEAR/sink_iact.h
+++ b/src/sink/GEAR/sink_iact.h
@@ -117,6 +117,39 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_sink(
}
}
+/**
+ * @brief Density interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param ri Comoving cut off radius of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param si First particle (sink).
+ * @param pj Second particle (gas, not updated).
+ * @param with_cosmology Are we doing a cosmological run?
+ * @param cosmo The cosmological model.
+ * @param grav_props The properties of the gravity scheme (softening, G, ...).
+ * @param sink_props the sink properties to use.
+ * @param ti_current Current integer time value (for random numbers).
+ * @param time current physical time in the simulation
+ */
+__attribute__((always_inline)) INLINE static void
+runner_iact_nonsym_sinks_gas_density(
+ const float r2, const float dx[3], const float ri, const float hj,
+ struct sink *si, const struct part *pj, const int with_cosmology,
+ const struct cosmology *cosmo, const struct gravity_props *grav_props,
+ const struct sink_props *sink_props, const integertime_t ti_current,
+ const double time) {
+
+ /* Get r. */
+ const float r = sqrtf(r2);
+
+ if (r < ri) {
+ /* Contribution to the number of neighbours in cutoff radius */
+ si->num_ngbs++;
+ }
+}
+
/**
* @brief Compute sink-sink swallow interaction (non-symmetric).
*
@@ -132,6 +165,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_sink(
* @param cosmo The cosmological parameters and properties.
* @param grav_props The gravity scheme parameters and properties.
* @param sink_props the sink properties to use.
+ * @param ti_current Current integer time value (for random numbers).
+ * @param time current physical time in the simulation
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_sinks_sink_swallow(
@@ -139,7 +174,8 @@ runner_iact_nonsym_sinks_sink_swallow(
struct sink *restrict si, struct sink *restrict sj,
const int with_cosmology, const struct cosmology *cosmo,
const struct gravity_props *grav_props,
- const struct sink_props *sink_properties) {
+ const struct sink_props *sink_properties, const integertime_t ti_current,
+ const double time) {
const float r = sqrtf(r2);
const float f_acc_r_acc_i = sink_properties->f_acc * ri;
@@ -260,16 +296,17 @@ runner_iact_nonsym_sinks_sink_swallow(
* @param cosmo The cosmological parameters and properties.
* @param grav_props The gravity scheme parameters and properties.
* @param sink_props the sink properties to use.
+ * @param ti_current Current integer time value (for random numbers).
+ * @param time current physical time in the simulation
*/
__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_sinks_gas_swallow(const float r2, const float dx[3],
- const float ri, const float hj,
- struct sink *restrict si,
- struct part *restrict pj,
- const int with_cosmology,
- const struct cosmology *cosmo,
- const struct gravity_props *grav_props,
- const struct sink_props *sink_properties) {
+runner_iact_nonsym_sinks_gas_swallow(
+ const float r2, const float dx[3], const float ri, const float hj,
+ struct sink *restrict si, struct part *restrict pj,
+ const int with_cosmology, const struct cosmology *cosmo,
+ const struct gravity_props *grav_props,
+ const struct sink_props *sink_properties, const integertime_t ti_current,
+ const double time) {
const float r = sqrtf(r2);
const float f_acc_r_acc = sink_properties->f_acc * ri;
diff --git a/src/sink/GEAR/sink_part.h b/src/sink/GEAR/sink_part.h
index 0d0ac23ac5b376c6f1cf1e9de12e73a45092231c..be7f4359972d38570d944f1be8cae187fb6d83f1 100644
--- a/src/sink/GEAR/sink_part.h
+++ b/src/sink/GEAR/sink_part.h
@@ -55,6 +55,9 @@ struct sink {
/*! Sink target mass. In Msun. */
float target_mass_Msun;
+ /*! Integer number of neighbours */
+ int num_ngbs;
+
/*! Mass of the sink before starting the star spawning loop */
float mass_tot_before_star_spawning;
diff --git a/src/space_recycle.c b/src/space_recycle.c
index cf842273027b22e0123f6efa7a937797e9797f26..1a7a42c830e3636815feea2b4ad5e1d70bdccb45 100644
--- a/src/space_recycle.c
+++ b/src/space_recycle.c
@@ -138,6 +138,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->stars.feedback = NULL;
c->stars.prepare1 = NULL;
c->stars.prepare2 = NULL;
+ c->sinks.density = NULL;
c->sinks.swallow = NULL;
c->sinks.do_sink_swallow = NULL;
c->sinks.do_gas_swallow = NULL;
@@ -169,6 +170,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->black_holes.black_holes_in = NULL;
c->black_holes.black_holes_out = NULL;
c->sinks.sink_in = NULL;
+ c->sinks.density_ghost = NULL;
c->sinks.sink_ghost1 = NULL;
c->sinks.sink_ghost2 = NULL;
c->sinks.sink_out = NULL;
diff --git a/src/task.c b/src/task.c
index 2b1147cc4b01930dcd1c89f187dc642ca672acd4..9a4c7f2cb2311bcd5e88e9f67b07078b4f52b067 100644
--- a/src/task.c
+++ b/src/task.c
@@ -114,6 +114,7 @@ const char *taskID_names[task_type_count] = {
"fof_attach_pair",
"neutrino_weight",
"sink_in",
+ "sink_density_ghost",
"sink_ghost1",
"sink_ghost2",
"sink_out",
@@ -160,6 +161,7 @@ const char *subtaskID_names[task_subtype_count] = {
"do_gas_swallow",
"do_bh_swallow",
"bh_feedback",
+ "sink_density",
"sink_do_sink_swallow",
"sink_swallow",
"sink_do_gas_swallow",
@@ -248,6 +250,7 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
break;
case task_type_drift_sink:
+ case task_type_sink_density_ghost:
return task_action_sink;
break;
@@ -293,6 +296,7 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
return task_action_bpart;
break;
+ case task_subtype_sink_density:
case task_subtype_sink_do_gas_swallow:
case task_subtype_sink_do_sink_swallow:
case task_subtype_sink_swallow:
@@ -575,7 +579,8 @@ void task_unlock(struct task *t) {
cell_gunlocktree(ci);
cell_munlocktree(ci);
#endif
- } else if ((subtype == task_subtype_sink_swallow) ||
+ } else if ((subtype == task_subtype_sink_density) ||
+ (subtype == task_subtype_sink_swallow) ||
(subtype == task_subtype_sink_do_gas_swallow)) {
cell_sink_unlocktree(ci);
cell_unlocktree(ci);
@@ -613,7 +618,8 @@ void task_unlock(struct task *t) {
cell_munlocktree(ci);
cell_munlocktree(cj);
#endif
- } else if ((subtype == task_subtype_sink_swallow) ||
+ } else if ((subtype == task_subtype_sink_density) ||
+ (subtype == task_subtype_sink_swallow) ||
(subtype == task_subtype_sink_do_gas_swallow)) {
cell_sink_unlocktree(ci);
cell_sink_unlocktree(cj);
@@ -807,7 +813,8 @@ int task_lock(struct task *t) {
return 0;
}
#endif
- } else if ((subtype == task_subtype_sink_swallow) ||
+ } else if ((subtype == task_subtype_sink_density) ||
+ (subtype == task_subtype_sink_swallow) ||
(subtype == task_subtype_sink_do_gas_swallow)) {
if (ci->sinks.hold) return 0;
if (ci->hydro.hold) return 0;
@@ -876,7 +883,8 @@ int task_lock(struct task *t) {
return 0;
}
#endif
- } else if ((subtype == task_subtype_sink_swallow) ||
+ } else if ((subtype == task_subtype_sink_density) ||
+ (subtype == task_subtype_sink_swallow) ||
(subtype == task_subtype_sink_do_gas_swallow)) {
if (ci->sinks.hold || cj->sinks.hold) return 0;
if (ci->hydro.hold || cj->hydro.hold) return 0;
@@ -1192,6 +1200,9 @@ void task_get_group_name(int type, int subtype, char *cluster) {
strcpy(cluster, "RTtransport");
}
break;
+ case task_subtype_sink_density:
+ strcpy(cluster, "SinkDensity");
+ break;
case task_subtype_sink_swallow:
strcpy(cluster, "SinkSwallow");
break;
@@ -1669,6 +1680,7 @@ enum task_categories task_get_category(const struct task *t) {
case task_type_star_formation_sink:
return task_category_star_formation;
+ case task_type_sink_density_ghost:
case task_type_sink_formation:
return task_category_sink;
@@ -1778,6 +1790,7 @@ enum task_categories task_get_category(const struct task *t) {
case task_subtype_bh_feedback:
return task_category_black_holes;
+ case task_subtype_sink_density:
case task_subtype_sink_swallow:
case task_subtype_sink_do_sink_swallow:
case task_subtype_sink_do_gas_swallow:
diff --git a/src/task.h b/src/task.h
index b718326384d12bbcba69e5d8f901193f7ffa565c..1e31478890accf1b3975ec692d43deec26fe99db 100644
--- a/src/task.h
+++ b/src/task.h
@@ -106,7 +106,8 @@ enum task_types {
task_type_fof_attach_self,
task_type_fof_attach_pair,
task_type_neutrino_weight,
- task_type_sink_in, /* Implicit */
+ task_type_sink_in, /* Implicit */
+ task_type_sink_density_ghost,
task_type_sink_ghost1, /* Implicit */
task_type_sink_ghost2, /* Implicit */
task_type_sink_out, /* Implicit */
@@ -156,6 +157,7 @@ enum task_subtypes {
task_subtype_do_gas_swallow,
task_subtype_do_bh_swallow,
task_subtype_bh_feedback,
+ task_subtype_sink_density,
task_subtype_sink_do_sink_swallow,
task_subtype_sink_swallow,
task_subtype_sink_do_gas_swallow,
diff --git a/src/timers.c b/src/timers.c
index e21b99c30fa97d644df15446d2a3de567ec0fa8c..f23ffbc1b9f3dc1fc3dc98c9f0a2f5522331b5c5 100644
--- a/src/timers.c
+++ b/src/timers.c
@@ -61,6 +61,7 @@ const char* timers_names[timer_count] = {
"doself_bh_swallow",
"doself_bh_feedback",
"doself_grav_pp",
+ "doself_sink_density",
"doself_sink_swallow",
"dopair_density",
"dopair_gradient",
@@ -73,6 +74,7 @@ const char* timers_names[timer_count] = {
"dopair_bh_feedback",
"dopair_grav_mm",
"dopair_grav_pp",
+ "dopair_sink_density",
"dopair_sink_swallow",
"dograv_external",
"dograv_down",
@@ -89,6 +91,7 @@ const char* timers_names[timer_count] = {
"dosub_self_bh_swallow",
"dosub_self_bh_feedback",
"dosub_self_grav",
+ "dosub_self_sink_density",
"dosub_self_sink_swallow",
"dosub_pair_density",
"dosub_pair_gradient",
@@ -100,6 +103,7 @@ const char* timers_names[timer_count] = {
"dosub_pair_bh_swallow",
"dosub_pair_bh_feedback",
"dosub_pair_grav",
+ "dosub_pair_sink_density",
"dosub_pair_sink_swallow",
"doself_subset",
"dopair_subset",
@@ -109,6 +113,7 @@ const char* timers_names[timer_count] = {
"do_extra_ghost",
"do_stars_ghost",
"do_black_holes_ghost",
+ "do_sinks_ghost",
"dorecv_part",
"dorecv_gpart",
"dorecv_spart",
diff --git a/src/timers.h b/src/timers.h
index f5546a3ce1a14a3e192bc2c3b7ee2363fc9562db..fef6a6ebf8b81bb776d07e89b8ffedd78caefdcf 100644
--- a/src/timers.h
+++ b/src/timers.h
@@ -61,6 +61,7 @@ enum {
timer_doself_bh_swallow,
timer_doself_bh_feedback,
timer_doself_grav_pp,
+ timer_doself_sink_density,
timer_doself_sink_swallow,
timer_dopair_density,
timer_dopair_gradient,
@@ -73,6 +74,7 @@ enum {
timer_dopair_bh_feedback,
timer_dopair_grav_mm,
timer_dopair_grav_pp,
+ timer_dopair_sink_density,
timer_dopair_sink_swallow,
timer_dograv_external,
timer_dograv_down,
@@ -89,6 +91,7 @@ enum {
timer_dosub_self_bh_swallow,
timer_dosub_self_bh_feedback,
timer_dosub_self_grav,
+ timer_dosub_self_sink_density,
timer_dosub_self_sink_swallow,
timer_dosub_pair_density,
timer_dosub_pair_gradient,
@@ -100,6 +103,7 @@ enum {
timer_dosub_pair_bh_swallow,
timer_dosub_pair_bh_feedback,
timer_dosub_pair_grav,
+ timer_dosub_pair_sink_density,
timer_dosub_pair_sink_swallow,
timer_doself_subset,
timer_dopair_subset,
@@ -109,6 +113,7 @@ enum {
timer_do_extra_ghost,
timer_do_stars_ghost,
timer_do_black_holes_ghost,
+ timer_do_sinks_ghost,
timer_dorecv_part,
timer_dorecv_gpart,
timer_dorecv_spart,
diff --git a/tools/plot_task_dependencies.py b/tools/plot_task_dependencies.py
index fade3d5de847cb33b4796b650b0e6f1aba4af79f..17013ca2f131ddf648d81f4fa70158df65bbd610 100755
--- a/tools/plot_task_dependencies.py
+++ b/tools/plot_task_dependencies.py
@@ -171,7 +171,7 @@ def task_is_black_holes(name):
name: str
Task name
"""
- if "bh" in name or "bpart" in name or "swallow" in name:
+ if "bh" in name or "bpart" in name or ("swallow" in name and not "sink" in name):
return True
return False
@@ -208,7 +208,7 @@ def task_is_hydro(name):
return True
if "_part" in name:
return True
- if "density" in name and "stars" not in name and "bh" not in name:
+ if "density" in name and "stars" not in name and "sink" not in name and "bh" not in name:
return True
if "rho" in name and "bpart" not in name:
return True