diff --git a/examples/analyse_tasks.py b/examples/analyse_tasks.py
index 78d03ae5436a1fc56abcd6f6b99a4920011fa6db..48a00a63ea5da5cb80ebd6f10187cc613c1a5ed5 100755
--- a/examples/analyse_tasks.py
+++ b/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"]
diff --git a/examples/plot_tasks.py b/examples/plot_tasks.py
index a8d8c49a7ca8676278ac1ab50c7b6cdf3755445c..ec184ca1feee5879366d31ef8a984879bd0cddaf 100755
--- a/examples/plot_tasks.py
+++ b/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"]
diff --git a/src/cell.c b/src/cell.c
index c7d22639974d0c4ae639db62f8d082fa2d26452d..9a3a71a5a254159862d4efc10789031c46f1d9e9 100644
--- a/src/cell.c
+++ b/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);
diff --git a/src/cell.h b/src/cell.h
index 333f35ea6343a9dc7df3cf6fddc4bf24ba4888e3..1f5b308c729ac1f87d3e6fe81ed122ecf4842057 100644
--- a/src/cell.h
+++ b/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. */
diff --git a/src/engine.c b/src/engine.c
index 1942bf9e4f48fe51164a3b86ed0e42ead200695d..57cfaf1067230c8415eb787e8bada535be13fd27 100644
--- a/src/engine.c
+++ b/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);
diff --git a/src/engine.h b/src/engine.h
index 36b9177beabb23f33784032ca756b42c5bac073b..012021a1712fb6748ee44cec50b61c70615c0906 100644
--- a/src/engine.h
+++ b/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;
diff --git a/src/runner.c b/src/runner.c
index 73dc5bdfd4f1f8b7788b87a778609b448fbb2e70..0b01546918b0abee58402cc6541355239605f08e 100644
--- a/src/runner.c
+++ b/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;
diff --git a/src/task.c b/src/task.c
index 61d096cd32b6bab04c33dc439d5abe721bd7369f..efcc681e9c44d72885709581c74d960a8346f30c 100644
--- a/src/task.c
+++ b/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] = {
diff --git a/src/task.h b/src/task.h
index 866465e4997770d3e0e48818c6eebc6e9abf67a2..e38cfe87ea2ff3b157308fdb74f4e5c11a0c1256 100644
--- a/src/task.h
+++ b/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,
diff --git a/src/timers.c b/src/timers.c
index 51d0e5f6dc4fd9e4e0567592750b8de45ecda06b..88496af8bca17ea8605c964aaa277a5123ffc0c3 100644
--- a/src/timers.c
+++ b/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",
diff --git a/src/timers.h b/src/timers.h
index aba6ae33e3b8268c088694863967b96851715153..4eb14242df5f9f2adccc253a3485d2443f9fc8a4 100644
--- a/src/timers.h
+++ b/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,