Feedback task

src/cell.c

@@ -2771,7 +2771,7 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
 
       /* Check whether there was too much particle motion, i.e. the
          cell neighbour conditions were violated. */
-      if (cell_need_rebuild_for_pair(ci, cj)) rebuild = 1;
+      if (cell_need_rebuild_for_hydro_pair(ci, cj)) rebuild = 1;
 
 #ifdef WITH_MPI
       /* Activate the send/recv tasks. */
@@ -3125,7 +3125,7 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s) {
 
       /* Check whether there was too much particle motion, i.e. the
          cell neighbour conditions were violated. */
-      if (cell_need_rebuild_for_pair(ci, cj)) rebuild = 1;
+      if (cell_need_rebuild_for_stars_pair(ci, cj)) rebuild = 1;
 
 #ifdef WITH_MPI
       error("MPI with stars not implemented");
@@ -3212,6 +3212,10 @@ 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)) {
 
+    /* Un-skip the feedback tasks involved with this cell. */
+    for (struct link *l = c->stars.feedback; l != NULL; l = l->next)
+      scheduler_activate(s, l->t);
+
     if (c->stars.ghost_in != NULL) scheduler_activate(s, c->stars.ghost_in);
     if (c->stars.ghost_out != NULL) scheduler_activate(s, c->stars.ghost_out);
     if (c->stars.ghost != NULL) scheduler_activate(s, c->stars.ghost);

src/cell.h

@@ -469,6 +469,9 @@ struct cell {
     /*! Linked list of the tasks computing this cell's star density. */
     struct link *density;
 
+    /*! Linked list of the tasks computing this cell's star feedback. */
+    struct link *feedback;
+
     /*! The task computing this cell's sorts. */
     struct task *sorts;
 
@@ -952,14 +955,15 @@ cell_can_split_self_gravity_task(const struct cell *c) {
 }
 
 /**
- * @brief Have particles in a pair of cells moved too much and require a rebuild
+ * @brief Have gas particles in a pair of cells moved too much and require a
+ * rebuild
  * ?
  *
  * @param ci The first #cell.
  * @param cj The second #cell.
  */
-__attribute__((always_inline)) INLINE static int cell_need_rebuild_for_pair(
-    const struct cell *ci, const struct cell *cj) {
+__attribute__((always_inline)) INLINE static int
+cell_need_rebuild_for_hydro_pair(const struct cell *ci, const struct cell *cj) {
 
   /* Is the cut-off radius plus the max distance the parts in both cells have */
   /* moved larger than the cell size ? */
@@ -969,6 +973,25 @@ __attribute__((always_inline)) INLINE static int cell_need_rebuild_for_pair(
           cj->dmin);
 }
 
+/**
+ * @brief Have star particles in a pair of cells moved too much and require a
+ * rebuild
+ * ?
+ *
+ * @param ci The first #cell.
+ * @param cj The second #cell.
+ */
+__attribute__((always_inline)) INLINE static int
+cell_need_rebuild_for_stars_pair(const struct cell *ci, const struct cell *cj) {
+
+  /* Is the cut-off radius plus the max distance the parts in both cells have */
+  /* moved larger than the cell size ? */
+  /* Note ci->dmin == cj->dmin */
+  return (kernel_gamma * max(ci->stars.h_max, cj->stars.h_max) +
+              ci->stars.dx_max_part + cj->stars.dx_max_part >
+          cj->dmin);
+}
+
 /**
  * @brief Add a unique tag to a cell, mostly for MPI communications.
  *

src/engine_maketasks.c

@@ -1015,6 +1015,8 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
         engine_addlink(e, &ci->grav.grav, t);
       } else if (t->subtype == task_subtype_stars_density) {
         engine_addlink(e, &ci->stars.density, t);
+      } else if (t->subtype == task_subtype_stars_feedback) {
+        engine_addlink(e, &ci->stars.feedback, t);
       }
 
       /* Link pair tasks to cells. */
@@ -1031,6 +1033,9 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
       } else if (t->subtype == task_subtype_stars_density) {
         engine_addlink(e, &ci->stars.density, t);
         engine_addlink(e, &cj->stars.density, t);
+      } else if (t->subtype == task_subtype_stars_feedback) {
+        engine_addlink(e, &ci->stars.feedback, t);
+        engine_addlink(e, &cj->stars.feedback, t);
       }
 #ifdef SWIFT_DEBUG_CHECKS
       else if (t_subtype == task_subtype_external_grav) {
@@ -1050,6 +1055,8 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
         engine_addlink(e, &ci->grav.grav, t);
       } else if (t->subtype == task_subtype_stars_density) {
         engine_addlink(e, &ci->stars.density, t);
+      } else if (t->subtype == task_subtype_stars_feedback) {
+        engine_addlink(e, &ci->stars.feedback, t);
       }
 
       /* Link sub-pair tasks to cells. */
@@ -1066,6 +1073,9 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
       } else if (t->subtype == task_subtype_stars_density) {
         engine_addlink(e, &ci->stars.density, t);
         engine_addlink(e, &cj->stars.density, t);
+      } else if (t->subtype == task_subtype_stars_feedback) {
+        engine_addlink(e, &ci->stars.feedback, t);
+        engine_addlink(e, &cj->stars.feedback, t);
       }
 #ifdef SWIFT_DEBUG_CHECKS
       else if (t_subtype == task_subtype_external_grav) {
@@ -1305,9 +1315,11 @@ static inline void engine_make_hydro_loops_dependencies(struct scheduler *sched,
  */
 static inline void engine_make_stars_loops_dependencies(struct scheduler *sched,
                                                         struct task *density,
+                                                        struct task *feedback,
                                                         struct cell *c) {
-  /* density loop --> ghost */
+  /* density loop --> ghost --> feedback loop*/
   scheduler_addunlock(sched, density, c->super->stars.ghost_in);
+  scheduler_addunlock(sched, c->super->stars.ghost_out, feedback);
 }
 
 /**
@@ -1563,14 +1575,17 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 }
 
 /**
- * @brief Creates all the task dependencies for the stars
+ * @brief Duplicates the first stars loop and construct all the
+ * dependencies for the stars part
  *
- * @param map_data The tasks
- * @param num_elements number of tasks
- * @param extra_data The #engine
+ * This is done by looping over all the previously constructed tasks
+ * and adding another task involving the same cells but this time
+ * corresponding to the second stars loop over neighbours.
+ * With all the relevant tasks for a given cell available, we construct
+ * all the dependencies for that cell.
  */
-void engine_link_stars_tasks_mapper(void *map_data, int num_elements,
-                                    void *extra_data) {
+void engine_make_extra_starsloop_tasks_mapper(void *map_data, int num_elements,
+                                              void *extra_data) {
 
   struct engine *e = (struct engine *)extra_data;
   struct scheduler *sched = &e->sched;
@@ -1579,88 +1594,114 @@ void engine_link_stars_tasks_mapper(void *map_data, int num_elements,
   for (int ind = 0; ind < num_elements; ind++) {
     struct task *t = &((struct task *)map_data)[ind];
 
-    /* Sort tasks depend on the drift of the cell. */
+    /* Sort tasks depend on the drift and gravity drift of the cell. */
     if (t->type == task_type_stars_sort && t->ci->nodeID == engine_rank) {
+      scheduler_addunlock(sched, t->ci->hydro.super->hydro.drift, t);
       scheduler_addunlock(sched, t->ci->super->grav.drift, t);
     }
 
     /* Self-interaction? */
-    if (t->type == task_type_self && t->subtype == task_subtype_stars_density) {
-
-      /* Make the self-density tasks depend on the drifts. */
-      scheduler_addunlock(sched, t->ci->super->hydro.drift, t);
+    else if (t->type == task_type_self &&
+             t->subtype == task_subtype_stars_density) {
 
+      /* Make the self-density tasks depend on the drift and gravity drift. */
+      scheduler_addunlock(sched, t->ci->hydro.super->hydro.drift, t);
       scheduler_addunlock(sched, t->ci->super->grav.drift, t);
 
+      /* Start by constructing the task for the second stars loop */
+      struct task *t2 =
+          scheduler_addtask(sched, task_type_self, task_subtype_stars_feedback,
+                            0, 0, t->ci, NULL);
+
+      /* Add the link between the new loop and the cell */
+      engine_addlink(e, &t->ci->stars.feedback, t2);
+
       /* Now, build all the dependencies for the stars */
-      engine_make_stars_loops_dependencies(sched, t, t->ci);
-      if (t->ci == t->ci->super)
-        scheduler_addunlock(sched, t->ci->super->stars.ghost_out,
-                            t->ci->super->end_force);
+      engine_make_stars_loops_dependencies(sched, t, t2, t->ci);
+
+      /* end_force depends on feedback tasks */
+      scheduler_addunlock(sched, t2, t->ci->super->end_force);
     }
 
     /* Otherwise, pair interaction? */
     else if (t->type == task_type_pair &&
              t->subtype == task_subtype_stars_density) {
 
-      /* Make all density tasks depend on the drift and the sorts. */
-      if (t->cj->nodeID == engine_rank)
-        scheduler_addunlock(sched, t->cj->super->hydro.drift, t);
-      scheduler_addunlock(sched, t->cj->super->hydro.sorts, t);
-
+      /* Make all stars density tasks depend on the hydro drift and sorts,
+       * gravity drift and star sorts. */
+      if (t->ci->nodeID == engine_rank)
+        scheduler_addunlock(sched, t->ci->super->hydro.drift, t);
+      scheduler_addunlock(sched, t->ci->super->hydro.sorts, t);
       if (t->cj->nodeID == engine_rank)
         scheduler_addunlock(sched, t->cj->super->grav.drift, t);
       scheduler_addunlock(sched, t->ci->super->stars.sorts, t);
 
       if (t->ci->super != t->cj->super) {
-        if (t->ci->nodeID == engine_rank)
-          scheduler_addunlock(sched, t->ci->super->hydro.drift, t);
-        scheduler_addunlock(sched, t->ci->super->hydro.sorts, t);
-
+        if (t->cj->nodeID == engine_rank)
+          scheduler_addunlock(sched, t->cj->super->hydro.drift, t);
+        scheduler_addunlock(sched, t->cj->super->hydro.sorts, t);
         if (t->ci->nodeID == engine_rank)
           scheduler_addunlock(sched, t->ci->super->grav.drift, t);
         scheduler_addunlock(sched, t->cj->super->stars.sorts, t);
       }
 
+      /* Start by constructing the task for the second stars loop */
+      struct task *t2 =
+          scheduler_addtask(sched, task_type_pair, task_subtype_stars_feedback,
+                            0, 0, t->ci, t->cj);
+
+      /* Add the link between the new loop and both cells */
+      engine_addlink(e, &t->ci->stars.feedback, t2);
+      engine_addlink(e, &t->cj->stars.feedback, t2);
+
       /* Now, build all the dependencies for the stars for the cells */
-      /* that are local and are not descendant of the same super-cells */
       if (t->ci->nodeID == nodeID) {
-        engine_make_stars_loops_dependencies(sched, t, t->ci);
+        engine_make_stars_loops_dependencies(sched, t, t2, t->ci);
+        scheduler_addunlock(sched, t2, t->ci->super->end_force);
       }
       if (t->cj->nodeID == nodeID) {
         if (t->ci->super != t->cj->super)
-          engine_make_stars_loops_dependencies(sched, t, t->cj);
+          engine_make_stars_loops_dependencies(sched, t, t2, t->cj);
+        if (t->ci->super != t->cj->super)
+          scheduler_addunlock(sched, t2, t->cj->super->end_force);
       }
-
     }
 
     /* Otherwise, sub-self interaction? */
     else if (t->type == task_type_sub_self &&
              t->subtype == task_subtype_stars_density) {
 
-      /* Make all density tasks depend on the drift and sorts. */
+      /* Make all stars density tasks depend on the hydro drift and sorts,
+       * gravity drift and star sorts. */
       scheduler_addunlock(sched, t->ci->super->hydro.drift, t);
       scheduler_addunlock(sched, t->ci->super->hydro.sorts, t);
       scheduler_addunlock(sched, t->ci->super->grav.drift, t);
       scheduler_addunlock(sched, t->ci->super->stars.sorts, t);
 
+      /* Start by constructing the task for the second stars loop */
+      struct task *t2 = scheduler_addtask(sched, task_type_sub_self,
+                                          task_subtype_stars_feedback, t->flags,
+                                          0, t->ci, t->cj);
+
+      /* Add the link between the new loop and the cell */
+      engine_addlink(e, &t->ci->stars.feedback, t2);
+
       /* Now, build all the dependencies for the stars for the cells */
-      /* that are local and are not descendant of the same super-cells */
       if (t->ci->nodeID == nodeID) {
-        engine_make_stars_loops_dependencies(sched, t, t->ci);
-      } else
-        error("oo");
+        engine_make_stars_loops_dependencies(sched, t, t2, t->ci);
+        scheduler_addunlock(sched, t2, t->ci->super->end_force);
+      }
     }
 
     /* Otherwise, sub-pair interaction? */
     else if (t->type == task_type_sub_pair &&
              t->subtype == task_subtype_stars_density) {
 
-      /* Make all density tasks depend on the drift. */
+      /* Make all stars density tasks depend on the hydro drift and sorts,
+       * gravity drift and star sorts. */
       if (t->cj->nodeID == engine_rank)
         scheduler_addunlock(sched, t->cj->super->hydro.drift, t);
       scheduler_addunlock(sched, t->cj->super->hydro.sorts, t);
-
       if (t->cj->nodeID == engine_rank)
         scheduler_addunlock(sched, t->cj->super->grav.drift, t);
       scheduler_addunlock(sched, t->ci->super->stars.sorts, t);
@@ -1669,20 +1710,30 @@ void engine_link_stars_tasks_mapper(void *map_data, int num_elements,
         if (t->ci->nodeID == engine_rank)
           scheduler_addunlock(sched, t->ci->super->hydro.drift, t);
         scheduler_addunlock(sched, t->ci->super->hydro.sorts, t);
-
         if (t->ci->nodeID == engine_rank)
           scheduler_addunlock(sched, t->ci->super->grav.drift, t);
         scheduler_addunlock(sched, t->cj->super->stars.sorts, t);
       }
 
+      /* Start by constructing the task for the second stars loop */
+      struct task *t2 = scheduler_addtask(sched, task_type_sub_pair,
+                                          task_subtype_stars_feedback, t->flags,
+                                          0, t->ci, t->cj);
+
+      /* Add the link between the new loop and both cells */
+      engine_addlink(e, &t->ci->stars.feedback, t2);
+      engine_addlink(e, &t->cj->stars.feedback, t2);
+
       /* Now, build all the dependencies for the stars for the cells */
-      /* that are local and are not descendant of the same super-cells */
       if (t->ci->nodeID == nodeID) {
-        engine_make_stars_loops_dependencies(sched, t, t->ci);
+        engine_make_stars_loops_dependencies(sched, t, t2, t->ci);
+        scheduler_addunlock(sched, t2, t->ci->super->end_force);
       }
       if (t->cj->nodeID == nodeID) {
         if (t->ci->super != t->cj->super)
-          engine_make_stars_loops_dependencies(sched, t, t->cj);
+          engine_make_stars_loops_dependencies(sched, t, t2, t->cj);
+        if (t->ci->super != t->cj->super)
+          scheduler_addunlock(sched, t2, t->cj->super->end_force);
       }
     }
   }
@@ -1729,9 +1780,10 @@ void engine_make_starsloop_tasks_mapper(void *map_data, int num_elements,
     if (ci->stars.count == 0 && ci->hydro.count == 0) continue;
 
     /* If the cells is local build a self-interaction */
-    if (ci->nodeID == nodeID)
+    if (ci->nodeID == nodeID) {
       scheduler_addtask(sched, task_type_self, task_subtype_stars_density, 0, 0,
                         ci, NULL);
+    }
 
     /* Now loop over all the neighbours of this cell */
     for (int ii = -1; ii < 2; ii++) {
@@ -2104,24 +2156,29 @@ void engine_maketasks(struct engine *e) {
 
   tic2 = getticks();
 
-  /* Add the dependencies for the gravity stuff */
-  if (e->policy & (engine_policy_self_gravity | engine_policy_external_gravity))
-    engine_link_gravity_tasks(e);
+  /* Run through the tasks and make stars feedback tasks for each stars density
+     task. Each stars feedback task depends on the stars ghosts and unlocks the
+     kick task of its super-cell. */
+  if (e->policy & engine_policy_stars)
+    threadpool_map(&e->threadpool, engine_make_extra_starsloop_tasks_mapper,
+                   sched->tasks, sched->nr_tasks, sizeof(struct task), 0, e);
 
   if (e->verbose)
-    message("Linking gravity tasks took %.3f %s.",
+    message("Making extra starsloop tasks took %.3f %s.",
             clocks_from_ticks(getticks() - tic2), clocks_getunit());
 
   tic2 = getticks();
 
-  if (e->policy & engine_policy_stars)
-    threadpool_map(&e->threadpool, engine_link_stars_tasks_mapper, sched->tasks,
-                   sched->nr_tasks, sizeof(struct task), 0, e);
+  /* Add the dependencies for the gravity stuff */
+  if (e->policy & (engine_policy_self_gravity | engine_policy_external_gravity))
+    engine_link_gravity_tasks(e);
 
   if (e->verbose)
-    message("Linking stars tasks took %.3f %s.",
+    message("Linking gravity tasks took %.3f %s.",
             clocks_from_ticks(getticks() - tic2), clocks_getunit());
 
+  tic2 = getticks();
+
 #ifdef WITH_MPI
   if (e->policy & engine_policy_feedback)
     error("Cannot run stellar feedback with MPI (yet).");

src/engine_marktasks.c

@@ -141,6 +141,23 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
         }
       }
 
+      /* Activate the star feedback */
+      else if (t_type == task_type_self &&
+               t_subtype == task_subtype_stars_feedback) {
+        if (cell_is_active_stars(ci, e)) {
+          scheduler_activate(s, t);
+        }
+      }
+
+      /* Store current values of dx_max and h_max. */
+      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);
+          cell_activate_subcell_stars_tasks(ci, NULL, s);
+        }
+      }
+
       /* Activate the gravity drift */
       else if (t_type == task_type_self && t_subtype == task_subtype_grav) {
         if (cell_is_active_gravity(ci, e)) {
@@ -222,7 +239,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
         }
       }
 
-      /* Stars */
+      /* Stars density */
       if (t_subtype == task_subtype_stars_density &&
           ((ci_active_stars && ci->nodeID == engine_rank) ||
            (cj_active_stars && cj->nodeID == engine_rank))) {
@@ -279,6 +296,14 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
         }
       }
 
+      /* Stars feedback */
+      if (t_subtype == task_subtype_stars_feedback &&
+          ((ci_active_stars && ci->nodeID == engine_rank) ||
+           (cj_active_stars && cj->nodeID == engine_rank))) {
+
+        scheduler_activate(s, t);
+      }
+
       /* Gravity */
       if ((t_subtype == task_subtype_grav) &&
           ((ci_active_gravity && ci_nodeID == nodeID) ||
@@ -303,7 +328,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
       if (t_subtype == task_subtype_density) {
 
         /* Too much particle movement? */
-        if (cell_need_rebuild_for_pair(ci, cj)) *rebuild_space = 1;
+        if (cell_need_rebuild_for_hydro_pair(ci, cj)) *rebuild_space = 1;
 
 #ifdef WITH_MPI
         /* Activate the send/recv tasks. */
@@ -399,7 +424,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
       if (t->subtype == task_subtype_stars_density) {