From 8a1c735c85f9924cdb41b56aef96d63194f0631d Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Sat, 4 May 2019 13:38:45 +0100
Subject: [PATCH] Activate star tasks also if we are running with SF on and the
 cell is active for hydro purposes.

---
 src/cell.c             | 86 ++++++++++++++++++++++++++++++------------
 src/cell.h             |  6 ++-
 src/engine_marktasks.c | 72 ++++++++++++++++++++++-------------
 src/runner.c           | 17 ++++++---
 4 files changed, 122 insertions(+), 59 deletions(-)

diff --git a/src/cell.c b/src/cell.c
index be593c7da4..ed56312126 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -2692,9 +2692,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. */
@@ -2712,21 +2714,24 @@ 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)
-      return;
+    if (!ci_active || ci->hydro.count == 0 || ci->stars.count == 0) return;
 
     /* Recurse? */
     if (cell_can_recurse_in_self_stars_task(ci)) {
       /* 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 {
@@ -2738,8 +2743,15 @@ void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
 
   /* Otherwise, pair interation */
   else {
+
+    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 (!cell_is_active_stars(ci, e) && !cell_is_active_stars(cj, e)) return;
+    if (!ci_active && !cj_active) return;
 
     /* Get the orientation of the pair. */
     double shift[3];
@@ -2754,13 +2766,14 @@ void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
         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);
@@ -2777,7 +2790,7 @@ 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);
@@ -3444,10 +3457,13 @@ 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;
@@ -3461,8 +3477,14 @@ 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;
+
+    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)));
+
 #ifdef WITH_MPI
     const int ci_nodeID = ci->nodeID;
     const int cj_nodeID = (cj != NULL) ? cj->nodeID : -1;
@@ -3522,8 +3544,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);
       }
     }
 
@@ -3606,8 +3632,14 @@ 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;
+
+    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)));
+
 #ifdef WITH_MPI
     const int ci_nodeID = ci->nodeID;
     const int cj_nodeID = (cj != NULL) ? cj->nodeID : -1;
@@ -3643,14 +3675,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;
diff --git a/src/cell.h b/src/cell.h
index 1cc6b81e35..461ed45028 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -844,7 +844,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);
@@ -860,7 +861,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,
diff --git a/src/engine_marktasks.c b/src/engine_marktasks.c
index e38e2efc41..5560714dbe 100644
--- a/src/engine_marktasks.c
+++ b/src/engine_marktasks.c
@@ -70,8 +70,9 @@ 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;
-#ifdef WITH_MPI
   const int with_star_formation = e->policy & engine_policy_star_formation;
+
+#ifdef WITH_MPI
   const int with_feedback = e->policy & engine_policy_feedback;
 #endif
 
@@ -88,11 +89,19 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
       /* Local pointer. */
       struct cell *ci = t->ci;
 
+      /* How active is this cell? */
+      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 && cell_is_active_hydro(ci, e));
+
       if (ci->nodeID != nodeID) error("Non-local self task found");
 
       /* 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);
@@ -102,7 +111,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);
@@ -110,37 +119,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);
@@ -150,28 +159,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);
@@ -181,7 +190,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);
         }
@@ -189,19 +198,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);
         }
@@ -210,7 +219,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);
         }
@@ -242,12 +251,16 @@ 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 && cell_is_active_hydro(ci, e));
+      const int cj_active_stars =
+          cell_is_active_stars(cj, e) ||
+          (with_star_formation && cell_is_active_hydro(cj, e));
+
       /* Only activate tasks that involve a local active cell. */
       if ((t_subtype == task_subtype_density ||
            t_subtype == task_subtype_gradient ||
@@ -342,7 +355,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);
         }
       }
 
@@ -825,7 +838,8 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
     /* logger tasks ? */
     else if (t->type == task_type_logger) {
       if (cell_is_active_hydro(t->ci, e) || cell_is_active_gravity(t->ci, e) ||
-          cell_is_active_stars(t->ci, e))
+          cell_is_active_stars(t->ci, e) ||
+          cell_is_active_black_holes(t->ci, e))
         scheduler_activate(s, t);
     }
 
@@ -862,12 +876,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/runner.c b/src/runner.c
index 8420153689..95e70aed72 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -2344,27 +2344,32 @@ 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 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) {
 
   /* Ignore empty cells. */
   if (c->stars.count == 0) return;
 
   /* Skip inactive cells. */
-  if (!cell_is_active_stars(c, e)) return;
+  if (!cell_is_active_stars(c, e) &&
+      (!with_star_formation || !cell_is_active_hydro(c, e)))
+    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);
 }
 
@@ -2436,6 +2441,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;
@@ -2453,7 +2459,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)
-- 
GitLab