diff --git a/examples/main.c b/examples/main.c
index 8ff44b2dd8d05109d516d27e4f9ff8074ce8ec8e..133b7a2e261ba633001144b748eb6769f35a0671 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -519,13 +519,12 @@ int main(int argc, char *argv[]) {
#ifdef WITH_MPI
if (with_mpole_reconstruction && nr_nodes > 1)
error("Cannot reconstruct m-poles every step over MPI (yet).");
- if (with_star_formation && with_feedback)
- error("Can't run with star formation and feedback over MPI (yet)");
if (with_limiter) error("Can't run with time-step limiter over MPI (yet)");
#endif
/* Temporary early aborts for modes not supported with hand-vec. */
-#if defined(WITH_VECTORIZATION) && !defined(CHEMISTRY_NONE)
+#if defined(WITH_VECTORIZATION) && defined(GADGET2_SPH) && \
+ !defined(CHEMISTRY_NONE)
error(
"Cannot run with chemistry and hand-vectorization (yet). "
"Use --disable-hand-vec at configure time.");
diff --git a/src/active.h b/src/active.h
index e7cd35fa570f15694ab8fa3011d6a02b65837886..f15a98fceeddc741b0f66009ff4a0c3b416b32f1 100644
--- a/src/active.h
+++ b/src/active.h
@@ -496,6 +496,28 @@ __attribute__((always_inline)) INLINE static int cell_is_starting_stars(
return (c->stars.ti_beg_max == e->ti_current);
}
+/**
+ * @brief Does a cell contain any b-particle starting their time-step now ?
+ *
+ * @param c The #cell.
+ * @param e The #engine containing information about the current time.
+ * @return 1 if the #cell contains at least an active particle, 0 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int cell_is_starting_black_holes(
+ const struct cell *c, const struct engine *e) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->black_holes.ti_beg_max > e->ti_current)
+ error(
+ "cell in an impossible time-zone! c->ti_beg_max=%lld (t=%e) and "
+ "e->ti_current=%lld (t=%e, a=%e)",
+ c->black_holes.ti_beg_max, c->black_holes.ti_beg_max * e->time_base,
+ e->ti_current, e->ti_current * e->time_base, e->cosmology->a);
+#endif
+
+ return (c->black_holes.ti_beg_max == e->ti_current);
+}
+
/**
* @brief Is this particle starting its time-step now ?
*
diff --git a/src/black_holes/EAGLE/black_holes.h b/src/black_holes/EAGLE/black_holes.h
index 909057bb962cd9a21af9a708d8796f5cd18d6a8d..94962b7deea4900598a565dfd6ea90cfc6e17798 100644
--- a/src/black_holes/EAGLE/black_holes.h
+++ b/src/black_holes/EAGLE/black_holes.h
@@ -233,7 +233,7 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
4. * M_PI * G * BH_mass * proton_mass / (epsilon_r * c * sigma_Thomson);
/* Limit the accretion rate to the Eddington fraction */
- const double accr_rate = max(Bondi_rate, f_Edd * Eddington_rate);
+ const double accr_rate = min(Bondi_rate, f_Edd * Eddington_rate);
bp->accretion_rate = accr_rate;
/* Factor in the radiative efficiency */
diff --git a/src/cell.c b/src/cell.c
index 4c6bc5b6fdb168c315f777600ea47ae19d021bb1..86a0bee9ee65499f7aa4fadd471e8f7a87ac0162 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -260,6 +260,7 @@ int cell_link_sparts(struct cell *c, struct spart *sparts) {
#endif
c->stars.parts = sparts;
+ c->stars.parts_rebuild = sparts;
/* Fill the progeny recursively, depth-first. */
if (c->split) {
@@ -1054,6 +1055,92 @@ int cell_unpack_multipoles(struct cell *restrict c,
#endif
}
+/**
+ * @brief Pack the counts for star formation of the given cell and all it's
+ * sub-cells.
+ *
+ * @param c The #cell.
+ * @param pcells (output) The multipole information we pack into
+ *
+ * @return The number of packed cells.
+ */
+int cell_pack_sf_counts(struct cell *restrict c,
+ struct pcell_sf *restrict pcells) {
+
+#ifdef WITH_MPI
+
+ /* Pack this cell's data. */
+ pcells[0].stars.delta_from_rebuild = c->stars.parts - c->stars.parts_rebuild;
+ pcells[0].stars.count = c->stars.count;
+ pcells[0].stars.dx_max_part = c->stars.dx_max_part;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->stars.parts_rebuild == NULL)
+ error("Star particles array at rebuild is NULL! c->depth=%d", c->depth);
+
+ if (pcells[0].stars.delta_from_rebuild < 0)
+ error("Stars part pointer moved in the wrong direction!");
+
+ if (pcells[0].stars.delta_from_rebuild > 0 && c->depth == 0)
+ error("Shifting the top-level pointer is not allowed!");
+#endif
+
+ /* Fill in the progeny, depth-first recursion. */
+ int count = 1;
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) {
+ count += cell_pack_sf_counts(c->progeny[k], &pcells[count]);
+ }
+
+ /* Return the number of packed values. */
+ return count;
+
+#else
+ error("SWIFT was not compiled with MPI support.");
+ return 0;
+#endif
+}
+
+/**
+ * @brief Unpack the counts for star formation of a given cell and its
+ * sub-cells.
+ *
+ * @param c The #cell
+ * @param pcells The multipole information to unpack
+ *
+ * @return The number of cells created.
+ */
+int cell_unpack_sf_counts(struct cell *restrict c,
+ struct pcell_sf *restrict pcells) {
+
+#ifdef WITH_MPI
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->stars.parts_rebuild == NULL)
+ error("Star particles array at rebuild is NULL!");
+#endif
+
+ /* Unpack this cell's data. */
+ c->stars.count = pcells[0].stars.count;
+ c->stars.parts = c->stars.parts_rebuild + pcells[0].stars.delta_from_rebuild;
+ c->stars.dx_max_part = pcells[0].stars.dx_max_part;
+
+ /* Fill in the progeny, depth-first recursion. */
+ int count = 1;
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) {
+ count += cell_unpack_sf_counts(c->progeny[k], &pcells[count]);
+ }
+
+ /* Return the number of packed values. */
+ return count;
+
+#else
+ error("SWIFT was not compiled with MPI support.");
+ return 0;
+#endif
+}
+
/**
* @brief Lock a cell for access to its array of #part and hold its parents.
*
@@ -1597,6 +1684,7 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
c->progeny[k]->stars.count = bucket_count[k];
c->progeny[k]->stars.count_total = c->progeny[k]->stars.count;
c->progeny[k]->stars.parts = &c->stars.parts[bucket_offset[k]];
+ c->progeny[k]->stars.parts_rebuild = c->progeny[k]->stars.parts;
}
/* Now do the same song and dance for the bparts. */
@@ -2455,7 +2543,9 @@ void cell_activate_hydro_sorts(struct cell *c, int sid, struct scheduler *s) {
* @brief Activate the sorts up a cell hierarchy.
*/
void cell_activate_stars_sorts_up(struct cell *c, struct scheduler *s) {
+
if (c == c->hydro.super) {
+
#ifdef SWIFT_DEBUG_CHECKS
if (c->stars.sorts == NULL)
error("Trying to activate un-existing c->stars.sorts");
@@ -2465,11 +2555,17 @@ void cell_activate_stars_sorts_up(struct cell *c, struct scheduler *s) {
cell_activate_drift_spart(c, s);
}
} else {
+
+ /* Climb up the tree and set the flags */
for (struct cell *parent = c->parent;
parent != NULL && !cell_get_flag(parent, cell_flag_do_stars_sub_sort);
parent = parent->parent) {
+
cell_set_flag(parent, cell_flag_do_stars_sub_sort);
+
+ /* Reached the super-level? Activate the task and abort */
if (parent == c->hydro.super) {
+
#ifdef SWIFT_DEBUG_CHECKS
if (parent->stars.sorts == NULL)
error("Trying to activate un-existing parents->stars.sorts");
@@ -2486,8 +2582,10 @@ void cell_activate_stars_sorts_up(struct cell *c, struct scheduler *s) {
* @brief Activate the sorts on a given cell, if needed.
*/
void cell_activate_stars_sorts(struct cell *c, int sid, struct scheduler *s) {
+
/* Do we need to re-sort? */
if (c->stars.dx_max_sort > space_maxreldx * c->dmin) {
+
/* Climb up the tree to active the sorts in that direction */
for (struct cell *finger = c; finger != NULL; finger = finger->parent) {
if (finger->stars.requires_sorts) {
@@ -2559,12 +2657,12 @@ void cell_activate_subcell_hydro_tasks(struct cell *ci, struct cell *cj,
/* Get the orientation of the pair. */
double shift[3];
- int sid = space_getsid(s->space, &ci, &cj, shift);
+ const int sid = space_getsid(s->space, &ci, &cj, shift);
/* recurse? */
if (cell_can_recurse_in_pair_hydro_task(ci) &&
cell_can_recurse_in_pair_hydro_task(cj)) {
- struct cell_split_pair *csp = &cell_split_pairs[sid];
+ const struct cell_split_pair *csp = &cell_split_pairs[sid];
for (int k = 0; k < csp->count; k++) {
const int pid = csp->pairs[k].pid;
const int pjd = csp->pairs[k].pjd;
@@ -2606,9 +2704,11 @@ void cell_activate_subcell_hydro_tasks(struct cell *ci, struct cell *cj,
* @param ci The first #cell we recurse in.
* @param cj The second #cell we recurse in.
* @param s The task #scheduler.
+ * @param with_star_formation Are we running with star formation switched on?
*/
void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
- struct scheduler *s) {
+ struct scheduler *s,
+ const int with_star_formation) {
const struct engine *e = s->space->e;
/* Store the current dx_max and h_max values. */
@@ -2626,9 +2726,13 @@ void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
/* Self interaction? */
if (cj == NULL) {
+
+ const int ci_active = cell_is_active_stars(ci, e) ||
+ (with_star_formation && cell_is_active_hydro(ci, e));
+
/* Do anything? */
- if (!cell_is_active_stars(ci, e) || ci->hydro.count == 0 ||
- ci->stars.count == 0)
+ if (!ci_active || ci->hydro.count == 0 ||
+ (!with_star_formation && ci->stars.count == 0))
return;
/* Recurse? */
@@ -2636,11 +2740,12 @@ void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
/* Loop over all progenies and pairs of progenies */
for (int j = 0; j < 8; j++) {
if (ci->progeny[j] != NULL) {
- cell_activate_subcell_stars_tasks(ci->progeny[j], NULL, s);
+ cell_activate_subcell_stars_tasks(ci->progeny[j], NULL, s,
+ with_star_formation);
for (int k = j + 1; k < 8; k++)
if (ci->progeny[k] != NULL)
cell_activate_subcell_stars_tasks(ci->progeny[j], ci->progeny[k],
- s);
+ s, with_star_formation);
}
}
} else {
@@ -2652,29 +2757,38 @@ void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
/* Otherwise, pair interation */
else {
- /* Should we even bother? */
- if (!cell_is_active_stars(ci, e) && !cell_is_active_stars(cj, e)) return;
/* Get the orientation of the pair. */
double shift[3];
- int sid = space_getsid(s->space, &ci, &cj, shift);
+ const int sid = space_getsid(s->space, &ci, &cj, shift);
+
+ const int ci_active = cell_is_active_stars(ci, e) ||
+ (with_star_formation && cell_is_active_hydro(ci, e));
+ const int cj_active = cell_is_active_stars(cj, e) ||
+ (with_star_formation && cell_is_active_hydro(cj, e));
+
+ /* Should we even bother? */
+ if (!ci_active && !cj_active) return;
/* recurse? */
if (cell_can_recurse_in_pair_stars_task(ci, cj) &&
cell_can_recurse_in_pair_stars_task(cj, ci)) {
- struct cell_split_pair *csp = &cell_split_pairs[sid];
+
+ const struct cell_split_pair *csp = &cell_split_pairs[sid];
for (int k = 0; k < csp->count; k++) {
const int pid = csp->pairs[k].pid;
const int pjd = csp->pairs[k].pjd;
if (ci->progeny[pid] != NULL && cj->progeny[pjd] != NULL)
cell_activate_subcell_stars_tasks(ci->progeny[pid], cj->progeny[pjd],
- s);
+ s, with_star_formation);
}
}
/* Otherwise, activate the sorts and drifts. */
else {
- if (cell_is_active_stars(ci, e)) {
+
+ if (ci_active) {
+
/* We are going to interact this pair, so store some values. */
atomic_or(&cj->hydro.requires_sorts, 1 << sid);
atomic_or(&ci->stars.requires_sorts, 1 << sid);
@@ -2691,7 +2805,8 @@ void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
cell_activate_stars_sorts(ci, sid, s);
}
- if (cell_is_active_stars(cj, e)) {
+ if (cj_active) {
+
/* We are going to interact this pair, so store some values. */
atomic_or(&cj->stars.requires_sorts, 1 << sid);
atomic_or(&ci->hydro.requires_sorts, 1 << sid);
@@ -2771,12 +2886,12 @@ void cell_activate_subcell_black_holes_tasks(struct cell *ci, struct cell *cj,
/* Get the orientation of the pair. */
double shift[3];
- int sid = space_getsid(s->space, &ci, &cj, shift);
+ const int sid = space_getsid(s->space, &ci, &cj, shift);
/* recurse? */
if (cell_can_recurse_in_pair_black_holes_task(ci, cj) &&
cell_can_recurse_in_pair_black_holes_task(cj, ci)) {
- struct cell_split_pair *csp = &cell_split_pairs[sid];
+ const struct cell_split_pair *csp = &cell_split_pairs[sid];
for (int k = 0; k < csp->count; k++) {
const int pid = csp->pairs[k].pid;
const int pjd = csp->pairs[k].pjd;
@@ -2964,6 +3079,11 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
struct engine *e = s->space->e;
const int nodeID = e->nodeID;
const int with_limiter = (e->policy & engine_policy_limiter);
+
+#ifdef WITH_MPI
+ const int with_star_formation = e->policy & engine_policy_star_formation;
+ const int with_feedback = e->policy & engine_policy_feedback;
+#endif
int rebuild = 0;
/* Un-skip the density tasks involved with this cell. */
@@ -3012,8 +3132,14 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
cell_activate_hydro_sorts(ci, t->flags, s);
cell_activate_hydro_sorts(cj, t->flags, s);
}
+
/* Store current values of dx_max and h_max. */
- else if (t->type == task_type_sub_pair || t->type == task_type_sub_self) {
+ else if (t->type == task_type_sub_self) {
+ cell_activate_subcell_hydro_tasks(ci, NULL, s);
+ }
+
+ /* Store current values of dx_max and h_max. */
+ else if (t->type == task_type_sub_pair) {
cell_activate_subcell_hydro_tasks(ci, cj, s);
}
}
@@ -3076,6 +3202,20 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
scheduler_activate_send(s, cj->mpi.send, task_subtype_tend_part,
ci_nodeID);
+ /* Propagating new star counts? */
+ if (with_star_formation && with_feedback) {
+ if (ci_active && ci->hydro.count > 0) {
+ scheduler_activate_recv(s, ci->mpi.recv, task_subtype_sf_counts);
+ scheduler_activate_recv(s, ci->mpi.recv, task_subtype_tend_spart);
+ }
+ if (cj_active && cj->hydro.count > 0) {
+ scheduler_activate_send(s, cj->mpi.send, task_subtype_sf_counts,
+ ci_nodeID);
+ scheduler_activate_send(s, cj->mpi.send, task_subtype_tend_spart,
+ ci_nodeID);
+ }
+ }
+
} else if (cj_nodeID != nodeID) {
/* If the local cell is active, receive data from the foreign cell. */
if (ci_active) {
@@ -3126,6 +3266,20 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
if (ci_active || with_limiter)
scheduler_activate_send(s, ci->mpi.send, task_subtype_tend_part,
cj_nodeID);
+
+ /* Propagating new star counts? */
+ if (with_star_formation && with_feedback) {
+ if (cj_active && cj->hydro.count > 0) {
+ scheduler_activate_recv(s, cj->mpi.recv, task_subtype_sf_counts);
+ scheduler_activate_recv(s, cj->mpi.recv, task_subtype_tend_spart);
+ }
+ if (ci_active && ci->hydro.count > 0) {
+ scheduler_activate_send(s, ci->mpi.send, task_subtype_sf_counts,
+ cj_nodeID);
+ scheduler_activate_send(s, ci->mpi.send, task_subtype_tend_spart,
+ cj_nodeID);
+ }
+ }
}
#endif
}
@@ -3325,16 +3479,23 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
*
* @param c the #cell.
* @param s the #scheduler.
+ * @param with_star_formation Are we running with star formation switched on?
*
* @return 1 If the space needs rebuilding. 0 otherwise.
*/
-int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s) {
+int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s,
+ const int with_star_formation) {
+
struct engine *e = s->space->e;
const int nodeID = e->nodeID;
int rebuild = 0;
- if (c->stars.drift != NULL && cell_is_active_stars(c, e)) {
- cell_activate_drift_spart(c, s);
+ if (c->stars.drift != NULL) {
+ if (cell_is_active_stars(c, e) ||
+ (with_star_formation && cell_is_active_hydro(c, e))) {
+
+ cell_activate_drift_spart(c, s);
+ }
}
/* Un-skip the density tasks involved with this cell. */
@@ -3342,8 +3503,6 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s) {
struct task *t = l->t;
struct cell *ci = t->ci;
struct cell *cj = t->cj;
- const int ci_active = cell_is_active_stars(ci, e);
- const int cj_active = (cj != NULL) ? cell_is_active_stars(cj, e) : 0;
#ifdef WITH_MPI
const int ci_nodeID = ci->nodeID;
const int cj_nodeID = (cj != NULL) ? cj->nodeID : -1;
@@ -3352,6 +3511,13 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s) {
const int cj_nodeID = nodeID;
#endif
+ const int ci_active = cell_is_active_stars(ci, e) ||
+ (with_star_formation && cell_is_active_hydro(ci, e));
+
+ const int cj_active =
+ (cj != NULL) && (cell_is_active_stars(cj, e) ||
+ (with_star_formation && cell_is_active_hydro(cj, e)));
+
/* Activate the drifts */
if (t->type == task_type_self && ci_active) {
cell_activate_drift_part(ci, s);
@@ -3403,8 +3569,12 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s) {
}
}
- else if (t->type == task_type_sub_pair || t->type == task_type_sub_self) {
- cell_activate_subcell_stars_tasks(ci, cj, s);
+ else if (t->type == task_type_sub_self) {
+ cell_activate_subcell_stars_tasks(ci, NULL, s, with_star_formation);
+ }
+
+ else if (t->type == task_type_sub_pair) {
+ cell_activate_subcell_stars_tasks(ci, cj, s, with_star_formation);
}
}
@@ -3487,8 +3657,6 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s) {
struct task *t = l->t;
struct cell *ci = t->ci;
struct cell *cj = t->cj;
- const int ci_active = cell_is_active_stars(ci, e);
- const int cj_active = (cj != NULL) ? cell_is_active_stars(cj, e) : 0;
#ifdef WITH_MPI
const int ci_nodeID = ci->nodeID;
const int cj_nodeID = (cj != NULL) ? cj->nodeID : -1;
@@ -3497,6 +3665,13 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s) {
const int cj_nodeID = nodeID;
#endif
+ const int ci_active = cell_is_active_stars(ci, e) ||
+ (with_star_formation && cell_is_active_hydro(ci, e));
+
+ const int cj_active =
+ (cj != NULL) && (cell_is_active_stars(cj, e) ||
+ (with_star_formation && cell_is_active_hydro(cj, e)));
+
if (t->type == task_type_self && ci_active) {
scheduler_activate(s, t);
}
@@ -3524,14 +3699,18 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s) {
}
/* Unskip all the other task types. */
- if (c->nodeID == nodeID && cell_is_active_stars(c, e)) {
- if (c->stars.ghost != NULL) scheduler_activate(s, c->stars.ghost);
- if (c->stars.stars_in != NULL) scheduler_activate(s, c->stars.stars_in);
- if (c->stars.stars_out != NULL) scheduler_activate(s, c->stars.stars_out);
- if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
- if (c->kick2 != NULL) scheduler_activate(s, c->kick2);
- if (c->timestep != NULL) scheduler_activate(s, c->timestep);
- if (c->logger != NULL) scheduler_activate(s, c->logger);
+ if (c->nodeID == nodeID) {
+ if (cell_is_active_stars(c, e) ||
+ (with_star_formation && cell_is_active_hydro(c, e))) {
+
+ if (c->stars.ghost != NULL) scheduler_activate(s, c->stars.ghost);
+ if (c->stars.stars_in != NULL) scheduler_activate(s, c->stars.stars_in);
+ if (c->stars.stars_out != NULL) scheduler_activate(s, c->stars.stars_out);
+ if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
+ if (c->kick2 != NULL) scheduler_activate(s, c->kick2);
+ if (c->timestep != NULL) scheduler_activate(s, c->timestep);
+ if (c->logger != NULL) scheduler_activate(s, c->logger);
+ }
}
return rebuild;
@@ -3599,7 +3778,11 @@ int cell_unskip_black_holes_tasks(struct cell *c, struct scheduler *s) {
}
}
- else if (t->type == task_type_sub_pair || t->type == task_type_sub_self) {
+ else if (t->type == task_type_sub_self) {
+ cell_activate_subcell_black_holes_tasks(ci, NULL, s);
+ }
+
+ else if (t->type == task_type_sub_pair) {
cell_activate_subcell_black_holes_tasks(ci, cj, s);
}
}
@@ -4556,19 +4739,68 @@ void cell_drift_multipole(struct cell *c, const struct engine *e) {
* @brief Resets all the sorting properties for the stars in a given cell
* hierarchy.
*
+ * The clear_unused_flags argument can be used to additionally clean up all
+ * the flags demanding a sort for the given cell. This should be used with
+ * caution as it will prevent the sort tasks from doing anything on that cell
+ * until these flags are reset.
+ *
* @param c The #cell to clean.
+ * @param clear_unused_flags Do we also clean the flags demanding a sort?
*/
-void cell_clear_stars_sort_flags(struct cell *c) {
+void cell_clear_stars_sort_flags(struct cell *c, const int clear_unused_flags) {
+
+ /* Clear the flags that have not been reset by the sort task? */
+ if (clear_unused_flags) {
+ c->stars.requires_sorts = 0;
+ c->stars.do_sort = 0;
+ cell_clear_flag(c, cell_flag_do_stars_sub_sort);
+ }
+
+ /* Indicate that the cell is not sorted and cancel the pointer sorting arrays.
+ */
+ c->stars.sorted = 0;
+ cell_free_stars_sorts(c);
+
/* Recurse if possible */
if (c->split) {
for (int k = 0; k < 8; k++)
- if (c->progeny[k] != NULL) cell_clear_stars_sort_flags(c->progeny[k]);
+ if (c->progeny[k] != NULL)
+ cell_clear_stars_sort_flags(c->progeny[k], clear_unused_flags);
+ }
+}
+
+/**
+ * @brief Resets all the sorting properties for the hydro in a given cell
+ * hierarchy.
+ *
+ * The clear_unused_flags argument can be used to additionally clean up all
+ * the flags demanding a sort for the given cell. This should be used with
+ * caution as it will prevent the sort tasks from doing anything on that cell
+ * until these flags are reset.
+ *
+ * @param c The #cell to clean.
+ * @param clear_unused_flags Do we also clean the flags demanding a sort?
+ */
+void cell_clear_hydro_sort_flags(struct cell *c, const int clear_unused_flags) {
+
+ /* Clear the flags that have not been reset by the sort task? */
+ if (clear_unused_flags) {
+ c->hydro.do_sort = 0;
+ c->hydro.requires_sorts = 0;
+ cell_clear_flag(c, cell_flag_do_hydro_sub_sort);
}
/* Indicate that the cell is not sorted and cancel the pointer sorting arrays.
*/
- c->stars.sorted = 0;
- cell_free_stars_sorts(c);
+ c->hydro.sorted = 0;
+ cell_free_hydro_sorts(c);
+
+ /* Recurse if possible */
+ if (c->split) {
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL)
+ cell_clear_hydro_sort_flags(c->progeny[k], clear_unused_flags);
+ }
}
/**
@@ -4635,6 +4867,36 @@ void cell_check_spart_pos(const struct cell *c,
#endif
}
+/**
+ * @brief Checks that a cell and all its progenies have cleared their sort
+ * flags.
+ *
+ * Should only be used for debugging purposes.
+ *
+ * @param c The #cell to check.
+ */
+void cell_check_sort_flags(const struct cell *c) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ const int do_hydro_sub_sort = cell_get_flag(c, cell_flag_do_hydro_sub_sort);
+ const int do_stars_sub_sort = cell_get_flag(c, cell_flag_do_stars_sub_sort);
+
+ if (do_hydro_sub_sort)
+ error("cell %d has a hydro sub_sort flag set. Node=%d depth=%d maxdepth=%d",
+ c->cellID, c->nodeID, c->depth, c->maxdepth);
+
+ if (do_stars_sub_sort)
+ error("cell %d has a stars sub_sort flag set. Node=%d depth=%d maxdepth=%d",
+ c->cellID, c->nodeID, c->depth, c->maxdepth);
+
+ if (c->split) {
+ for (int k = 0; k < 8; ++k) {
+ if (c->progeny[k] != NULL) cell_check_sort_flags(c->progeny[k]);
+ }
+ }
+#endif
+}
+
/**
* @brief Recursively update the pointer and counter for #spart after the
* addition of a new particle.
@@ -5024,6 +5286,11 @@ struct spart *cell_convert_part_to_spart(struct engine *e, struct cell *c,
/* Did we run out of free spart slots? */
if (sp == NULL) return NULL;
+ /* Copy over the distance since rebuild */
+ sp->x_diff[0] = xp->x_diff[0];
+ sp->x_diff[1] = xp->x_diff[1];
+ sp->x_diff[2] = xp->x_diff[2];
+
/* Destroy the gas particle and get it's gpart friend */
struct gpart *gp = cell_convert_part_to_gpart(e, c, p, xp);
diff --git a/src/cell.h b/src/cell.h
index a98a1f6379e78ab23d5901269f8e90defbb1df09..4f8bed5fa5e337842bf37a0b382e3915a806a474 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -253,6 +253,26 @@ struct pcell_step_black_holes {
float dx_max_part;
};
+/**
+ * @brief Cell information to propagate the new counts of star particles.
+ */
+struct pcell_sf {
+
+ /*! Stars variables */
+ struct {
+
+ /* Distance by which the stars pointer has moved since the last rebuild */
+ ptrdiff_t delta_from_rebuild;
+
+ /* Number of particles in the cell */
+ int count;
+
+ /*! Maximum part movement in this cell since last construction. */
+ float dx_max_part;
+
+ } stars;
+};
+
/** Bitmasks for the cell flags. Beware when adding flags that you don't exceed
the size of the flags variable in the struct cell. */
enum cell_flags {
@@ -521,6 +541,9 @@ struct cell {
/*! Pointer to the #spart data. */
struct spart *parts;
+ /*! Pointer to the #spart data at rebuild time. */
+ struct spart *parts_rebuild;
+
/*! The star ghost task itself */
struct task *ghost;
@@ -802,6 +825,8 @@ int cell_unpack_end_step_black_holes(struct cell *c,
struct pcell_step_black_holes *pcell);
int cell_pack_multipoles(struct cell *c, struct gravity_tensors *m);
int cell_unpack_multipoles(struct cell *c, struct gravity_tensors *m);
+int cell_pack_sf_counts(struct cell *c, struct pcell_sf *pcell);
+int cell_unpack_sf_counts(struct cell *c, struct pcell_sf *pcell);
int cell_getsize(struct cell *c);
int cell_link_parts(struct cell *c, struct part *parts);
int cell_link_gparts(struct cell *c, struct gpart *gparts);
@@ -822,7 +847,8 @@ void cell_check_spart_drift_point(struct cell *c, void *data);
void cell_check_multipole_drift_point(struct cell *c, void *data);
void cell_reset_task_counters(struct cell *c);
int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s);
-int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s);
+int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s,
+ const int with_star_formation);
int cell_unskip_black_holes_tasks(struct cell *c, struct scheduler *s);
int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s);
void cell_drift_part(struct cell *c, const struct engine *e, int force);
@@ -838,7 +864,8 @@ void cell_activate_subcell_hydro_tasks(struct cell *ci, struct cell *cj,
void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
struct scheduler *s);
void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
- struct scheduler *s);
+ struct scheduler *s,
+ const int with_star_formation);
void cell_activate_subcell_black_holes_tasks(struct cell *ci, struct cell *cj,
struct scheduler *s);
void cell_activate_subcell_external_grav_tasks(struct cell *ci,
@@ -856,7 +883,9 @@ void cell_clear_limiter_flags(struct cell *c, void *data);
void cell_set_super_mapper(void *map_data, int num_elements, void *extra_data);
void cell_check_spart_pos(const struct cell *c,
const struct spart *global_sparts);
-void cell_clear_stars_sort_flags(struct cell *c);
+void cell_check_sort_flags(const struct cell *c);
+void cell_clear_stars_sort_flags(struct cell *c, const int unused_flags);
+void cell_clear_hydro_sort_flags(struct cell *c, const int unused_flags);
int cell_has_tasks(struct cell *c);
void cell_remove_part(const struct engine *e, struct cell *c, struct part *p,
struct xpart *xp);
diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c
index 0b5fcde02e4005c1d8c6a3bc4ff830e435b35ad6..3fda0c80040b6ad139b9100064b937e695e6030f 100644
--- a/src/cooling/EAGLE/cooling.c
+++ b/src/cooling/EAGLE/cooling.c
@@ -316,11 +316,14 @@ INLINE static double bisection_iter(
eagle_cooling_rate(log10(u_next_cgs), redshift, n_H_cgs,
abundance_ratio, n_H_index, d_n_H,
He_index, d_He, cooling);
-
- // Debugging
+#ifdef SWIFT_DEBUG_CHECKS
if (u_next_cgs <= 0)
- error("u_next_cgs %.5e u_upper %.5e u_lower %.5e Lambda %.5e", u_next_cgs,
- u_upper_cgs, u_lower_cgs, LambdaNet_cgs);
+ error(
+ "Got negative energy! u_next_cgs=%.5e u_upper=%.5e u_lower=%.5e "
+ "Lambda=%.5e",
+ u_next_cgs, u_upper_cgs, u_lower_cgs, LambdaNet_cgs);
+#endif
+
/* Where do we go next? */
if (u_next_cgs - u_ini_cgs - LambdaNet_cgs * ratefact_cgs * dt_cgs > 0.0) {
u_upper_cgs = u_next_cgs;
diff --git a/src/engine.c b/src/engine.c
index fe4be2a97d70802f37d6317109a0383a7a7f4e3e..253d8657ad0ac5231860452ff17cf68662d9a272 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2117,7 +2117,8 @@ void engine_allocate_foreign_particles(struct engine *e) {
}
/* For stars, we just use the numbers in the top-level cells */
- count_sparts_in += e->proxies[k].cells_in[j]->stars.count;
+ count_sparts_in +=
+ e->proxies[k].cells_in[j]->stars.count + space_extra_sparts;
/* For black holes, we just use the numbers in the top-level cells */
count_bparts_in += e->proxies[k].cells_in[j]->black_holes.count;
@@ -2208,7 +2209,8 @@ void engine_allocate_foreign_particles(struct engine *e) {
/* For stars, we just use the numbers in the top-level cells */
cell_link_sparts(e->proxies[k].cells_in[j], sparts);
- sparts = &sparts[e->proxies[k].cells_in[j]->stars.count];
+ sparts =
+ &sparts[e->proxies[k].cells_in[j]->stars.count + space_extra_sparts];
/* For black holes, we just use the numbers in the top-level cells */
cell_link_bparts(e->proxies[k].cells_in[j], bparts);
@@ -2563,6 +2565,8 @@ void engine_rebuild(struct engine *e, int repartitioned,
for (int k = 0; k < e->s->nr_local_cells; k++)
cell_check_foreign_multipole(&e->s->cells_top[e->s->local_cells_top[k]]);
}
+
+ space_check_sort_flags(e->s);
#endif
/* Run through the tasks and mark as skip or not. */
@@ -3333,7 +3337,8 @@ void engine_skip_force_and_kick(struct engine *e) {
t->subtype == task_subtype_tend_gpart ||
t->subtype == task_subtype_tend_spart ||
t->subtype == task_subtype_tend_bpart ||
- t->subtype == task_subtype_rho || t->subtype == task_subtype_gpart)
+ t->subtype == task_subtype_rho || t->subtype == task_subtype_gpart ||
+ t->subtype == task_subtype_sf_counts)
t->skip = 1;
}
@@ -3850,6 +3855,7 @@ void engine_step(struct engine *e) {
#ifdef SWIFT_DEBUG_CHECKS
/* Make sure all woken-up particles have been processed */
space_check_limiter(e->s);
+ space_check_sort_flags(e->s);
#endif
/* Collect information about the next time-step */
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index a9a9371a471ec58ec7f258dc87e72e94c6b45582..e1b1dcbbb2ceddbcf2dd26cc0ead72d694094990 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -232,11 +232,14 @@ void engine_addtasks_send_hydro(struct engine *e, struct cell *ci,
* @param ci The sending #cell.
* @param cj Dummy cell containing the nodeID of the receiving node.
* @param t_feedback The send_feed #task, if it has already been created.
+ * @param t_sf_counts The send_sf_counts, if it has been created.
* @param t_ti The recv_ti_end #task, if it has already been created.
+ * @param with_star_formation Are we running with star formation on?
*/
void engine_addtasks_send_stars(struct engine *e, struct cell *ci,
struct cell *cj, struct task *t_feedback,
- struct task *t_ti) {
+ struct task *t_sf_counts, struct task *t_ti,
+ const int with_star_formation) {
#ifdef WITH_MPI
@@ -244,6 +247,19 @@ void engine_addtasks_send_stars(struct engine *e, struct cell *ci,
struct scheduler *s = &e->sched;
const int nodeID = cj->nodeID;
+ if (t_sf_counts == NULL && with_star_formation && ci->hydro.count > 0) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ci->depth != 0)
+ error(
+ "Attaching a sf_count task at a non-top level c->depth=%d "
+ "c->count=%d",
+ ci->depth, ci->hydro.count);
+#endif
+ t_sf_counts = scheduler_addtask(s, task_type_send, task_subtype_sf_counts,
+ ci->mpi.tag, 0, ci, cj);
+ scheduler_addunlock(s, ci->hydro.star_formation, t_sf_counts);
+ }
+
/* Check if any of the density tasks are for the target node. */
for (l = ci->stars.density; l != NULL; l = l->next)
if (l->t->ci->nodeID == nodeID ||
@@ -279,13 +295,17 @@ void engine_addtasks_send_stars(struct engine *e, struct cell *ci,
engine_addlink(e, &ci->mpi.send, t_feedback);
engine_addlink(e, &ci->mpi.send, t_ti);
+ if (with_star_formation && ci->hydro.count > 0) {
+ engine_addlink(e, &ci->mpi.send, t_sf_counts);
+ }
}
/* Recurse? */
if (ci->split)
for (int k = 0; k < 8; k++)
if (ci->progeny[k] != NULL)
- engine_addtasks_send_stars(e, ci->progeny[k], cj, t_feedback, t_ti);
+ engine_addtasks_send_stars(e, ci->progeny[k], cj, t_feedback,
+ t_sf_counts, t_ti, with_star_formation);
#else
error("SWIFT was not compiled with MPI support.");
@@ -459,14 +479,30 @@ void engine_addtasks_recv_hydro(struct engine *e, struct cell *c,
* @param e The #engine.
* @param c The foreign #cell.
* @param t_feedback The recv_feed #task, if it has already been created.
+ * @param t_sf_counts The recv_sf_counts, if it has been created.
* @param t_ti The recv_ti_end #task, if it has already been created.
+ * @param with_star_formation Are we running with star formation on?
*/
void engine_addtasks_recv_stars(struct engine *e, struct cell *c,
- struct task *t_feedback, struct task *t_ti) {
+ struct task *t_feedback,
+ struct task *t_sf_counts, struct task *t_ti,
+ const int with_star_formation) {
#ifdef WITH_MPI
struct scheduler *s = &e->sched;
+ if (t_sf_counts == NULL && with_star_formation && c->hydro.count > 0) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->depth != 0)
+ error(
+ "Attaching a sf_count task at a non-top level c->depth=%d "
+ "c->count=%d",
+ c->depth, c->hydro.count);
+#endif
+ t_sf_counts = scheduler_addtask(s, task_type_recv, task_subtype_sf_counts,
+ c->mpi.tag, 0, c, NULL);
+ }
+
/* Have we reached a level where there are any stars tasks ? */
if (t_feedback == NULL && c->stars.density != NULL) {
@@ -481,11 +517,20 @@ void engine_addtasks_recv_stars(struct engine *e, struct cell *c,
t_ti = scheduler_addtask(s, task_type_recv, task_subtype_tend_spart,
c->mpi.tag, 0, c, NULL);
+
+ if (with_star_formation && c->hydro.count > 0) {
+
+ /* Receive the stars only once the counts have been received */
+ scheduler_addunlock(s, t_sf_counts, t_feedback);
+ }
}
if (t_feedback != NULL) {
engine_addlink(e, &c->mpi.recv, t_feedback);
engine_addlink(e, &c->mpi.recv, t_ti);
+ if (with_star_formation && c->hydro.count > 0) {
+ engine_addlink(e, &c->mpi.recv, t_sf_counts);
+ }
#ifdef SWIFT_DEBUG_CHECKS
if (c->nodeID == e->nodeID) error("Local cell!");
@@ -507,7 +552,8 @@ void engine_addtasks_recv_stars(struct engine *e, struct cell *c,
if (c->split)
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL)
- engine_addtasks_recv_stars(e, c->progeny[k], t_feedback, t_ti);
+ engine_addtasks_recv_stars(e, c->progeny[k], t_feedback, t_sf_counts,
+ t_ti, with_star_formation);
#else
error("SWIFT was not compiled with MPI support.");
@@ -614,6 +660,7 @@ void engine_addtasks_recv_gravity(struct engine *e, struct cell *c,
scheduler_addunlock(s, l->t, t_ti);
}
}
+
/* Recurse? */
if (c->split)
for (int k = 0; k < 8; k++)
@@ -2410,7 +2457,7 @@ void engine_addtasks_send_mapper(void *map_data, int num_elements,
void *extra_data) {
struct engine *e = (struct engine *)extra_data;
- // const int with_limiter = (e->policy & engine_policy_limiter);
+ const int with_star_formation = (e->policy & engine_policy_star_formation);
struct cell_type_pair *cell_type_pairs = (struct cell_type_pair *)map_data;
for (int k = 0; k < num_elements; k++) {
@@ -2428,7 +2475,9 @@ void engine_addtasks_send_mapper(void *map_data, int num_elements,
/* Add the send tasks for the cells in the proxy that have a stars
* connection. */
if ((e->policy & engine_policy_feedback) && (type & proxy_cell_type_hydro))
- engine_addtasks_send_stars(e, ci, cj, /*t_feedback=*/NULL, /*t_ti=*/NULL);
+ engine_addtasks_send_stars(e, ci, cj, /*t_feedback=*/NULL,
+ /*t_sf_counts=*/NULL, /*t_ti=*/NULL,
+ with_star_formation);
/* Add the send tasks for the cells in the proxy that have a black holes
* connection. */
@@ -2449,7 +2498,7 @@ void engine_addtasks_recv_mapper(void *map_data, int num_elements,
void *extra_data) {
struct engine *e = (struct engine *)extra_data;
- // const int with_limiter = (e->policy & engine_policy_limiter);
+ const int with_star_formation = (e->policy & engine_policy_star_formation);
struct cell_type_pair *cell_type_pairs = (struct cell_type_pair *)map_data;
for (int k = 0; k < num_elements; k++) {
@@ -2459,12 +2508,15 @@ void engine_addtasks_recv_mapper(void *map_data, int num_elements,
/* Add the recv tasks for the cells in the proxy that have a hydro
* connection. */
if ((e->policy & engine_policy_hydro) && (type & proxy_cell_type_hydro))
- engine_addtasks_recv_hydro(e, ci, NULL, NULL, NULL, NULL);
+ engine_addtasks_recv_hydro(e, ci, /*t_xv=*/NULL, /*t_rho=*/NULL,
+ /*t_gradient=*/NULL, /*t_ti=*/NULL);
/* Add the recv tasks for the cells in the proxy that have a stars
* connection. */
if ((e->policy & engine_policy_feedback) && (type & proxy_cell_type_hydro))
- engine_addtasks_recv_stars(e, ci, NULL, NULL);
+ engine_addtasks_recv_stars(e, ci, /*t_feedback=*/NULL,
+ /*t_sf_counts=*/NULL, /*t_ti=*/NULL,
+ with_star_formation);
/* Add the recv tasks for the cells in the proxy that have a black holes
* connection. */
@@ -2476,7 +2528,7 @@ void engine_addtasks_recv_mapper(void *map_data, int num_elements,
* connection. */
if ((e->policy & engine_policy_self_gravity) &&
(type & proxy_cell_type_gravity))
- engine_addtasks_recv_gravity(e, ci, NULL, NULL);
+ engine_addtasks_recv_gravity(e, ci, /*t_grav=*/NULL, /*t_ti=*/NULL);
}
}
diff --git a/src/engine_marktasks.c b/src/engine_marktasks.c
index ba52dd48de43b1388735a440f2a6f842dfbc22c6..0971c95bab862326e6c45cfb664c5cccca9cc0f2 100644
--- a/src/engine_marktasks.c
+++ b/src/engine_marktasks.c
@@ -70,6 +70,10 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
struct engine *e = (struct engine *)((size_t *)extra_data)[0];
const int nodeID = e->nodeID;
const int with_limiter = e->policy & engine_policy_limiter;
+ const int with_star_formation = e->policy & engine_policy_star_formation;
+#ifdef WITH_MPI
+ const int with_feedback = e->policy & engine_policy_feedback;
+#endif
for (int ind = 0; ind < num_elements; ind++) {
@@ -84,11 +88,18 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
/* Local pointer. */
struct cell *ci = t->ci;
+#ifdef SWIFT_DEBUG_CHECKS
if (ci->nodeID != nodeID) error("Non-local self task found");
+#endif
+ const int ci_active_hydro = cell_is_active_hydro(ci, e);
+ const int ci_active_gravity = cell_is_active_gravity(ci, e);
+ const int ci_active_black_holes = cell_is_active_black_holes(ci, e);
+ const int ci_active_stars = cell_is_active_stars(ci, e) ||
+ (with_star_formation && ci_active_hydro);
/* Activate the hydro drift */
if (t_type == task_type_self && t_subtype == task_subtype_density) {
- if (cell_is_active_hydro(ci, e)) {
+ if (ci_active_hydro) {
scheduler_activate(s, t);
cell_activate_drift_part(ci, s);
if (with_limiter) cell_activate_limiter(ci, s);
@@ -98,7 +109,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
/* Store current values of dx_max and h_max. */
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_density) {
- if (cell_is_active_hydro(ci, e)) {
+ if (ci_active_hydro) {
scheduler_activate(s, t);
cell_activate_subcell_hydro_tasks(ci, NULL, s);
if (with_limiter) cell_activate_limiter(ci, s);
@@ -106,37 +117,37 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
}
else if (t_type == task_type_self && t_subtype == task_subtype_force) {
- if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
+ if (ci_active_hydro) scheduler_activate(s, t);
}
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_force) {
- if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
+ if (ci_active_hydro) scheduler_activate(s, t);
}
else if (t->type == task_type_self &&
t->subtype == task_subtype_limiter) {
- if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
+ if (ci_active_hydro) scheduler_activate(s, t);
}
else if (t->type == task_type_sub_self &&
t->subtype == task_subtype_limiter) {
- if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
+ if (ci_active_hydro) scheduler_activate(s, t);
}
else if (t_type == task_type_self && t_subtype == task_subtype_gradient) {
- if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
+ if (ci_active_hydro) scheduler_activate(s, t);
}
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_gradient) {
- if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
+ if (ci_active_hydro) scheduler_activate(s, t);
}
/* Activate the star density */
else if (t_type == task_type_self &&
t_subtype == task_subtype_stars_density) {
- if (cell_is_active_stars(ci, e)) {
+ if (ci_active_stars) {
scheduler_activate(s, t);
cell_activate_drift_part(ci, s);
cell_activate_drift_spart(ci, s);
@@ -146,28 +157,28 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
/* Store current values of dx_max and h_max. */
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_stars_density) {
- if (cell_is_active_stars(ci, e)) {
+ if (ci_active_stars) {
scheduler_activate(s, t);
- cell_activate_subcell_stars_tasks(ci, NULL, s);
+ cell_activate_subcell_stars_tasks(ci, NULL, s, with_star_formation);
}
}
else if (t_type == task_type_self &&
t_subtype == task_subtype_stars_feedback) {
- if (cell_is_active_stars(ci, e)) {
+ if (ci_active_stars) {
scheduler_activate(s, t);
}
}
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_stars_feedback) {
- if (cell_is_active_stars(ci, e)) scheduler_activate(s, t);
+ if (ci_active_stars) scheduler_activate(s, t);
}
/* Activate the black hole density */
else if (t_type == task_type_self &&
t_subtype == task_subtype_bh_density) {
- if (cell_is_active_black_holes(ci, e)) {
+ if (ci_active_black_holes) {
scheduler_activate(s, t);
cell_activate_drift_part(ci, s);
cell_activate_drift_bpart(ci, s);
@@ -177,7 +188,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
/* Store current values of dx_max and h_max. */
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_bh_density) {
- if (cell_is_active_black_holes(ci, e)) {
+ if (ci_active_black_holes) {
scheduler_activate(s, t);
cell_activate_subcell_black_holes_tasks(ci, NULL, s);
}
@@ -185,19 +196,19 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
else if (t_type == task_type_self &&
t_subtype == task_subtype_bh_feedback) {
- if (cell_is_active_black_holes(ci, e)) {
+ if (ci_active_black_holes) {
scheduler_activate(s, t);
}
}
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_bh_feedback) {
- if (cell_is_active_black_holes(ci, e)) scheduler_activate(s, t);
+ if (ci_active_black_holes) scheduler_activate(s, t);
}
/* Activate the gravity drift */
else if (t_type == task_type_self && t_subtype == task_subtype_grav) {
- if (cell_is_active_gravity(ci, e)) {
+ if (ci_active_gravity) {
scheduler_activate(s, t);
cell_activate_subcell_grav_tasks(t->ci, NULL, s);
}
@@ -206,7 +217,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
/* Activate the gravity drift */
else if (t_type == task_type_self &&
t_subtype == task_subtype_external_grav) {
- if (cell_is_active_gravity(ci, e)) {
+ if (ci_active_gravity) {
scheduler_activate(s, t);
cell_activate_drift_gpart(t->ci, s);
}
@@ -238,12 +249,14 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
const int ci_active_gravity = cell_is_active_gravity(ci, e);
const int cj_active_gravity = cell_is_active_gravity(cj, e);
- const int ci_active_stars = cell_is_active_stars(ci, e);
- const int cj_active_stars = cell_is_active_stars(cj, e);
-
const int ci_active_black_holes = cell_is_active_black_holes(ci, e);
const int cj_active_black_holes = cell_is_active_black_holes(cj, e);
+ const int ci_active_stars = cell_is_active_stars(ci, e) ||
+ (with_star_formation && ci_active_hydro);
+ const int cj_active_stars = cell_is_active_stars(cj, e) ||
+ (with_star_formation && cj_active_hydro);
+
/* Only activate tasks that involve a local active cell. */
if ((t_subtype == task_subtype_density ||
t_subtype == task_subtype_gradient ||
@@ -338,7 +351,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
/* Store current values of dx_max and h_max. */
else if (t_type == task_type_sub_pair &&
t_subtype == task_subtype_stars_density) {
- cell_activate_subcell_stars_tasks(ci, cj, s);
+ cell_activate_subcell_stars_tasks(ci, cj, s, with_star_formation);
}
}
@@ -492,6 +505,20 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
scheduler_activate_send(s, cj->mpi.send, task_subtype_tend_part,
ci_nodeID);
+ /* Propagating new star counts? */
+ if (with_star_formation && with_feedback) {
+ if (ci_active_hydro && ci->hydro.count > 0) {
+ scheduler_activate_recv(s, ci->mpi.recv, task_subtype_sf_counts);
+ scheduler_activate_recv(s, ci->mpi.recv, task_subtype_tend_spart);
+ }
+ if (cj_active_hydro && cj->hydro.count > 0) {
+ scheduler_activate_send(s, cj->mpi.send, task_subtype_sf_counts,
+ ci_nodeID);
+ scheduler_activate_send(s, cj->mpi.send, task_subtype_tend_spart,
+ ci_nodeID);
+ }
+ }
+
} else if (cj_nodeID != nodeID) {
/* If the local cell is active, receive data from the foreign cell. */
@@ -538,6 +565,20 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
if (ci_active_hydro)
scheduler_activate_send(s, ci->mpi.send, task_subtype_tend_part,
cj_nodeID);
+
+ /* Propagating new star counts? */
+ if (with_star_formation && with_feedback) {
+ if (cj_active_hydro && cj->hydro.count > 0) {
+ scheduler_activate_recv(s, cj->mpi.recv, task_subtype_sf_counts);
+ scheduler_activate_recv(s, cj->mpi.recv, task_subtype_tend_spart);
+ }
+ if (ci_active_hydro && ci->hydro.count > 0) {
+ scheduler_activate_send(s, ci->mpi.send, task_subtype_sf_counts,
+ cj_nodeID);
+ scheduler_activate_send(s, ci->mpi.send, task_subtype_tend_spart,
+ cj_nodeID);
+ }
+ }
}
#endif
}
@@ -830,12 +871,16 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
/* Star ghost tasks ? */
else if (t_type == task_type_stars_ghost) {
- if (cell_is_active_stars(t->ci, e)) scheduler_activate(s, t);
+ if (cell_is_active_stars(t->ci, e) ||
+ (with_star_formation && cell_is_active_hydro(t->ci, e)))
+ scheduler_activate(s, t);
}
/* Feedback implicit tasks? */
else if (t_type == task_type_stars_in || t_type == task_type_stars_out) {
- if (cell_is_active_stars(t->ci, e)) scheduler_activate(s, t);
+ if (cell_is_active_stars(t->ci, e) ||
+ (with_star_formation && cell_is_active_hydro(t->ci, e)))
+ scheduler_activate(s, t);
}
/* Black hole ghost tasks ? */
diff --git a/src/feedback/EAGLE/feedback.c b/src/feedback/EAGLE/feedback.c
index c6bffce2cae2ca6b536cf7046076d38ef3ab544a..f7e714b21efdf6987e9a683a0873211eea8e13b7 100644
--- a/src/feedback/EAGLE/feedback.c
+++ b/src/feedback/EAGLE/feedback.c
@@ -138,14 +138,18 @@ INLINE static void compute_SNII_feedback(
const float ngb_gas_mass, const struct feedback_props* feedback_props) {
/* Time after birth considered for SNII feedback (internal units) */
- const float SNII_wind_delay = feedback_props->SNII_wind_delay;
+ const double SNII_wind_delay = feedback_props->SNII_wind_delay;
/* Are we doing feedback this step? */
if (star_age <= SNII_wind_delay && (star_age + dt) > SNII_wind_delay) {
+ if (sp->f_E != -1.f) {
#ifdef SWIFT_DEBUG_CHECKS
- if (sp->f_E != -1.f) error("Star has already done feedback!");
+ message("Star has already done feedback! sp->id=%lld age=%e d=%e", sp->id,
+ star_age, dt);
#endif
+ return;
+ }
/* Properties of the model (all in internal units) */
const double delta_T =
@@ -685,8 +689,8 @@ INLINE static void evolve_AGB(const float log10_min_mass, float log10_max_mass,
*/
void compute_stellar_evolution(const struct feedback_props* feedback_props,
const struct cosmology* cosmo, struct spart* sp,
- const struct unit_system* us, const float age,
- const float dt) {
+ const struct unit_system* us, const double age,
+ const double dt) {
TIMER_TIC;
@@ -843,9 +847,9 @@ void feedback_props_init(struct feedback_props* fp,
/* Properties of the SNII energy feedback model ------------------------- */
/* Set the delay time before SNII occur */
- const float Gyr_in_cgs = 1e9 * 365 * 24 * 3600;
+ const double Gyr_in_cgs = 1.0e9 * 365. * 24. * 3600.;
fp->SNII_wind_delay =
- parser_get_param_float(params, "EAGLEFeedback:SNII_wind_delay_Gyr") *
+ parser_get_param_double(params, "EAGLEFeedback:SNII_wind_delay_Gyr") *
Gyr_in_cgs / units_cgs_conversion_factor(us, UNIT_CONV_TIME);
/* Read the temperature change to use in stochastic heating */
diff --git a/src/feedback/EAGLE/feedback.h b/src/feedback/EAGLE/feedback.h
index 4f341a00fcc333190210ab446c738090647855ed..d6237adae9f87ce75b1a4f56f3d358fdc75d103d 100644
--- a/src/feedback/EAGLE/feedback.h
+++ b/src/feedback/EAGLE/feedback.h
@@ -30,8 +30,8 @@
void compute_stellar_evolution(const struct feedback_props* feedback_props,
const struct cosmology* cosmo, struct spart* sp,
- const struct unit_system* us, const float age,
- const float dt);
+ const struct unit_system* us, const double age,
+ const double dt);
/**
* @brief Should we do feedback for this star?
diff --git a/src/feedback/EAGLE/feedback_iact.h b/src/feedback/EAGLE/feedback_iact.h
index 63d5332f4d7eed815d7e32119d3ce6db67eb3af6..f1096beec2e1beb92fcdd92ed3f47da49d9fe439 100644
--- a/src/feedback/EAGLE/feedback_iact.h
+++ b/src/feedback/EAGLE/feedback_iact.h
@@ -116,7 +116,8 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx,
#ifdef SWIFT_DEBUG_CHECKS
if (Omega_frac < 0. || Omega_frac > 1.)
- error("Invalid fraction of material to dsitribute.");
+ error("Invalid fraction of material to distribute. Omega_frac=%e",
+ Omega_frac);
#endif
/* Update particle mass */
diff --git a/src/feedback/EAGLE/feedback_properties.h b/src/feedback/EAGLE/feedback_properties.h
index 3e6a90fb954ffaf2401bdd5bfe41a571b47254a7..931d5146b6d648a6b083830662bd6b1e57635b2a 100644
--- a/src/feedback/EAGLE/feedback_properties.h
+++ b/src/feedback/EAGLE/feedback_properties.h
@@ -201,7 +201,7 @@ struct feedback_props {
float num_SNII_per_msun;
/*! Wind delay time for SNII */
- float SNII_wind_delay;
+ double SNII_wind_delay;
/*! Temperature increase induced by SNe feedback */
float SNe_deltaT_desired;
diff --git a/src/runner.c b/src/runner.c
index 19b837812bde3f73a291e44c8978b637c39c3a39..061698da118e3afb50467020a1d3494d9418ea16 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -171,6 +171,11 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
TIMER_TIC;
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->nodeID != e->nodeID)
+ error("Running the star ghost on a foreign node!");
+#endif
+
/* Anything to do here? */
if (c->stars.count == 0) return;
if (!cell_is_active_stars(c, e)) return;
@@ -1095,6 +1100,9 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
/* Did we get a star? (Or did we run out of spare ones?) */
if (sp != NULL) {
+ /* message("We formed a star id=%lld cellID=%d", sp->id,
+ * c->cellID); */
+
/* Copy the properties of the gas particle to the star particle */
star_formation_copy_properties(p, xp, sp, e, sf_props, cosmo,
with_cosmology, phys_const,
@@ -1128,7 +1136,7 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
if (with_feedback && (c == c->top) &&
(current_stars_count != c->stars.count)) {
- cell_clear_stars_sort_flags(c);
+ cell_clear_stars_sort_flags(c, /*clear_unused_flags=*/0);
runner_do_all_stars_sort(r, c);
}
@@ -1303,15 +1311,25 @@ void runner_do_hydro_sort(struct runner *r, struct cell *c, int flags,
float dx_max_sort = 0.0f;
float dx_max_sort_old = 0.0f;
for (int k = 0; k < 8; k++) {
- if (c->progeny[k] != NULL && c->progeny[k]->hydro.count > 0) {
- /* Only propagate cleanup if the progeny is stale. */
- runner_do_hydro_sort(r, c->progeny[k], flags,
- cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
- space_maxreldx * c->progeny[k]->dmin),
- 0);
- dx_max_sort = max(dx_max_sort, c->progeny[k]->hydro.dx_max_sort);
- dx_max_sort_old =
- max(dx_max_sort_old, c->progeny[k]->hydro.dx_max_sort_old);
+ if (c->progeny[k] != NULL) {
+
+ if (c->progeny[k]->hydro.count > 0) {
+
+ /* Only propagate cleanup if the progeny is stale. */
+ runner_do_hydro_sort(
+ r, c->progeny[k], flags,
+ cleanup && (c->progeny[k]->hydro.dx_max_sort_old >
+ space_maxreldx * c->progeny[k]->dmin),
+ 0);
+ dx_max_sort = max(dx_max_sort, c->progeny[k]->hydro.dx_max_sort);
+ dx_max_sort_old =
+ max(dx_max_sort_old, c->progeny[k]->hydro.dx_max_sort_old);
+ } else {
+
+ /* We need to clean up the unused flags that were in case the
+ number of particles in the cell would change */
+ cell_clear_hydro_sort_flags(c->progeny[k], /*clear_unused_flags=*/1);
+ }
}
}
c->hydro.dx_max_sort = dx_max_sort;
@@ -1526,15 +1544,24 @@ void runner_do_stars_sort(struct runner *r, struct cell *c, int flags,
float dx_max_sort = 0.0f;
float dx_max_sort_old = 0.0f;
for (int k = 0; k < 8; k++) {
- if (c->progeny[k] != NULL && c->progeny[k]->stars.count > 0) {
- /* Only propagate cleanup if the progeny is stale. */
- const int cleanup_prog =
- cleanup && (c->progeny[k]->stars.dx_max_sort_old >
- space_maxreldx * c->progeny[k]->dmin);
- runner_do_stars_sort(r, c->progeny[k], flags, cleanup_prog, 0);
- dx_max_sort = max(dx_max_sort, c->progeny[k]->stars.dx_max_sort);
- dx_max_sort_old =
- max(dx_max_sort_old, c->progeny[k]->stars.dx_max_sort_old);
+ if (c->progeny[k] != NULL) {
+
+ if (c->progeny[k]->stars.count > 0) {
+
+ /* Only propagate cleanup if the progeny is stale. */
+ const int cleanup_prog =
+ cleanup && (c->progeny[k]->stars.dx_max_sort_old >
+ space_maxreldx * c->progeny[k]->dmin);
+ runner_do_stars_sort(r, c->progeny[k], flags, cleanup_prog, 0);
+ dx_max_sort = max(dx_max_sort, c->progeny[k]->stars.dx_max_sort);
+ dx_max_sort_old =
+ max(dx_max_sort_old, c->progeny[k]->stars.dx_max_sort_old);
+ } else {
+
+ /* We need to clean up the unused flags that were in case the
+ number of particles in the cell would change */
+ cell_clear_stars_sort_flags(c->progeny[k], /*clear_unused_flags=*/1);
+ }
}
}
c->stars.dx_max_sort = dx_max_sort;
@@ -1694,6 +1721,8 @@ void runner_do_all_hydro_sort(struct runner *r, struct cell *c) {
if (c->nodeID != engine_rank) error("Function called on a foreign cell!");
#endif
+ if (!cell_is_active_hydro(c, r->e)) return;
+
/* Shall we sort at this level? */
if (c->hydro.super == c) {
@@ -2337,27 +2366,36 @@ static void runner_do_unskip_hydro(struct cell *c, struct engine *e) {
*
* @param c The cell.
* @param e The engine.
+ * @param with_star_formation Are we running with star formation switched on?
*/
-static void runner_do_unskip_stars(struct cell *c, struct engine *e) {
+static void runner_do_unskip_stars(struct cell *c, struct engine *e,
+ const int with_star_formation) {
+
+ const int non_empty =
+ c->stars.count > 0 || (with_star_formation && c->hydro.count > 0);
/* Ignore empty cells. */
- if (c->stars.count == 0) return;
+ if (!non_empty) return;
+
+ const int ci_active = cell_is_active_stars(c, e) ||
+ (with_star_formation && cell_is_active_hydro(c, e));
/* Skip inactive cells. */
- if (!cell_is_active_stars(c, e)) return;
+ if (!ci_active) return;
/* Recurse */
if (c->split) {
for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL) {
struct cell *cp = c->progeny[k];
- runner_do_unskip_stars(cp, e);
+ runner_do_unskip_stars(cp, e, with_star_formation);
}
}
}
/* Unskip any active tasks. */
- const int forcerebuild = cell_unskip_stars_tasks(c, &e->sched);
+ const int forcerebuild =
+ cell_unskip_stars_tasks(c, &e->sched, with_star_formation);
if (forcerebuild) atomic_inc(&e->forcerebuild);
}
@@ -2429,6 +2467,7 @@ void runner_do_unskip_mapper(void *map_data, int num_elements,
void *extra_data) {
struct engine *e = (struct engine *)extra_data;
+ const int with_star_formation = e->policy & engine_policy_star_formation;
const int nodeID = e->nodeID;
struct space *s = e->s;
int *local_cells = (int *)map_data;
@@ -2446,7 +2485,8 @@ void runner_do_unskip_mapper(void *map_data, int num_elements,
runner_do_unskip_gravity(c, e);
/* Stars tasks */
- if (e->policy & engine_policy_stars) runner_do_unskip_stars(c, e);
+ if (e->policy & engine_policy_stars)
+ runner_do_unskip_stars(c, e, with_star_formation);
/* Black hole tasks */
if (e->policy & engine_policy_black_holes)
@@ -2547,7 +2587,7 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
/* Anything to do here? */
if (!cell_is_starting_hydro(c, e) && !cell_is_starting_gravity(c, e) &&
- !cell_is_starting_stars(c, e))
+ !cell_is_starting_stars(c, e) && !cell_is_starting_black_holes(c, e))
return;
/* Recurse? */
@@ -2731,7 +2771,7 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
/* Anything to do here? */
if (!cell_is_active_hydro(c, e) && !cell_is_active_gravity(c, e) &&
- !cell_is_active_stars(c, e))
+ !cell_is_active_stars(c, e) && !cell_is_active_black_holes(c, e))
return;
/* Recurse? */
@@ -3818,14 +3858,15 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int clear_sorts,
ti_stars_end_min =
min(ti_stars_end_min, c->progeny[k]->stars.ti_end_min);
ti_stars_end_max =
- max(ti_stars_end_max, c->progeny[k]->grav.ti_end_max);
+ max(ti_stars_end_max, c->progeny[k]->stars.ti_end_max);
h_max = max(h_max, c->progeny[k]->stars.h_max);
}
}
}
#ifdef SWIFT_DEBUG_CHECKS
- if (ti_stars_end_min < ti_current)
+ if (ti_stars_end_min < ti_current &&
+ !(r->e->policy & engine_policy_star_formation))
error(
"Received a cell at an incorrect time c->ti_end_min=%lld, "
"e->ti_current=%lld.",
@@ -3904,7 +3945,7 @@ void runner_do_recv_bpart(struct runner *r, struct cell *c, int clear_sorts,
ti_black_holes_end_min =
min(ti_black_holes_end_min, c->progeny[k]->black_holes.ti_end_min);
ti_black_holes_end_max =
- max(ti_black_holes_end_max, c->progeny[k]->grav.ti_end_max);
+ max(ti_black_holes_end_max, c->progeny[k]->black_holes.ti_end_max);
h_max = max(h_max, c->progeny[k]->black_holes.h_max);
}
}
@@ -4162,14 +4203,13 @@ void *runner_main(void *data) {
case task_type_send:
if (t->subtype == task_subtype_tend_part) {
free(t->buff);
- }
- if (t->subtype == task_subtype_tend_gpart) {
+ } else if (t->subtype == task_subtype_tend_gpart) {
free(t->buff);
- }
- if (t->subtype == task_subtype_tend_spart) {
+ } else if (t->subtype == task_subtype_tend_spart) {
free(t->buff);
- }
- if (t->subtype == task_subtype_tend_bpart) {
+ } else if (t->subtype == task_subtype_tend_bpart) {
+ free(t->buff);
+ } else if (t->subtype == task_subtype_sf_counts) {
free(t->buff);
}
break;
@@ -4187,6 +4227,10 @@ void *runner_main(void *data) {
cell_unpack_end_step_black_holes(
ci, (struct pcell_step_black_holes *)t->buff);
free(t->buff);
+ } else if (t->subtype == task_subtype_sf_counts) {
+ cell_unpack_sf_counts(ci, (struct pcell_sf *)t->buff);
+ cell_clear_stars_sort_flags(ci, /*clear_unused_flags=*/0);
+ free(t->buff);
} else if (t->subtype == task_subtype_xv) {
runner_do_recv_part(r, ci, 1, 1);
} else if (t->subtype == task_subtype_rho) {
diff --git a/src/scheduler.c b/src/scheduler.c
index 2c46fd87e4076e9628792504d6ce1ca5d0b738e7..6a58d657105a28e6882275c21b81125253db392c 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -1665,6 +1665,13 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
err = MPI_Irecv(t->buff, t->ci->mpi.pcell_size, multipole_mpi_type,
t->ci->nodeID, t->flags, subtaskMPI_comms[t->subtype],
&t->req);
+ } else if (t->subtype == task_subtype_sf_counts) {
+ t->buff = (struct pcell_sf *)malloc(sizeof(struct pcell_sf) *
+ t->ci->mpi.pcell_size);
+ err = MPI_Irecv(t->buff,
+ t->ci->mpi.pcell_size * sizeof(struct pcell_sf),
+ MPI_BYTE, t->ci->nodeID, t->flags,
+ subtaskMPI_comms[t->subtype], &t->req);
} else {
error("Unknown communication sub-type");
}
@@ -1800,6 +1807,14 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
err = MPI_Isend(t->buff, t->ci->mpi.pcell_size, multipole_mpi_type,
t->cj->nodeID, t->flags, subtaskMPI_comms[t->subtype],
&t->req);
+ } else if (t->subtype == task_subtype_sf_counts) {
+ t->buff = (struct pcell_sf *)malloc(sizeof(struct pcell_sf) *
+ t->ci->mpi.pcell_size);
+ cell_pack_sf_counts(t->ci, (struct pcell_sf *)t->buff);
+ err = MPI_Isend(t->buff,
+ t->ci->mpi.pcell_size * sizeof(struct pcell_sf),
+ MPI_BYTE, t->cj->nodeID, t->flags,
+ subtaskMPI_comms[t->subtype], &t->req);
} else {
error("Unknown communication sub-type");
}
diff --git a/src/space.c b/src/space.c
index 1712e512e5b65bfbe71231fc2cbaba20447a5a8e..1c5023c3b79aab472b40756072c01870eb50dcd4 100644
--- a/src/space.c
+++ b/src/space.c
@@ -266,6 +266,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->hydro.xparts = NULL;
c->grav.parts = NULL;
c->stars.parts = NULL;
+ c->stars.parts_rebuild = NULL;
c->black_holes.parts = NULL;
c->flags = 0;
c->hydro.ti_end_min = -1;
@@ -1874,6 +1875,9 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
c->stars.parts = sfinger;
c->black_holes.parts = bfinger;
+ /* Store the state at rebuild time */
+ c->stars.parts_rebuild = c->stars.parts;
+
c->hydro.count_total = c->hydro.count + space_extra_parts;
c->grav.count_total = c->grav.count + space_extra_gparts;
c->stars.count_total = c->stars.count + space_extra_sparts;
@@ -3700,8 +3704,8 @@ void space_split_mapper(void *map_data, int num_cells, void *extra_data) {
void space_recycle(struct space *s, struct cell *c) {
/* Clear the cell. */
- if (lock_destroy(&c->lock) != 0 || lock_destroy(&c->grav.plock) != 0 ||
- lock_destroy(&c->mlock) != 0 || lock_destroy(&c->stars.lock) != 0 ||
+ if (lock_destroy(&c->hydro.lock) != 0 || lock_destroy(&c->grav.plock) != 0 ||
+ lock_destroy(&c->grav.mlock) != 0 || lock_destroy(&c->stars.lock) != 0 ||
lock_destroy(&c->black_holes.lock) != 0 ||
lock_destroy(&c->stars.star_formation_lock))
error("Failed to destroy spinlocks.");
@@ -3751,8 +3755,10 @@ void space_recycle_list(struct space *s, struct cell *cell_list_begin,
/* Clean up the list of cells. */
for (struct cell *c = cell_list_begin; c != NULL; c = c->next) {
/* Clear the cell. */
- if (lock_destroy(&c->lock) != 0 || lock_destroy(&c->grav.plock) != 0 ||
- lock_destroy(&c->mlock) != 0 || lock_destroy(&c->stars.lock) != 0 ||
+ if (lock_destroy(&c->hydro.lock) != 0 ||
+ lock_destroy(&c->grav.plock) != 0 ||
+ lock_destroy(&c->grav.mlock) != 0 ||
+ lock_destroy(&c->stars.lock) != 0 ||
lock_destroy(&c->black_holes.lock) != 0 ||
lock_destroy(&c->stars.star_formation_lock))
error("Failed to destroy spinlocks.");
@@ -3899,6 +3905,7 @@ void space_list_useful_top_level_cells(struct space *s) {
const int has_particles =
(c->hydro.count > 0) || (c->grav.count > 0) || (c->stars.count > 0) ||
+ (c->black_holes.count > 0) ||
(c->grav.multipole != NULL && c->grav.multipole->m_pole.M_000 > 0.f);
if (has_particles) {
@@ -5185,6 +5192,41 @@ void space_check_limiter(struct space *s) {
#endif
}
+void space_check_sort_flags_mapper(void *map_data, int nr_cells,
+ void *extra_data) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+ const struct space *s = (struct space *)extra_data;
+ int *local_cells_top = map_data;
+
+ for (int ind = 0; ind < nr_cells; ++ind) {
+ const struct cell *c = &s->cells_top[local_cells_top[ind]];
+
+ cell_check_sort_flags(c);
+ }
+
+#endif
+}
+
+/**
+ * @brief Checks that all cells have cleared their sort flags.
+ *
+ * Should only be used for debugging purposes.
+ *
+ * @param s The #space to check.
+ */
+void space_check_sort_flags(struct space *s) {
+#ifdef SWIFT_DEBUG_CHECKS
+
+ threadpool_map(&s->e->threadpool, space_check_sort_flags_mapper,
+ s->local_cells_with_tasks_top, s->nr_local_cells_with_tasks,
+ sizeof(int), 1, s);
+#else
+ error("Calling debugging code without debugging flag activated.");
+#endif
+}
+
/**
* @brief Resets all the individual cell task counters to 0.
*
diff --git a/src/space.h b/src/space.h
index 4a2d5d8ce92d49ef129fc32c9332bc811e67f795..88f38209e7d8d82b021cd518dfa9806e4ec844eb 100644
--- a/src/space.h
+++ b/src/space.h
@@ -354,6 +354,7 @@ void space_check_top_multipoles_drift_point(struct space *s,
integertime_t ti_drift);
void space_check_timesteps(struct space *s);
void space_check_limiter(struct space *s);
+void space_check_sort_flags(struct space *s);
void space_replicate(struct space *s, int replicate, int verbose);
void space_generate_gas(struct space *s, const struct cosmology *cosmo,
int periodic, const double dim[3], int verbose);
diff --git a/src/task.c b/src/task.c
index ee2f9eb6b9cc74a28d581c4e4c42157d38b2e2ba..77731ae02f8c84597565717319241997a77fadf7 100644
--- a/src/task.c
+++ b/src/task.c
@@ -101,8 +101,8 @@ const char *subtaskID_names[task_subtype_count] = {
"limiter", "grav", "external_grav", "tend_part",
"tend_gpart", "tend_spart", "tend_bpart", "xv",
"rho", "gpart", "multipole", "spart",
- "stars_density", "stars_feedback", "bpart", "bh_density",
- "bh_feedback"};
+ "stars_density", "stars_feedback", "sf_count", "bpart",
+ "bh_density", "bh_feedback"};
#ifdef WITH_MPI
/* MPI communicators for the subtypes. */
@@ -407,9 +407,8 @@ void task_unlock(struct task *t) {
if (subtype == task_subtype_grav) {
cell_gunlocktree(ci);
cell_munlocktree(ci);
- } else if (subtype == task_subtype_stars_density) {
- cell_sunlocktree(ci);
- } else if (subtype == task_subtype_stars_feedback) {
+ } else if ((subtype == task_subtype_stars_density) ||
+ (subtype == task_subtype_stars_feedback)) {
cell_sunlocktree(ci);
cell_unlocktree(ci);
} else {
@@ -424,10 +423,8 @@ void task_unlock(struct task *t) {
cell_gunlocktree(cj);
cell_munlocktree(ci);
cell_munlocktree(cj);
- } else if (subtype == task_subtype_stars_density) {
- cell_sunlocktree(ci);
- cell_sunlocktree(cj);
- } else if (subtype == task_subtype_stars_feedback) {
+ } else if ((subtype == task_subtype_stars_density) ||
+ (subtype == task_subtype_stars_feedback)) {
cell_sunlocktree(ci);
cell_sunlocktree(cj);
cell_unlocktree(ci);
@@ -544,10 +541,8 @@ int task_lock(struct task *t) {
cell_gunlocktree(ci);
return 0;
}
- } else if (subtype == task_subtype_stars_density) {
- if (ci->stars.hold) return 0;
- if (cell_slocktree(ci) != 0) return 0;
- } else if (subtype == task_subtype_stars_feedback) {
+ } else if ((subtype == task_subtype_stars_density) ||
+ (subtype == task_subtype_stars_feedback)) {
if (ci->stars.hold) return 0;
if (ci->hydro.hold) return 0;
if (cell_slocktree(ci) != 0) return 0;
@@ -580,14 +575,8 @@ int task_lock(struct task *t) {
cell_munlocktree(ci);
return 0;
}
- } else if (subtype == task_subtype_stars_density) {
- if (ci->stars.hold || cj->stars.hold) return 0;
- if (cell_slocktree(ci) != 0) return 0;
- if (cell_slocktree(cj) != 0) {
- cell_sunlocktree(ci);
- return 0;
- }
- } else if (subtype == task_subtype_stars_feedback) {
+ } else if ((subtype == task_subtype_stars_density) ||
+ (subtype == task_subtype_stars_feedback)) {
/* Lock the stars and the gas particles in both cells */
if (ci->stars.hold || cj->stars.hold) return 0;
if (ci->hydro.hold || cj->hydro.hold) return 0;
diff --git a/src/task.h b/src/task.h
index 14d1f8518dde3fa574d01a3d500b9dedf53b8825..eef0aca5fb075725fb1da51ccaf4033a3608e078 100644
--- a/src/task.h
+++ b/src/task.h
@@ -114,6 +114,7 @@ enum task_subtypes {
task_subtype_spart,
task_subtype_stars_density,
task_subtype_stars_feedback,
+ task_subtype_sf_counts,
task_subtype_bpart,
task_subtype_bh_density,
task_subtype_bh_feedback,
diff --git a/tools/plot_task_dependencies.py b/tools/plot_task_dependencies.py
index 9a8b63ab315215c2973187a026fc1a5e182cd47b..b835b04261cf54b361f347f88bbb27f230a6be16 100644
--- a/tools/plot_task_dependencies.py
+++ b/tools/plot_task_dependencies.py
@@ -162,6 +162,9 @@ def taskIsStars(name):
"""
if "stars" in name or "spart" in name:
return True
+ if "sf_count" in name:
+ return True
+
return False
def taskIsHydro(name):
diff --git a/tools/task_plots/analyse_tasks.py b/tools/task_plots/analyse_tasks.py
index 463f4f1e9d792ae2f4103d80b5dc45761e3755e4..fd07340940b8903bf09fe807f2741d9e65261129 100755
--- a/tools/task_plots/analyse_tasks.py
+++ b/tools/task_plots/analyse_tasks.py
@@ -131,6 +131,7 @@ SUBTYPES = [
"spart",
"stars_density",
"stars_feedback",
+ "sf_counts"
"bpart",
"bh_density",
"bh_feedback",
diff --git a/tools/task_plots/plot_tasks.py b/tools/task_plots/plot_tasks.py
index ff12940e1788099dbd67cc11196f7a3e1261ab3c..070ae2f3c509d6ae8fef8bad9ffec6c316060a7c 100755
--- a/tools/task_plots/plot_tasks.py
+++ b/tools/task_plots/plot_tasks.py
@@ -216,6 +216,7 @@ SUBTYPES = [
"spart",
"stars_density",
"stars_feedback",
+ "sf_counts",
"bpart",
"bh_density",
"bh_feedback",
@@ -258,6 +259,8 @@ FULLTYPES = [
"send/gpart",
"recv/spart",
"send/spart",
+ "send/sf_counts",
+ "recv/sf_counts",
"recv/bpart",
"send/bpart",
"self/stars_density",