From 1568b371cfe697006936d12f59754d6be184a08c Mon Sep 17 00:00:00 2001
From: Jonathan Davies <j.j.davies@ljmu.ac.uk>
Date: Wed, 20 Nov 2024 11:42:29 +0000
Subject: [PATCH] Sink density tasks
---
AUTHORS | 3 +-
.../NewOption/sink_adding_new_scheme.rst | 13 ++
src/cell.c | 1 +
src/cell_sinks.h | 6 +
src/cell_unskip.c | 46 +++++--
src/engine_maketasks.c | 127 ++++++++++++------
src/runner.h | 1 +
src/runner_doiact_functions_sinks.h | 12 +-
src/runner_doiact_sinks.c | 6 +
src/runner_ghost.c | 86 ++++++++++++
src/runner_main.c | 17 +++
src/scheduler.c | 15 ++-
src/sink/Default/sink.h | 48 ++++++-
src/sink/Default/sink_iact.h | 58 ++++++--
src/sink/Default/sink_part.h | 3 +
src/sink/GEAR/sink.h | 45 +++++++
src/sink/GEAR/sink_iact.h | 55 ++++++--
src/sink/GEAR/sink_part.h | 3 +
src/space_recycle.c | 2 +
src/task.c | 21 ++-
src/task.h | 4 +-
src/timers.c | 5 +
src/timers.h | 5 +
tools/plot_task_dependencies.py | 4 +-
24 files changed, 493 insertions(+), 93 deletions(-)
diff --git a/AUTHORS b/AUTHORS
index ff77660e38..21268f93b6 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 b25829a18a..023446e1fd 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 a610135c87..4f58054830 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 ac1d9d5920..6bd2f9a9af 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 75c51f3f13..3ce7d791a7 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 227cdfec99..5592e384ec 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 6beae93b3f..73682b306d 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 959d79eaa0..5128f785fd 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 5b02344c52..7e65564731 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 8d0bcc544f..38aaf9fd52 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 157d6943c4..0ada1a7617 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 cfec9e728b..358377b691 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 666583c5ec..f284b9d0f9 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 72373da729..bae7b6129f 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 7df2a99a7d..25a7b2e765 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 977c580ed2..907d3d38ae 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 1c859cb696..fc75ba8842 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 0d0ac23ac5..be7f435997 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 cf84227302..1a7a42c830 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 2b1147cc4b..9a4c7f2cb2 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 b718326384..1e31478890 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 e21b99c30f..f23ffbc1b9 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 f5546a3ce1..fef6a6ebf8 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 fade3d5de8..17013ca2f1 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
--
GitLab