Added star formation task.

examples/analyse_tasks.py

@@ -54,7 +54,7 @@ infile = args.input
 TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair",
              "init_grav", "init_grav_out", "ghost_in", "ghost", "ghost_out", "extra_ghost", "drift_part", "drift_gpart",
              "end_force", "kick1", "kick2", "timestep", "send", "recv", "grav_long_range", "grav_mm", "grav_down_in", 
-             "grav_down", "grav_mesh", "cooling", "sourceterms",
+             "grav_down", "grav_mesh", "cooling", "star_formation", "sourceterms",
              "stars_ghost_in", "stars_ghost",   "stars_ghost_out",
              "count"]
 

examples/plot_tasks.py

@@ -112,7 +112,7 @@ pl.rcParams.update(PLOT_PARAMS)
 TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair",
              "init_grav", "init_grav_out", "ghost_in", "ghost", "ghost_out", "extra_ghost", "drift_part", "drift_gpart",
              "end_force", "kick1", "kick2", "timestep", "send", "recv", "grav_long_range", "grav_mm", "grav_down_in", 
-             "grav_down", "grav_mesh", "cooling", "sourceterms",
+             "grav_down", "grav_mesh", "cooling", "star_formation", "sourceterms",
              "stars_ghost_in", "stars_ghost",   "stars_ghost_out",
              "count"]
 

src/cell.c

@@ -2380,6 +2380,8 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
     if (c->timestep != NULL) scheduler_activate(s, c->timestep);
     if (c->hydro.end_force != NULL) scheduler_activate(s, c->hydro.end_force);
     if (c->hydro.cooling != NULL) scheduler_activate(s, c->hydro.cooling);
+    if (c->hydro.star_formation != NULL)
+      scheduler_activate(s, c->hydro.star_formation);
     if (c->sourceterms != NULL) scheduler_activate(s, c->sourceterms);
   }
 
@@ -2527,6 +2529,9 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
     if (c->hydro.end_force != NULL) scheduler_activate(s, c->hydro.end_force);
     if ((e->policy & engine_policy_cooling) && c->hydro.cooling != NULL)
       scheduler_activate(s, c->hydro.cooling);
+    if ((e->policy & engine_policy_star_formation) &&
+        c->hydro.star_formation != NULL)
+      scheduler_activate(s, c->hydro.star_formation);
     if (c->grav.down != NULL) scheduler_activate(s, c->grav.down);
     if (c->grav.down_in != NULL) scheduler_activate(s, c->grav.down_in);
     if (c->grav.mesh != NULL) scheduler_activate(s, c->grav.mesh);

src/cell.h

@@ -326,6 +326,9 @@ struct cell {
     /*! Task for cooling */
     struct task *cooling;
 
+    /*! Task for star formation */
+    struct task *star_formation;
+
 #ifdef SWIFT_DEBUG_CHECKS
 
     /*! Last (integer) time the cell's sort arrays were updated. */

src/engine.c

@@ -214,6 +214,7 @@ void engine_make_hierarchical_tasks_common(struct engine *e, struct cell *c) {
 
   struct scheduler *s = &e->sched;
   const int is_with_cooling = (e->policy & engine_policy_cooling);
+  const int is_with_star_formation = (e->policy & engine_policy_star_formation);
 
   /* Are we in a super-cell ? */
   if (c->super == c) {
@@ -247,7 +248,18 @@ void engine_make_hierarchical_tasks_common(struct engine *e, struct cell *c) {
       } else {
         scheduler_addunlock(s, c->hydro.end_force, c->kick2);
       }
-      scheduler_addunlock(s, c->kick2, c->timestep);
+
+      if (is_with_star_formation) {
+
+        c->hydro.star_formation = scheduler_addtask(
+            s, task_type_star_formation, task_subtype_none, 0, 0, c, NULL);
+
+        scheduler_addunlock(s, c->kick2, c->hydro.star_formation);
+        scheduler_addunlock(s, c->hydro.star_formation, c->timestep);
+
+      } else {
+        scheduler_addunlock(s, c->kick2, c->timestep);
+      }
       scheduler_addunlock(s, c->timestep, c->kick1);
     }
 
@@ -4204,6 +4216,9 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
     else if (t_type == task_type_cooling) {
       if (cell_is_active_hydro(t->ci, e) || cell_is_active_gravity(t->ci, e))
         scheduler_activate(s, t);
+    } else if (t_type == task_type_star_formation) {
+      if (cell_is_active_hydro(t->ci, e) || cell_is_active_gravity(t->ci, e))
+        scheduler_activate(s, t);
     }
   }
 }
@@ -4419,6 +4434,12 @@ void engine_rebuild(struct engine *e, int clean_smoothing_length_values) {
   e->total_nr_gparts = num_particles[1];
   e->total_nr_sparts = num_particles[2];
 
+#ifdef SWIFT_DEBUG_CHECKS
+  e->count_inhibited_parts = 0;
+  e->count_inhibited_gparts = 0;
+  e->count_inhibited_sparts = 0;
+#endif
+
   /* Re-compute the mesh forces */
   if ((e->policy & engine_policy_self_gravity) && e->s->periodic)
     pm_mesh_compute_potential(e->mesh, e->s, &e->threadpool, e->verbose);

src/engine.h

@@ -73,9 +73,10 @@ enum engine_policy {
   engine_policy_sourceterms = (1 << 14),
   engine_policy_stars = (1 << 15),
   engine_policy_structure_finding = (1 << 16),
-  engine_policy_feedback = (1 << 17)
+  engine_policy_star_formation = (1 << 17),
+  engine_policy_feedback = (1 << 18)
 };
-#define engine_maxpolicy 17
+#define engine_maxpolicy 18
 extern const char *engine_policy_names[];
 
 /**
@@ -204,6 +205,12 @@ struct engine {
   /* Total numbers of particles in the system. */
   long long total_nr_parts, total_nr_gparts, total_nr_sparts;
 
+#ifdef SWIFT_DEBUG_CHECKS
+  /* Total number of particles removed from the system since the last rebuild */
+  long long count_inhibited_parts, count_inhibited_gparts,
+      count_inhibited_sparts;
+#endif
+
   /* Total mass in the simulation */
   double total_mass;
 

src/runner.c

@@ -486,6 +486,28 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
   if (timer) TIMER_TOC(timer_do_cooling);
 }
 
+/**
+ *
+ */
+void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
+
+  const struct engine *e = r->e;
+
+  TIMER_TIC;
+
+  /* Anything to do here? */
+  if (!cell_is_active_hydro(c, e)) return;
+
+  /* Recurse? */
+  if (c->split) {
+    for (int k = 0; k < 8; k++)
+      if (c->progeny[k] != NULL) runner_do_star_formation(r, c->progeny[k], 0);
+  } else {
+  }
+
+  if (timer) TIMER_TOC(timer_do_star_formation);
+}
+
 /**
  * @brief Sort the entries in ascending order using QuickSort.
  *
@@ -2013,7 +2035,7 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
   const int gcount = c->grav.count;
   const int scount = c->stars.count;
   struct part *restrict parts = c->hydro.parts;
-  // struct xpart *restrict xparts = c->hydro.xparts;
+  struct xpart *restrict xparts = c->hydro.xparts;
   struct gpart *restrict gparts = c->grav.parts;
   struct spart *restrict sparts = c->stars.parts;
   const int periodic = s->periodic;
@@ -2043,7 +2065,7 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
 
       /* Get a handle on the part. */
       struct part *restrict p = &parts[k];
-      // struct xpart *restrict xp = &xparts[k];
+      struct xpart *restrict xp = &xparts[k];
 
       if (part_is_active(p, e)) {
 
@@ -2052,9 +2074,9 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
 
         // MATTHIEU: Temporary star-formation law
         if (p->rho > 1.5e7 && e->step > 2) {
-          // message("Removing particle id=%lld rho=%e", p->id, p->rho);
+          message("Removing particle id=%lld rho=%e", p->id, p->rho);
           // cell_convert_part_to_gpart(e, c, p, xp);
-          // cell_remove_part(e,c,p,xp);
+          cell_remove_part(e, c, p, xp);
         }
       }
     }
@@ -2109,7 +2131,8 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
 
           /* Check that this gpart has interacted with all the other
            * particles (via direct or multipoles) in the box */
-          if (gp->num_interacted != e->total_nr_gparts) {
+          if (gp->num_interacted !=
+              e->total_nr_gparts - e->count_inhibited_gparts) {
 
             /* Get the ID of the gpart */
             long long my_id = 0;
@@ -2126,9 +2149,9 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
                 "g-particle (id=%lld, type=%s) did not interact "
                 "gravitationally with all other gparts "
                 "gp->num_interacted=%lld, total_gparts=%lld (local "
-                "num_gparts=%zd)",
+                "num_gparts=%zd inhibited_gparts=%lld)",
                 my_id, part_type_names[gp->type], gp->num_interacted,
-                e->total_nr_gparts, e->s->nr_gparts);
+                e->total_nr_gparts, e->s->nr_gparts, e->count_inhibited_gparts);
           }
         }
 #endif
@@ -2594,6 +2617,9 @@ void *runner_main(void *data) {
         case task_type_cooling:
           runner_do_cooling(r, t->ci, 1);
           break;
+        case task_type_star_formation:
+          runner_do_star_formation(r, t->ci, 1);
+          break;
         case task_type_sourceterms:
           runner_do_sourceterms(r, t->ci, 1);
           break;

src/task.c

@@ -48,16 +48,17 @@
 
 /* Task type names. */
 const char *taskID_names[task_type_count] = {
-    "none",           "sort",          "self",
-    "pair",           "sub_self",      "sub_pair",
-    "init_grav",      "init_grav_out", "ghost_in",
-    "ghost",          "ghost_out",     "extra_ghost",
-    "drift_part",     "drift_gpart",   "end_force",
-    "kick1",          "kick2",         "timestep",
-    "send",           "recv",          "grav_long_range",
-    "grav_mm",        "grav_down_in",  "grav_down",
-    "grav_mesh",      "cooling",       "sourceterms",
-    "stars_ghost_in", "stars_ghost",   "stars_ghost_out"};
+    "none",           "sort",           "self",
+    "pair",           "sub_self",       "sub_pair",
+    "init_grav",      "init_grav_out",  "ghost_in",
+    "ghost",          "ghost_out",      "extra_ghost",
+    "drift_part",     "drift_gpart",    "end_force",
+    "kick1",          "kick2",          "timestep",
+    "send",           "recv",           "grav_long_range",
+    "grav_mm",        "grav_down_in",   "grav_down",
+    "grav_mesh",      "cooling",        "star_formation",
+    "sourceterms",    "stars_ghost_in", "stars_ghost",
+    "stars_ghost_out"};
 
 /* Sub-task type names. */
 const char *subtaskID_names[task_subtype_count] = {

src/task.h

@@ -65,6 +65,7 @@ enum task_types {
   task_type_grav_down,
   task_type_grav_mesh,
   task_type_cooling,
+  task_type_star_formation,
   task_type_sourceterms,
   task_type_stars_ghost_in,
   task_type_stars_ghost,

src/timers.c

@@ -80,6 +80,7 @@ const char* timers_names[timer_count] = {
     "dorecv_gpart",
     "dorecv_spart",
     "do_cooling",
+    "do_star_formation",
     "gettask",
     "qget",
     "qsteal",

src/timers.h

@@ -81,6 +81,7 @@ enum {
   timer_dorecv_gpart,
   timer_dorecv_spart,
   timer_do_cooling,
+  timer_do_star_formation,
   timer_gettask,
   timer_qget,
   timer_qsteal,