diff --git a/README b/README
index 272188b3f7926b562ba993da3d24ae547c6a0397..b51abc121f7cc7c1b4baa851c02045b5f4614bbb 100644
--- a/README
+++ b/README
@@ -36,7 +36,8 @@ Parameters:
-s, --hydro Run with hydrodynamics.
-S, --stars Run with stars.
-x, --velociraptor Run with structure finding.
-
+ --limiter Run with time-step limiter.
+
Control options:
-a, --pin Pin runners using processor affinity.
diff --git a/README.md b/README.md
index 7a3c1287c79922a751595840295063a8ca347ef7..29415f27ee62f154b01dcd6a65414d7288a0a63f 100644
--- a/README.md
+++ b/README.md
@@ -84,6 +84,7 @@ Parameters:
-s, --hydro Run with hydrodynamics.
-S, --stars Run with stars.
-x, --velociraptor Run with structure finding.
+ --limiter Run with time-step limiter.
Control options:
diff --git a/doc/RTD/source/CommandLineOptions/index.rst b/doc/RTD/source/CommandLineOptions/index.rst
index bd58f031e622272d0245599621fc635891588a8f..e2603532b4ed4e64c86887f2a4f7c35f80cb08bf 100644
--- a/doc/RTD/source/CommandLineOptions/index.rst
+++ b/doc/RTD/source/CommandLineOptions/index.rst
@@ -31,6 +31,7 @@ can be found by typing ``./swift -h``::
-s, --hydro Run with hydrodynamics.
-S, --stars Run with stars.
-x, --velociraptor Run with structure finding.
+ --limiter Run with time-step limiter.
Control options:
diff --git a/examples/SedovBlast_1D/run.sh b/examples/SedovBlast_1D/run.sh
index ba479214961c5957a2b19d6aa118e0f0e7ee0f63..e5674dc15e8fac1b36f43da07b829720c0ecd5f1 100755
--- a/examples/SedovBlast_1D/run.sh
+++ b/examples/SedovBlast_1D/run.sh
@@ -8,7 +8,7 @@ then
fi
# Run SWIFT
-../swift --hydro --threads=1 sedov.yml 2>&1 | tee output.log
+../swift --hydro --limiter --threads=1 sedov.yml 2>&1 | tee output.log
# Plot the solution
python plotSolution.py 5
diff --git a/examples/SedovBlast_1D/sedov.yml b/examples/SedovBlast_1D/sedov.yml
index b4912a95e797440dc6eb0c9f48806a5954adbc41..16b5cf9bc5351931d0ae318a5b690c262eff176d 100644
--- a/examples/SedovBlast_1D/sedov.yml
+++ b/examples/SedovBlast_1D/sedov.yml
@@ -11,7 +11,7 @@ TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 5e-2 # The end time of the simulation (in internal units).
dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
- dt_max: 1e-5 # The maximal time-step size of the simulation (in internal units).
+ dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
@@ -21,7 +21,7 @@ Snapshots:
# Parameters governing the conserved quantities statistics
Statistics:
- delta_time: 1e-5 # Time between statistics output
+ delta_time: 1e-2 # Time between statistics output
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/examples/main.c b/examples/main.c
index db8d8f0d4137d67273233956f1d9c53988100148..aca3cf8b1f85c19d7522259f34094271ec9998c9 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -153,6 +153,7 @@ int main(int argc, char *argv[]) {
int with_stars = 0;
int with_star_formation = 0;
int with_feedback = 0;
+ int with_limiter = 0;
int with_fp_exceptions = 0;
int with_drift_all = 0;
int with_mpole_reconstruction = 0;
@@ -202,6 +203,8 @@ int main(int argc, char *argv[]) {
OPT_BOOLEAN('S', "stars", &with_stars, "Run with stars.", NULL, 0, 0),
OPT_BOOLEAN('x', "velociraptor", &with_structure_finding,
"Run with structure finding.", NULL, 0, 0),
+ OPT_BOOLEAN(0, "limiter", &with_limiter, "Run with time-step limiter.",
+ NULL, 0, 0),
OPT_GROUP(" Control options:\n"),
OPT_BOOLEAN('a', "pin", &with_aff,
@@ -456,6 +459,7 @@ int main(int argc, char *argv[]) {
if (with_feedback) error("Can't run with feedback over MPI (yet).");
if (with_star_formation)
error("Can't run with star formation over MPI (yet)");
+ if (with_limiter) error("Can't run with time-step limiter over MPI (yet)");
#endif
#if defined(WITH_MPI) && defined(HAVE_VELOCIRAPTOR)
@@ -878,6 +882,7 @@ int main(int argc, char *argv[]) {
engine_policies |= engine_policy_external_gravity;
if (with_cosmology) engine_policies |= engine_policy_cosmology;
if (with_temperature) engine_policies |= engine_policy_temperature;
+ if (with_limiter) engine_policies |= engine_policy_limiter;
if (with_cooling) engine_policies |= engine_policy_cooling;
if (with_stars) engine_policies |= engine_policy_stars;
if (with_star_formation) engine_policies |= engine_policy_star_formation;
diff --git a/src/Makefile.am b/src/Makefile.am
index d9b629fca5cb7227d1e183581744fefa9f1728dc..db8cc00ed7dbb7d0b3b6812f06ce3daab8b9920e 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -75,7 +75,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
gravity_iact.h kernel_long_gravity.h vector.h cache.h runner_doiact.h runner_doiact_vec.h runner_doiact_grav.h \
runner_doiact_nosort.h runner_doiact_stars.h units.h intrinsics.h minmax.h kick.h timestep.h drift.h \
adiabatic_index.h io_properties.h dimension.h part_type.h periodic.h memswap.h dump.h logger.h sign.h \
- logger_io.h \
+ logger_io.h timestep_limiter.h \
gravity.h gravity_io.h gravity_cache.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
diff --git a/src/cell.c b/src/cell.c
index 10423e881f6140b6407e88deb543cf5b9ba944a9..e3ac72519ffe457449dc6b7751a861c32738d12b 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1232,8 +1232,11 @@ void cell_clean_links(struct cell *c, void *data) {
c->hydro.density = NULL;
c->hydro.gradient = NULL;
c->hydro.force = NULL;
+ c->hydro.limiter = NULL;
c->grav.grav = NULL;
c->grav.mm = NULL;
+ c->stars.density = NULL;
+ c->stars.feedback = NULL;
}
/**
@@ -1599,6 +1602,14 @@ void cell_clear_drift_flags(struct cell *c, void *data) {
c->grav.do_sub_drift = 0;
}
+/**
+ * @brief Clear the limiter flags on the given cell.
+ */
+void cell_clear_limiter_flags(struct cell *c, void *data) {
+ c->hydro.do_limiter = 0;
+ c->hydro.do_sub_limiter = 0;
+}
+
/**
* @brief Activate the #part drifts on the given cell.
*/
@@ -1686,6 +1697,34 @@ void cell_activate_drift_spart(struct cell *c, struct scheduler *s) {
cell_activate_drift_gpart(c, s);
}
+/**
+ * @brief Activate the drifts on the given cell.
+ */
+void cell_activate_limiter(struct cell *c, struct scheduler *s) {
+
+ /* If this cell is already marked for drift, quit early. */
+ if (c->hydro.do_limiter) return;
+
+ /* Mark this cell for drifting. */
+ c->hydro.do_limiter = 1;
+
+ /* Set the do_sub_limiter all the way up and activate the super limiter
+ if this has not yet been done. */
+ if (c == c->super) {
+ scheduler_activate(s, c->timestep_limiter);
+ } else {
+ for (struct cell *parent = c->parent;
+ parent != NULL && !parent->hydro.do_sub_limiter;
+ parent = parent->parent) {
+ parent->hydro.do_sub_limiter = 1;
+ if (parent == c->super) {
+ scheduler_activate(s, parent->timestep_limiter);
+ break;
+ }
+ }
+ }
+}
+
/**
* @brief Activate the sorts up a cell hierarchy.
*/
@@ -1816,6 +1855,7 @@ void cell_activate_stars_sorts(struct cell *c, int sid, struct scheduler *s) {
void cell_activate_subcell_hydro_tasks(struct cell *ci, struct cell *cj,
struct scheduler *s) {
const struct engine *e = s->space->e;
+ const int with_limiter = (e->policy & engine_policy_limiter);
/* Store the current dx_max and h_max values. */
ci->hydro.dx_max_part_old = ci->hydro.dx_max_part;
@@ -1849,6 +1889,7 @@ void cell_activate_subcell_hydro_tasks(struct cell *ci, struct cell *cj,
/* We have reached the bottom of the tree: activate drift */
cell_activate_drift_part(ci, s);
+ if (with_limiter) cell_activate_limiter(ci, s);
}
}
@@ -2154,6 +2195,12 @@ void cell_activate_subcell_hydro_tasks(struct cell *ci, struct cell *cj,
if (ci->nodeID == engine_rank) cell_activate_drift_part(ci, s);
if (cj->nodeID == engine_rank) cell_activate_drift_part(cj, s);
+ /* Also activate the time-step limiter */
+ if (ci->nodeID == engine_rank && with_limiter)
+ cell_activate_limiter(ci, s);
+ if (cj->nodeID == engine_rank && with_limiter)
+ cell_activate_limiter(cj, s);
+
/* Do we need to sort the cells? */
cell_activate_hydro_sorts(ci, sid, s);
cell_activate_hydro_sorts(cj, sid, s);
@@ -2718,6 +2765,7 @@ 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);
int rebuild = 0;
/* Un-skip the density tasks involved with this cell. */
@@ -2743,6 +2791,7 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
/* Activate hydro drift */
if (t->type == task_type_self) {
if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
+ if (ci->nodeID == nodeID && with_limiter) cell_activate_limiter(ci, s);
}
/* Set the correct sorting flags and activate hydro drifts */
@@ -2757,6 +2806,10 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
+ /* Activate the limiter tasks. */
+ if (ci->nodeID == nodeID && with_limiter) cell_activate_limiter(ci, s);
+ if (cj->nodeID == nodeID && with_limiter) cell_activate_limiter(cj, s);
+
/* Check the sorts and activate them if needed. */
cell_activate_hydro_sorts(ci, t->flags, s);
cell_activate_hydro_sorts(cj, t->flags, s);
@@ -2791,7 +2844,11 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
}
/* If the foreign cell is active, we want its ti_end values. */
- if (ci_active) scheduler_activate(s, ci->mpi.recv_ti);
+ if (ci_active || with_limiter) scheduler_activate(s, ci->mpi.recv_ti);
+
+ if (with_limiter) scheduler_activate(s, ci->mpi.limiter.recv);
+ if (with_limiter)
+ scheduler_activate_send(s, cj->mpi.limiter.send, ci->nodeID);
/* Is the foreign cell active and will need stuff from us? */
if (ci_active) {
@@ -2801,6 +2858,7 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
/* Drift the cell which will be sent; note that not all sent
particles will be drifted, only those that are needed. */
cell_activate_drift_part(cj, s);
+ if (with_limiter) cell_activate_limiter(cj, s);
/* If the local cell is also active, more stuff will be needed. */
if (cj_active) {
@@ -2813,7 +2871,8 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
}
/* If the local cell is active, send its ti_end values. */
- if (cj_active) scheduler_activate_send(s, cj->mpi.send_ti, ci_nodeID);
+ if (cj_active || with_limiter)
+ scheduler_activate_send(s, cj->mpi.send_ti, ci_nodeID);
} else if (cj_nodeID != nodeID) {
@@ -2830,7 +2889,11 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
}
/* If the foreign cell is active, we want its ti_end values. */
- if (cj_active) scheduler_activate(s, cj->mpi.recv_ti);
+ if (cj_active || with_limiter) scheduler_activate(s, cj->mpi.recv_ti);
+
+ if (with_limiter) scheduler_activate(s, cj->mpi.limiter.recv);
+ if (with_limiter)
+ scheduler_activate_send(s, ci->mpi.limiter.send, cj->nodeID);
/* Is the foreign cell active and will need stuff from us? */
if (cj_active) {
@@ -2840,6 +2903,7 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
/* Drift the cell which will be sent; note that not all sent
particles will be drifted, only those that are needed. */
cell_activate_drift_part(ci, s);
+ if (with_limiter) cell_activate_limiter(ci, s);
/* If the local cell is also active, more stuff will be needed. */
if (ci_active) {
@@ -2853,7 +2917,8 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
}
/* If the local cell is active, send its ti_end values. */
- if (ci_active) scheduler_activate_send(s, ci->mpi.send_ti, cj_nodeID);
+ if (ci_active || with_limiter)
+ scheduler_activate_send(s, ci->mpi.send_ti, cj_nodeID);
}
#endif
}
@@ -2866,6 +2931,8 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
scheduler_activate(s, l->t);
for (struct link *l = c->hydro.force; l != NULL; l = l->next)
scheduler_activate(s, l->t);
+ for (struct link *l = c->hydro.limiter; l != NULL; l = l->next)
+ scheduler_activate(s, l->t);
if (c->hydro.extra_ghost != NULL)
scheduler_activate(s, c->hydro.extra_ghost);
@@ -2879,7 +2946,6 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
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);
if (c->logger != NULL) scheduler_activate(s, c->logger);
}
diff --git a/src/cell.h b/src/cell.h
index 9452accbab6312f235764e689522ac3edbaafe61..8b00be8165a359f43bf13dc370d530f5b9046dc5 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -263,6 +263,9 @@ struct cell {
/*! Linked list of the tasks computing this cell's hydro forces. */
struct link *force;
+ /*! Linked list of the tasks computing this cell's limiter. */
+ struct link *limiter;
+
/*! Dependency implicit task for the ghost (in->ghost->out)*/
struct task *ghost_in;
@@ -348,6 +351,12 @@ struct cell {
/*! Do any of this cell's sub-cells need to be sorted? */
char do_sub_sort;
+ /*! Does this cell need to be limited? */
+ char do_limiter;
+
+ /*! Do any of this cell's sub-cells need to be limited? */
+ char do_sub_limiter;
+
#ifdef SWIFT_DEBUG_CHECKS
/*! Last (integer) time the cell's sort arrays were updated. */
@@ -570,6 +579,15 @@ struct cell {
struct link *send;
} grav;
+ struct {
+
+ /* Task receiving gpart data. */
+ struct task *recv;
+
+ /* Linked list for sending gpart data. */
+ struct link *send;
+ } limiter;
+
/* Task receiving data (time-step). */
struct task *recv_ti;
@@ -603,8 +621,8 @@ struct cell {
/*! The task to compute time-steps */
struct task *timestep;
- /*! Task for source terms */
- struct task *sourceterms;
+ /*! The task to limit the time-step of inactive particles */
+ struct task *timestep_limiter;
/*! The logger task */
struct task *logger;
@@ -705,7 +723,9 @@ void cell_activate_drift_gpart(struct cell *c, struct scheduler *s);
void cell_activate_drift_spart(struct cell *c, struct scheduler *s);
void cell_activate_hydro_sorts(struct cell *c, int sid, struct scheduler *s);
void cell_activate_stars_sorts(struct cell *c, int sid, struct scheduler *s);
+void cell_activate_limiter(struct cell *c, struct scheduler *s);
void cell_clear_drift_flags(struct cell *c, void *data);
+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);
diff --git a/src/const.h b/src/const.h
index e417b8ca3827ef87396706c56df36bb9bd3aed75..613a48920e6f26c209faf6e354b82c2ed5be0bf1 100644
--- a/src/const.h
+++ b/src/const.h
@@ -33,6 +33,9 @@
/* Time integration constants. */
#define const_max_u_change 0.1f
+/* Time-step limiter maximal difference in signal velocity */
+#define const_limiter_max_v_sig_ratio 4.1f
+
/* Type of gradients to use (GIZMO_SPH only) */
/* If no option is chosen, no gradients are used (first order scheme) */
//#define GRADIENTS_SPH
diff --git a/src/engine.c b/src/engine.c
index ce9153e8c856330b2d60c97a01b36f5c7af670a4..2734fe6a1d572f19c6763931c314535170cc8b8c 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -113,7 +113,8 @@ const char *engine_policy_names[] = {"none",
"stars",
"structure finding",
"star formation",
- "feedback"};
+ "feedback",
+ "time-step limiter"};
/** The rank of the engine as a global variable (for messages). */
int engine_rank;
@@ -1907,6 +1908,10 @@ int engine_estimate_nr_tasks(struct engine *e) {
#endif
#endif
}
+ if (e->policy & engine_policy_limiter) {
+ n1 += 18;
+ n2 += 1;
+ }
if (e->policy & engine_policy_self_gravity) {
n1 += 125;
n2 += 8;
@@ -2585,8 +2590,11 @@ void engine_skip_force_and_kick(struct engine *e) {
/* Skip everything that updates the particles */
if (t->type == task_type_drift_part || t->type == task_type_drift_gpart ||
t->type == task_type_kick1 || t->type == task_type_kick2 ||
- t->type == task_type_timestep || t->subtype == task_subtype_force ||
- t->subtype == task_subtype_grav || t->type == task_type_end_force ||
+ t->type == task_type_timestep ||
+ t->type == task_type_timestep_limiter ||
+ t->subtype == task_subtype_force ||
+ t->subtype == task_subtype_limiter || t->subtype == task_subtype_grav ||
+ t->type == task_type_end_force ||
t->type == task_type_grav_long_range || t->type == task_type_grav_mm ||
t->type == task_type_grav_down || t->type == task_type_cooling)
t->skip = 1;
@@ -2594,6 +2602,7 @@ void engine_skip_force_and_kick(struct engine *e) {
/* Run through the cells and clear some flags. */
space_map_cells_pre(e->s, 1, cell_clear_drift_flags, NULL);
+ space_map_cells_pre(e->s, 1, cell_clear_limiter_flags, NULL);
}
/**
@@ -2803,6 +2812,11 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
gravity_exact_force_check(e->s, e, 1e-1);
#endif
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Make sure all woken-up particles have been processed */
+ space_check_limiter(e->s);
+#endif
+
/* Recover the (integer) end of the next time-step */
engine_collect_end_of_step(e, 1);
@@ -3060,6 +3074,11 @@ void engine_step(struct engine *e) {
gravity_exact_force_check(e->s, e, 1e-1);
#endif
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Make sure all woken-up particles have been processed */
+ space_check_limiter(e->s);
+#endif
+
/* Collect information about the next time-step */
engine_collect_end_of_step(e, 1);
e->forcerebuild = e->collect_group1.forcerebuild;
diff --git a/src/engine.h b/src/engine.h
index fea874435cd8160a4dfd37f6e7530f36318646d7..ed73ff6b0bb4e8ac0eb2a8a5ad63f5e3427e4eeb 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -74,9 +74,10 @@ enum engine_policy {
engine_policy_stars = (1 << 15),
engine_policy_structure_finding = (1 << 16),
engine_policy_star_formation = (1 << 17),
- engine_policy_feedback = (1 << 18)
+ engine_policy_feedback = (1 << 18),
+ engine_policy_limiter = (1 << 19)
};
-#define engine_maxpolicy 19
+#define engine_maxpolicy 20
extern const char *engine_policy_names[engine_maxpolicy + 1];
/**
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index f280c2a6106d9da59e58f79d9acb9e8e199040b4..1c9ae1caf4966c9551ba22d3bede72c81e5ebe06 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -210,9 +210,13 @@ 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_ti The send_ti #task, if it has already been created.
+ * @param t_limiter The send_limiter #task, if already created.
+ * @param with_limiter Are we running with the time-step limiter?
*/
void engine_addtasks_send_timestep(struct engine *e, struct cell *ci,
- struct cell *cj, struct task *t_ti) {
+ struct cell *cj, struct task *t_ti,
+ struct task *t_limiter,
+ const int with_limiter) {
#ifdef WITH_MPI
struct link *l = NULL;
@@ -244,19 +248,31 @@ void engine_addtasks_send_timestep(struct engine *e, struct cell *ci,
t_ti = scheduler_addtask(s, task_type_send, task_subtype_tend,
ci->mpi.tag, 0, ci, cj);
+ if (with_limiter)
+ t_limiter = scheduler_addtask(s, task_type_send, task_subtype_limiter,
+ ci->mpi.tag, 0, ci, cj);
+
/* The super-cell's timestep task should unlock the send_ti task. */
scheduler_addunlock(s, ci->super->timestep, t_ti);
+ if (with_limiter) scheduler_addunlock(s, t_limiter, ci->super->timestep);
+ if (with_limiter)
+ scheduler_addunlock(s, t_limiter, ci->super->timestep_limiter);
+ if (with_limiter) scheduler_addunlock(s, ci->super->kick2, t_limiter);
+ if (with_limiter)
+ scheduler_addunlock(s, ci->super->timestep_limiter, t_ti);
}
/* Add them to the local cell. */
engine_addlink(e, &ci->mpi.send_ti, t_ti);
+ if (with_limiter) engine_addlink(e, &ci->mpi.limiter.send, t_limiter);
}
/* Recurse? */
if (ci->split)
for (int k = 0; k < 8; k++)
if (ci->progeny[k] != NULL)
- engine_addtasks_send_timestep(e, ci->progeny[k], cj, t_ti);
+ engine_addtasks_send_timestep(e, ci->progeny[k], cj, t_ti, t_limiter,
+ with_limiter);
#else
error("SWIFT was not compiled with MPI support.");
@@ -380,9 +396,12 @@ void engine_addtasks_recv_gravity(struct engine *e, struct cell *c,
* @param e The #engine.
* @param c The foreign #cell.
* @param t_ti The recv_ti #task, if already been created.
+ * @param t_limiter The recv_limiter #task, if already created.
+ * @param with_limiter Are we running with the time-step limiter?
*/
void engine_addtasks_recv_timestep(struct engine *e, struct cell *c,
- struct task *t_ti) {
+ struct task *t_ti, struct task *t_limiter,
+ const int with_limiter) {
#ifdef WITH_MPI
struct scheduler *s = &e->sched;
@@ -397,21 +416,42 @@ void engine_addtasks_recv_timestep(struct engine *e, struct cell *c,
t_ti = scheduler_addtask(s, task_type_recv, task_subtype_tend, c->mpi.tag,
0, c, NULL);
+
+ if (with_limiter)
+ t_limiter = scheduler_addtask(s, task_type_recv, task_subtype_limiter,
+ c->mpi.tag, 0, c, NULL);
}
c->mpi.recv_ti = t_ti;
- for (struct link *l = c->grav.grav; l != NULL; l = l->next)
+ for (struct link *l = c->grav.grav; l != NULL; l = l->next) {
scheduler_addunlock(s, l->t, t_ti);
+ }
- for (struct link *l = c->hydro.force; l != NULL; l = l->next)
- scheduler_addunlock(s, l->t, t_ti);
+ if (with_limiter) {
+
+ for (struct link *l = c->hydro.force; l != NULL; l = l->next) {
+ scheduler_addunlock(s, l->t, t_limiter);
+ }
+
+ for (struct link *l = c->hydro.limiter; l != NULL; l = l->next) {
+ scheduler_addunlock(s, t_limiter, l->t);
+ scheduler_addunlock(s, l->t, t_ti);
+ }
+
+ } else {
+
+ for (struct link *l = c->hydro.force; l != NULL; l = l->next) {
+ scheduler_addunlock(s, l->t, t_ti);
+ }
+ }
/* Recurse? */
if (c->split)
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL)
- engine_addtasks_recv_timestep(e, c->progeny[k], t_ti);
+ engine_addtasks_recv_timestep(e, c->progeny[k], t_ti, t_limiter,
+ with_limiter);
#else
error("SWIFT was not compiled with MPI support.");
@@ -435,6 +475,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);
+ const int with_limiter = (e->policy & engine_policy_limiter);
/* Are we in a super-cell ? */
if (c->super == c) {
@@ -489,6 +530,18 @@ void engine_make_hierarchical_tasks_common(struct engine *e, struct cell *c) {
scheduler_addunlock(s, c->timestep, c->kick1);
+ /* Time-step limiting */
+ if (with_limiter) {
+ c->timestep_limiter = scheduler_addtask(
+ s, task_type_timestep_limiter, task_subtype_none, 0, 0, c, NULL);
+
+ /* Make sure it is not run before kick2 */
+ scheduler_addunlock(s, c->timestep, c->timestep_limiter);
+ scheduler_addunlock(s, c->timestep_limiter, c->kick1);
+ } else {
+ scheduler_addunlock(s, c->kick2, c->timestep);
+ }
+
#if defined(WITH_LOGGER)
scheduler_addunlock(s, c->kick1, c->logger);
#endif
@@ -1274,7 +1327,8 @@ void engine_link_gravity_tasks(struct engine *e) {
*/
static inline void engine_make_hydro_loops_dependencies(
struct scheduler *sched, struct task *density, struct task *gradient,
- struct task *force, struct cell *c, int with_cooling) {
+ struct task *force, struct task *limiter, struct cell *c, int with_cooling,
+ int with_limiter) {
/* density loop --> ghost --> gradient loop --> extra_ghost */
/* extra_ghost --> force loop */
@@ -1292,17 +1346,24 @@ static inline void engine_make_hydro_loops_dependencies(
* @param sched The #scheduler.
* @param density The density task to link.
* @param force The force task to link.
+ * @param limiter The limiter task to link.
* @param c The cell.
- * @param with_cooling Are we running with cooling switched on ?
+ * @param with_cooling Are we running with cooling switched on?
+ * @param with_limiter Are we running with limiter switched on?
*/
-static inline void engine_make_hydro_loops_dependencies(struct scheduler *sched,
- struct task *density,
- struct task *force,
- struct cell *c,
- int with_cooling) {
+static inline void engine_make_hydro_loops_dependencies(
+ struct scheduler *sched, struct task *density, struct task *force,
+ struct task *limiter, struct cell *c, int with_cooling, int with_limiter) {
/* density loop --> ghost --> force loop */
scheduler_addunlock(sched, density, c->hydro.super->hydro.ghost_in);
scheduler_addunlock(sched, c->hydro.super->hydro.ghost_out, force);
+
+ if (with_limiter) {
+ scheduler_addunlock(sched, c->super->kick2, limiter);
+ scheduler_addunlock(sched, limiter, c->super->timestep);
+ if (limiter->type != task_type_self)
+ scheduler_addunlock(sched, limiter, c->super->timestep_limiter);
+ }
}
#endif
@@ -1340,6 +1401,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
struct scheduler *sched = &e->sched;
const int nodeID = e->nodeID;
const int with_cooling = (e->policy & engine_policy_cooling);
+ const int with_limiter = (e->policy & engine_policy_limiter);
+#ifdef EXTRA_HYDRO_LOOP
+ struct task *t_gradient = NULL;
+#endif
+ struct task *t_force = NULL;
+ struct task *t_limiter = NULL;
for (int ind = 0; ind < num_elements; ind++) {
struct task *t = &((struct task *)map_data)[ind];
@@ -1357,31 +1424,45 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
#ifdef EXTRA_HYDRO_LOOP
/* Start by constructing the task for the second and third hydro loop. */
- struct task *t2 = scheduler_addtask(
- sched, task_type_self, task_subtype_gradient, 0, 0, t->ci, NULL);
- struct task *t3 = scheduler_addtask(
- sched, task_type_self, task_subtype_force, 0, 0, t->ci, NULL);
+ t_gradient = scheduler_addtask(sched, task_type_self,
+ task_subtype_gradient, 0, 0, t->ci, NULL);
+ t_force = scheduler_addtask(sched, task_type_self, task_subtype_force, 0,
+ 0, t->ci, NULL);
+
+ /* and the task for the time-step limiter */
+ if (with_limiter)
+ t_limiter = scheduler_addtask(sched, task_type_self,
+ task_subtype_limiter, 0, 0, t->ci, NULL);
/* Add the link between the new loops and the cell */
- engine_addlink(e, &t->ci->hydro.gradient, t2);
- engine_addlink(e, &t->ci->hydro.force, t3);
+ engine_addlink(e, &t->ci->hydro.gradient, t_gradient);
+ engine_addlink(e, &t->ci->hydro.force, t_force);
+ if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
/* Now, build all the dependencies for the hydro */
- engine_make_hydro_loops_dependencies(sched, t, t2, t3, t->ci,
- with_cooling);
- scheduler_addunlock(sched, t3, t->ci->super->end_force);
+ engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
+ t_limiter, t->ci, with_cooling,
+ with_limiter);
+ scheduler_addunlock(sched, t_force, t->ci->super->end_force);
#else
/* Start by constructing the task for the second hydro loop */
- struct task *t2 = scheduler_addtask(
- sched, task_type_self, task_subtype_force, 0, 0, t->ci, NULL);
+ t_force = scheduler_addtask(sched, task_type_self, task_subtype_force, 0,
+ 0, t->ci, NULL);
+
+ /* and the task for the time-step limiter */
+ if (with_limiter)
+ t_limiter = scheduler_addtask(sched, task_type_self,
+ task_subtype_limiter, 0, 0, t->ci, NULL);
/* Add the link between the new loop and the cell */
- engine_addlink(e, &t->ci->hydro.force, t2);
+ engine_addlink(e, &t->ci->hydro.force, t_force);
+ if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
/* Now, build all the dependencies for the hydro */
- engine_make_hydro_loops_dependencies(sched, t, t2, t->ci, with_cooling);
- scheduler_addunlock(sched, t2, t->ci->super->end_force);
+ engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter, t->ci,
+ with_cooling, with_limiter);
+ scheduler_addunlock(sched, t_force, t->ci->super->end_force);
#endif
}
@@ -1400,54 +1481,71 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
#ifdef EXTRA_HYDRO_LOOP
/* Start by constructing the task for the second and third hydro loop */
- struct task *t2 = scheduler_addtask(
- sched, task_type_pair, task_subtype_gradient, 0, 0, t->ci, t->cj);
- struct task *t3 = scheduler_addtask(
- sched, task_type_pair, task_subtype_force, 0, 0, t->ci, t->cj);
+ t_gradient = scheduler_addtask(sched, task_type_pair,
+ task_subtype_gradient, 0, 0, t->ci, t->cj);
+ t_force = scheduler_addtask(sched, task_type_pair, task_subtype_force, 0,
+ 0, t->ci, t->cj);
+
+ /* and the task for the time-step limiter */
+ if (with_limiter)
+ t_limiter = scheduler_addtask(sched, task_type_pair,
+ task_subtype_limiter, 0, 0, t->ci, t->cj);
/* Add the link between the new loop and both cells */
- engine_addlink(e, &t->ci->hydro.gradient, t2);
- engine_addlink(e, &t->cj->hydro.gradient, t2);
- engine_addlink(e, &t->ci->hydro.force, t3);
- engine_addlink(e, &t->cj->hydro.force, t3);
+ engine_addlink(e, &t->ci->hydro.gradient, t_gradient);
+ engine_addlink(e, &t->cj->hydro.gradient, t_gradient);
+ engine_addlink(e, &t->ci->hydro.force, t_force);
+ engine_addlink(e, &t->cj->hydro.force, t_force);
+ if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
+ if (with_limiter) engine_addlink(e, &t->cj->hydro.limiter, t_limiter);
/* Now, build all the dependencies for the hydro for the cells */
/* that are local and are not descendant of the same super_hydro-cells */
if (t->ci->nodeID == nodeID) {
- engine_make_hydro_loops_dependencies(sched, t, t2, t3, t->ci,
- with_cooling);
- scheduler_addunlock(sched, t3, t->ci->super->end_force);
+ engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
+ t_limiter, t->ci, with_cooling,
+ with_limiter);
+ scheduler_addunlock(sched, t_force, t->ci->super->end_force);
}
if (t->cj->nodeID == nodeID) {
if (t->ci->hydro.super != t->cj->hydro.super)
- engine_make_hydro_loops_dependencies(sched, t, t2, t3, t->cj,
- with_cooling);
+ engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
+ t_limiter, t->cj, with_cooling,
+ with_limiter);
if (t->ci->super != t->cj->super)
- scheduler_addunlock(sched, t3, t->cj->super->end_force);
+ scheduler_addunlock(sched, t_force, t->cj->super->end_force);
}
#else
/* Start by constructing the task for the second hydro loop */
- struct task *t2 = scheduler_addtask(
- sched, task_type_pair, task_subtype_force, 0, 0, t->ci, t->cj);
+ t_force = scheduler_addtask(sched, task_type_pair, task_subtype_force, 0,
+ 0, t->ci, t->cj);
+
+ /* and the task for the time-step limiter */
+ if (with_limiter)
+ t_limiter = scheduler_addtask(sched, task_type_pair,
+ task_subtype_limiter, 0, 0, t->ci, t->cj);
/* Add the link between the new loop and both cells */
- engine_addlink(e, &t->ci->hydro.force, t2);
- engine_addlink(e, &t->cj->hydro.force, t2);
+ engine_addlink(e, &t->ci->hydro.force, t_force);
+ engine_addlink(e, &t->cj->hydro.force, t_force);
+ if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
+ if (with_limiter) engine_addlink(e, &t->cj->hydro.limiter, t_limiter);
/* Now, build all the dependencies for the hydro for the cells */
/* that are local and are not descendant of the same super_hydro-cells */
if (t->ci->nodeID == nodeID) {
- engine_make_hydro_loops_dependencies(sched, t, t2, t->ci, with_cooling);
- scheduler_addunlock(sched, t2, t->ci->super->end_force);
+ engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter,
+ t->ci, with_cooling, with_limiter);
+ scheduler_addunlock(sched, t_force, t->ci->super->end_force);
}
if (t->cj->nodeID == nodeID) {
if (t->ci->hydro.super != t->cj->hydro.super)
- engine_make_hydro_loops_dependencies(sched, t, t2, t->cj,
- with_cooling);
+ engine_make_hydro_loops_dependencies(
+ sched, t, t_force, t_limiter, t->cj, with_cooling, with_limiter);
if (t->ci->super != t->cj->super)
- scheduler_addunlock(sched, t2, t->cj->super->end_force);
+ scheduler_addunlock(sched, t_force, t->cj->super->end_force);
}
#endif
@@ -1465,39 +1563,53 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
#ifdef EXTRA_HYDRO_LOOP
/* Start by constructing the task for the second and third hydro loop */
- struct task *t2 =
+ t_gradient =
scheduler_addtask(sched, task_type_sub_self, task_subtype_gradient,
- t->flags, 0, t->ci, t->cj);
- struct task *t3 =
- scheduler_addtask(sched, task_type_sub_self, task_subtype_force,
- t->flags, 0, t->ci, t->cj);
+ t->flags, 0, t->ci, NULL);
+ t_force = scheduler_addtask(sched, task_type_sub_self, task_subtype_force,
+ t->flags, 0, t->ci, NULL);
+
+ /* and the task for the time-step limiter */
+ if (with_limiter)
+ t_limiter =
+ scheduler_addtask(sched, task_type_sub_self, task_subtype_limiter,
+ t->flags, 0, t->ci, NULL);
/* Add the link between the new loop and the cell */
- engine_addlink(e, &t->ci->hydro.gradient, t2);
- engine_addlink(e, &t->ci->hydro.force, t3);
+ engine_addlink(e, &t->ci->hydro.gradient, t_gradient);
+ engine_addlink(e, &t->ci->hydro.force, t_force);
+ if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
/* Now, build all the dependencies for the hydro for the cells */
/* that are local and are not descendant of the same super_hydro-cells */
if (t->ci->nodeID == nodeID) {
- engine_make_hydro_loops_dependencies(sched, t, t2, t3, t->ci,
- with_cooling);
- scheduler_addunlock(sched, t3, t->ci->super->end_force);
+ engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
+ t_limiter, t->ci, with_cooling,
+ with_limiter);
+ scheduler_addunlock(sched, t_force, t->ci->super->end_force);
}
#else
/* Start by constructing the task for the second hydro loop */
- struct task *t2 =
- scheduler_addtask(sched, task_type_sub_self, task_subtype_force,
- t->flags, 0, t->ci, t->cj);
+ t_force = scheduler_addtask(sched, task_type_sub_self, task_subtype_force,
+ t->flags, 0, t->ci, NULL);
+
+ /* and the task for the time-step limiter */
+ if (with_limiter)
+ t_limiter =
+ scheduler_addtask(sched, task_type_sub_self, task_subtype_limiter,
+ t->flags, 0, t->ci, NULL);
/* Add the link between the new loop and the cell */
- engine_addlink(e, &t->ci->hydro.force, t2);
+ engine_addlink(e, &t->ci->hydro.force, t_force);
+ if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
/* Now, build all the dependencies for the hydro for the cells */
/* that are local and are not descendant of the same super_hydro-cells */
if (t->ci->nodeID == nodeID) {
- engine_make_hydro_loops_dependencies(sched, t, t2, t->ci, with_cooling);
- scheduler_addunlock(sched, t2, t->ci->super->end_force);
+ engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter,
+ t->ci, with_cooling, with_limiter);
+ scheduler_addunlock(sched, t_force, t->ci->super->end_force);
}
#endif
}
@@ -1519,56 +1631,73 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
#ifdef EXTRA_HYDRO_LOOP
/* Start by constructing the task for the second and third hydro loop */
- struct task *t2 =
+ t_gradient =
scheduler_addtask(sched, task_type_sub_pair, task_subtype_gradient,
t->flags, 0, t->ci, t->cj);
- struct task *t3 =
- scheduler_addtask(sched, task_type_sub_pair, task_subtype_force,
- t->flags, 0, t->ci, t->cj);
+ t_force = scheduler_addtask(sched, task_type_sub_pair, task_subtype_force,
+ t->flags, 0, t->ci, t->cj);
+
+ /* and the task for the time-step limiter */
+ if (with_limiter)
+ t_limiter =
+ scheduler_addtask(sched, task_type_sub_pair, task_subtype_limiter,
+ t->flags, 0, t->ci, t->cj);
/* Add the link between the new loop and both cells */
- engine_addlink(e, &t->ci->hydro.gradient, t2);
- engine_addlink(e, &t->cj->hydro.gradient, t2);
- engine_addlink(e, &t->ci->hydro.force, t3);
- engine_addlink(e, &t->cj->hydro.force, t3);
+ engine_addlink(e, &t->ci->hydro.gradient, t_gradient);
+ engine_addlink(e, &t->cj->hydro.gradient, t_gradient);
+ engine_addlink(e, &t->ci->hydro.force, t_force);
+ engine_addlink(e, &t->cj->hydro.force, t_force);
+ if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
+ if (with_limiter) engine_addlink(e, &t->cj->hydro.limiter, t_limiter);
/* Now, build all the dependencies for the hydro for the cells */
/* that are local and are not descendant of the same super_hydro-cells */
if (t->ci->nodeID == nodeID) {
- engine_make_hydro_loops_dependencies(sched, t, t2, t3, t->ci,
- with_cooling);
- scheduler_addunlock(sched, t3, t->ci->super->end_force);
+ engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
+ t_limiter, t->ci, with_cooling,
+ with_limiter);
+ scheduler_addunlock(sched, t_force, t->ci->super->end_force);
}
if (t->cj->nodeID == nodeID) {
if (t->ci->hydro.super != t->cj->hydro.super)
- engine_make_hydro_loops_dependencies(sched, t, t2, t3, t->cj,
- with_cooling);
+ engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
+ t_limiter, t->cj, with_cooling,
+ with_limiter);
if (t->ci->super != t->cj->super)
- scheduler_addunlock(sched, t3, t->cj->super->end_force);
+ scheduler_addunlock(sched, t_force, t->cj->super->end_force);
}
#else
/* Start by constructing the task for the second hydro loop */
- struct task *t2 =
- scheduler_addtask(sched, task_type_sub_pair, task_subtype_force,
- t->flags, 0, t->ci, t->cj);
+ t_force = scheduler_addtask(sched, task_type_sub_pair, task_subtype_force,
+ t->flags, 0, t->ci, t->cj);
+
+ /* and the task for the time-step limiter */
+ if (with_limiter)
+ t_limiter =
+ scheduler_addtask(sched, task_type_sub_pair, task_subtype_limiter,
+ t->flags, 0, t->ci, t->cj);
/* Add the link between the new loop and both cells */
- engine_addlink(e, &t->ci->hydro.force, t2);
- engine_addlink(e, &t->cj->hydro.force, t2);
+ engine_addlink(e, &t->ci->hydro.force, t_force);
+ engine_addlink(e, &t->cj->hydro.force, t_force);
+ if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
+ if (with_limiter) engine_addlink(e, &t->cj->hydro.limiter, t_limiter);
/* Now, build all the dependencies for the hydro for the cells */
/* that are local and are not descendant of the same super_hydro-cells */
if (t->ci->nodeID == nodeID) {
- engine_make_hydro_loops_dependencies(sched, t, t2, t->ci, with_cooling);
- scheduler_addunlock(sched, t2, t->ci->super->end_force);
+ engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter,
+ t->ci, with_cooling, with_limiter);
+ scheduler_addunlock(sched, t_force, t->ci->super->end_force);
}
if (t->cj->nodeID == nodeID) {
if (t->ci->hydro.super != t->cj->hydro.super)
- engine_make_hydro_loops_dependencies(sched, t, t2, t->cj,
- with_cooling);
+ engine_make_hydro_loops_dependencies(
+ sched, t, t_force, t_limiter, t->cj, with_cooling, with_limiter);
if (t->ci->super != t->cj->super)
- scheduler_addunlock(sched, t2, t->cj->super->end_force);
+ scheduler_addunlock(sched, t_force, t->cj->super->end_force);
}
#endif
}
@@ -1952,6 +2081,7 @@ struct cell_type_pair {
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);
struct cell_type_pair *cell_type_pairs = (struct cell_type_pair *)map_data;
for (int k = 0; k < num_elements; k++) {
@@ -1960,7 +2090,7 @@ void engine_addtasks_send_mapper(void *map_data, int num_elements,
const int type = cell_type_pairs[k].type;
/* Add the send task for the particle timesteps. */
- engine_addtasks_send_timestep(e, ci, cj, NULL);
+ engine_addtasks_send_timestep(e, ci, cj, NULL, NULL, with_limiter);
/* Add the send tasks for the cells in the proxy that have a hydro
* connection. */
@@ -1979,6 +2109,7 @@ void engine_addtasks_send_mapper(void *map_data, int num_elements,
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);
struct cell_type_pair *cell_type_pairs = (struct cell_type_pair *)map_data;
for (int k = 0; k < num_elements; k++) {
@@ -1986,7 +2117,7 @@ void engine_addtasks_recv_mapper(void *map_data, int num_elements,
const int type = cell_type_pairs[k].type;
/* Add the recv task for the particle timesteps. */
- engine_addtasks_recv_timestep(e, ci, NULL);
+ engine_addtasks_recv_timestep(e, ci, NULL, NULL, with_limiter);
/* Add the recv tasks for the cells in the proxy that have a hydro
* connection. */
@@ -2076,6 +2207,7 @@ void engine_maketasks(struct engine *e) {
const size_t self_grav_tasks_per_cell = 125;
const size_t ext_grav_tasks_per_cell = 1;
const size_t stars_tasks_per_cell = 27;
+ const size_t limiter_tasks_per_cell = 27;
if (e->policy & engine_policy_hydro)
e->size_links += s->tot_cells * hydro_tasks_per_cell;
@@ -2085,6 +2217,8 @@ void engine_maketasks(struct engine *e) {
e->size_links += s->tot_cells * self_grav_tasks_per_cell;
if (e->policy & engine_policy_stars)
e->size_links += s->tot_cells * stars_tasks_per_cell;
+ if (e->policy & engine_policy_limiter)
+ e->size_links += s->tot_cells * limiter_tasks_per_cell;
/* Allocate the new link list */
if ((e->links = (struct link *)malloc(sizeof(struct link) * e->size_links)) ==
diff --git a/src/engine_marktasks.c b/src/engine_marktasks.c
index 9c7a783c2547899816842cf9a05163e75d329aa8..d9b1a4cd345b725c442018836487b77c30169325 100644
--- a/src/engine_marktasks.c
+++ b/src/engine_marktasks.c
@@ -69,6 +69,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
struct scheduler *s = (struct scheduler *)(((size_t *)extra_data)[2]);
struct engine *e = (struct engine *)((size_t *)extra_data)[0];
const int nodeID = e->nodeID;
+ const int with_limiter = e->policy & engine_policy_limiter;
for (int ind = 0; ind < num_elements; ind++) {
@@ -90,6 +91,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
if (cell_is_active_hydro(ci, e)) {
scheduler_activate(s, t);
cell_activate_drift_part(ci, s);
+ if (with_limiter) cell_activate_limiter(ci, s);
}
}
@@ -99,6 +101,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
if (cell_is_active_hydro(ci, e)) {
scheduler_activate(s, t);
cell_activate_subcell_hydro_tasks(ci, NULL, s);
+ if (with_limiter) cell_activate_limiter(ci, s);
}
}
@@ -111,6 +114,16 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
if (cell_is_active_hydro(ci, e)) 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);
+ }
+
+ else if (t->type == task_type_sub_self &&
+ t->subtype == task_subtype_limiter) {
+ if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
+ }
+
#ifdef EXTRA_HYDRO_LOOP
else if (t_type == task_type_self && t_subtype == task_subtype_gradient) {
if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
@@ -207,6 +220,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
/* Only activate tasks that involve a local active cell. */
if ((t_subtype == task_subtype_density ||
t_subtype == task_subtype_gradient ||
+ t_subtype == task_subtype_limiter ||
t_subtype == task_subtype_force) &&
((ci_active_hydro && ci_nodeID == nodeID) ||
(cj_active_hydro && cj_nodeID == nodeID))) {
@@ -226,6 +240,12 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
+ /* And the limiter */
+ if (ci->nodeID == nodeID && with_limiter)
+ cell_activate_limiter(ci, s);
+ if (cj->nodeID == nodeID && with_limiter)
+ cell_activate_limiter(cj, s);
+
/* Check the sorts and activate them if needed. */
cell_activate_hydro_sorts(ci, t->flags, s);
cell_activate_hydro_sorts(cj, t->flags, s);
@@ -560,6 +580,11 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
scheduler_activate(s, t);
}
+ /* Time-step limiter? */
+ else if (t->type == task_type_timestep) {
+ if (cell_is_active_hydro(t->ci, e)) scheduler_activate(s, t);
+ }
+
/* Subgrid tasks */
else if (t_type == task_type_cooling) {
if (cell_is_active_hydro(t->ci, e) || cell_is_active_gravity(t->ci, e))
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index e16d41a8b4fc50d07f5cac7dd39e8d245bed9923..84d5ee5a19f6c6e6552dfff523dfaf5dc0f8dbc2 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -750,6 +750,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
struct part *restrict p, struct xpart *restrict xp) {
p->time_bin = 0;
+ p->wakeup = time_bin_not_awake;
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2];
diff --git a/src/hydro/Gadget2/hydro_debug.h b/src/hydro/Gadget2/hydro_debug.h
index d0642a03a4c4eecb2da80fdae473948e460c5e31..aeb43ee5d68930debfa867dc856465ac9d22902a 100644
--- a/src/hydro/Gadget2/hydro_debug.h
+++ b/src/hydro/Gadget2/hydro_debug.h
@@ -27,14 +27,14 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"h=%.3e, wcount=%.3f, wcount_dh=%.3e, m=%.3e, dh_drho=%.3e, rho=%.3e, "
"P=%.3e, P_over_rho2=%.3e, S=%.3e, dS/dt=%.3e, c=%.3e\n"
"divV=%.3e, rotV=[%.3e,%.3e,%.3e], balsara=%.3e \n "
- "v_sig=%e dh/dt=%.3e time_bin=%d\n",
+ "v_sig=%e dh/dt=%.3e time_bin=%d wakeup=%d\n",
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
p->h, p->density.wcount, p->density.wcount_dh, p->mass, p->density.rho_dh,
p->rho, hydro_get_comoving_pressure(p), p->force.P_over_rho2, p->entropy,
p->entropy_dt, p->force.soundspeed, p->density.div_v, p->density.rot_v[0],
p->density.rot_v[1], p->density.rot_v[2], p->force.balsara,
- p->force.v_sig, p->force.h_dt, p->time_bin);
+ p->force.v_sig, p->force.h_dt, p->time_bin, p->wakeup);
}
#endif /* SWIFT_GADGET2_HYDRO_DEBUG_H */
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index a3c5e21dbdf8df60b25b01c0326c33c3a10d1bce..14af63309e9c9dc27034769beea2d311e7a1115e 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -1051,4 +1051,34 @@ runner_iact_nonsym_2_vec_force(
#endif
+/**
+ * @brief Timestep limiter loop
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_limiter(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, float a, float H) {
+
+ /* Nothing to do here if both particles are active */
+}
+
+/**
+ * @brief Timestep limiter loop (non-symmetric version)
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_limiter(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, float a, float H) {
+
+ /* Wake up the neighbour? */
+ if (pi->force.v_sig > const_limiter_max_v_sig_ratio * pj->force.v_sig) {
+
+ pj->wakeup = time_bin_awake;
+
+ // MATTHIEU
+ // if (pj->wakeup == time_bin_not_awake)
+ // pj->wakeup = time_bin_awake;
+ // else if (pj->wakeup > 0)
+ // pj->wakeup = -pj->wakeup;
+ }
+}
+
#endif /* SWIFT_GADGET2_HYDRO_IACT_H */
diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h
index c3170414a375a5dfd54029cdc500c1cbcb764f6f..099f3b769d1ffa946c5e83f7190f325c5a42074c 100644
--- a/src/hydro/Gadget2/hydro_part.h
+++ b/src/hydro/Gadget2/hydro_part.h
@@ -150,6 +150,9 @@ struct part {
/* Time-step length */
timebin_t time_bin;
+ /* Need waking-up ? */
+ char wakeup;
+
#ifdef SWIFT_DEBUG_CHECKS
/* Time of the last drift */
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 85f88d418bd46354f7a1cd3dd89b0e77b556b7d9..314bf0dad7b4377aebed52fd9234dd13957ce1f3 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -262,6 +262,8 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
p->viscosity.alpha_min);
io_write_attribute_f(h_grpsph, "Viscosity decay length", p->viscosity.length);
io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
+ io_write_attribute_f(h_grpsph, "Max v_sig ratio (limiter)",
+ const_limiter_max_v_sig_ratio);
}
#endif
diff --git a/src/kick.h b/src/kick.h
index e85c9de40d2084304bde108e6f5fa9c776fd3e8f..7fb2afc1fa251dff30e81c5b0093b6313d8c4b7c 100644
--- a/src/kick.h
+++ b/src/kick.h
@@ -85,8 +85,8 @@ __attribute__((always_inline)) INLINE static void kick_part(
if (p->ti_kick != ti_start)
error(
"particle has not been kicked to the current time p->ti_kick=%lld, "
- "ti_start=%lld, ti_end=%lld id=%lld",
- p->ti_kick, ti_start, ti_end, p->id);
+ "ti_start=%lld, ti_end=%lld id=%lld time_bin=%d wakeup=%d",
+ p->ti_kick, ti_start, ti_end, p->id, p->time_bin, p->wakeup);
p->ti_kick = ti_end;
#endif
diff --git a/src/runner.c b/src/runner.c
index 7dea47aa0820f375f85d574a11e3dfa3700d36c6..01df46ed6eb8128cb9fc67a4d3a9a9c1e77ee8fa 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -64,6 +64,7 @@
#include "task.h"
#include "timers.h"
#include "timestep.h"
+#include "timestep_limiter.h"
#include "tracers.h"
#define TASK_LOOP_DENSITY 0
@@ -94,6 +95,13 @@
#undef FUNCTION
#undef FUNCTION_TASK_LOOP
+/* Import the limiter loop functions. */
+#define FUNCTION limiter
+#define FUNCTION_TASK_LOOP TASK_LOOP_LIMITER
+#include "runner_doiact.h"
+#undef FUNCTION
+#undef FUNCTION_TASK_LOOP
+
/* Import the gravity loop functions. */
#include "runner_doiact_grav.h"
@@ -1673,19 +1681,26 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
/* If particle needs to be kicked */
if (part_is_starting(p, e)) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (p->wakeup == time_bin_awake)
+ error("Woken-up particle that has not been processed in kick1");
+#endif
+
+ /* Skip particles that have been woken up and treated by the limiter. */
+ if (p->wakeup != time_bin_not_awake) continue;
+
const integertime_t ti_step = get_integer_timestep(p->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(ti_current + 1, p->time_bin);
#ifdef SWIFT_DEBUG_CHECKS
- const integertime_t ti_end =
- get_integer_time_end(ti_current + 1, p->time_bin);
+ const integertime_t ti_end = ti_begin + ti_step;
if (ti_begin != ti_current)
error(
"Particle in wrong time-bin, ti_end=%lld, ti_begin=%lld, "
- "ti_step=%lld time_bin=%d ti_current=%lld",
- ti_end, ti_begin, ti_step, p->time_bin, ti_current);
+ "ti_step=%lld time_bin=%d wakeup=%d ti_current=%lld",
+ ti_end, ti_begin, ti_step, p->time_bin, p->wakeup, ti_current);
#endif
/* Time interval for this half-kick */
@@ -1847,39 +1862,60 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
/* If particle needs to be kicked */
if (part_is_active(p, e)) {
- const integertime_t ti_step = get_integer_timestep(p->time_bin);
- const integertime_t ti_begin =
- get_integer_time_begin(ti_current, p->time_bin);
+ integertime_t ti_begin, ti_end, ti_step;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (p->wakeup == time_bin_awake)
+ error("Woken-up particle that has not been processed in kick1");
+#endif
+
+ if (p->wakeup == time_bin_not_awake) {
+
+ /* Time-step from a regular kick */
+ ti_step = get_integer_timestep(p->time_bin);
+ ti_begin = get_integer_time_begin(ti_current, p->time_bin);
+ ti_end = ti_begin + ti_step;
+
+ } else {
+
+ /* Time-step that follows a wake-up call */
+ ti_begin = get_integer_time_begin(ti_current, p->wakeup);
+ ti_end = get_integer_time_end(ti_current, p->time_bin);
+ ti_step = ti_end - ti_begin;
+
+ /* Reset the flag. Everything is back to normal from now on. */
+ p->wakeup = time_bin_awake;
+ }
#ifdef SWIFT_DEBUG_CHECKS
if (ti_begin + ti_step != ti_current)
error(
"Particle in wrong time-bin, ti_begin=%lld, ti_step=%lld "
- "time_bin=%d ti_current=%lld",
- ti_begin, ti_step, p->time_bin, ti_current);
+ "time_bin=%d wakeup=%d ti_current=%lld",
+ ti_begin, ti_step, p->time_bin, p->wakeup, ti_current);
#endif
/* Time interval for this half-kick */
double dt_kick_grav, dt_kick_hydro, dt_kick_therm, dt_kick_corr;
if (with_cosmology) {
dt_kick_hydro = cosmology_get_hydro_kick_factor(
- cosmo, ti_begin + ti_step / 2, ti_begin + ti_step);
+ cosmo, ti_begin + ti_step / 2, ti_end);
dt_kick_grav = cosmology_get_grav_kick_factor(
- cosmo, ti_begin + ti_step / 2, ti_begin + ti_step);
+ cosmo, ti_begin + ti_step / 2, ti_end);
dt_kick_therm = cosmology_get_therm_kick_factor(
- cosmo, ti_begin + ti_step / 2, ti_begin + ti_step);
+ cosmo, ti_begin + ti_step / 2, ti_end);
dt_kick_corr = cosmology_get_corr_kick_factor(
- cosmo, ti_begin + ti_step / 2, ti_begin + ti_step);
+ cosmo, ti_begin + ti_step / 2, ti_end);
} else {
- dt_kick_hydro = (ti_step / 2) * time_base;
- dt_kick_grav = (ti_step / 2) * time_base;
- dt_kick_therm = (ti_step / 2) * time_base;
- dt_kick_corr = (ti_step / 2) * time_base;
+ dt_kick_hydro = (ti_end - (ti_begin + ti_step / 2)) * time_base;
+ dt_kick_grav = (ti_end - (ti_begin + ti_step / 2)) * time_base;
+ dt_kick_therm = (ti_end - (ti_begin + ti_step / 2)) * time_base;
+ dt_kick_corr = (ti_end - (ti_begin + ti_step / 2)) * time_base;
}
/* Finish the time-step with a second half-kick */
kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm,
dt_kick_corr, cosmo, hydro_props, ti_begin + ti_step / 2,
- ti_begin + ti_step);
+ ti_end);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that kick and the drift are synchronized */
@@ -2281,6 +2317,132 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
if (timer) TIMER_TOC(timer_timestep);
}
+/**
+ * @brief Apply the time-step limiter to all awaken particles in a cell
+ * hierarchy.
+ *
+ * @param r The task #runner.
+ * @param c The #cell.
+ * @param force Limit the particles irrespective of the #cell flags.
+ * @param timer Are we timing this ?
+ */
+void runner_do_limiter(struct runner *r, struct cell *c, int force, int timer) {
+
+ const struct engine *e = r->e;
+ const integertime_t ti_current = e->ti_current;
+ const int count = c->hydro.count;
+ struct part *restrict parts = c->hydro.parts;
+ struct xpart *restrict xparts = c->hydro.xparts;
+
+ TIMER_TIC;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that we only limit local cells. */
+ if (c->nodeID != engine_rank) error("Limiting dt of a foreign cell is nope.");
+#endif
+
+ integertime_t ti_hydro_end_min = max_nr_timesteps, ti_hydro_end_max = 0,
+ ti_hydro_beg_max = 0;
+ integertime_t ti_gravity_end_min = max_nr_timesteps, ti_gravity_end_max = 0,
+ ti_gravity_beg_max = 0;
+
+ /* Limit irrespective of cell flags? */
+ force |= c->hydro.do_limiter;
+
+ /* Loop over the progeny ? */
+ if (c->split && (force || c->hydro.do_sub_limiter)) {
+ for (int k = 0; k < 8; k++) {
+ if (c->progeny[k] != NULL) {
+ struct cell *restrict cp = c->progeny[k];
+
+ /* Recurse */
+ runner_do_limiter(r, cp, force, 0);
+
+ /* And aggregate */
+ ti_hydro_end_min = min(cp->hydro.ti_end_min, ti_hydro_end_min);
+ ti_hydro_end_max = max(cp->hydro.ti_end_max, ti_hydro_end_max);
+ ti_hydro_beg_max = max(cp->hydro.ti_beg_max, ti_hydro_beg_max);
+ ti_gravity_end_min = min(cp->grav.ti_end_min, ti_gravity_end_min);
+ ti_gravity_end_max = max(cp->grav.ti_end_max, ti_gravity_end_max);
+ ti_gravity_beg_max = max(cp->grav.ti_beg_max, ti_gravity_beg_max);
+ }
+ }
+
+ /* Store the updated values */
+ c->hydro.ti_end_min = min(c->hydro.ti_end_min, ti_hydro_end_min);
+ c->hydro.ti_end_max = max(c->hydro.ti_end_max, ti_hydro_end_max);
+ c->hydro.ti_beg_max = max(c->hydro.ti_beg_max, ti_hydro_beg_max);
+ c->grav.ti_end_min = min(c->grav.ti_end_min, ti_gravity_end_min);
+ c->grav.ti_end_max = max(c->grav.ti_end_max, ti_gravity_end_max);
+ c->grav.ti_beg_max = max(c->grav.ti_beg_max, ti_gravity_beg_max);
+
+ } else if (!c->split && force) {
+
+ ti_hydro_end_min = c->hydro.ti_end_min;
+ ti_hydro_end_max = c->hydro.ti_end_max;
+ ti_hydro_beg_max = c->hydro.ti_beg_max;
+ ti_gravity_end_min = c->grav.ti_end_min;
+ ti_gravity_end_max = c->grav.ti_end_max;
+ ti_gravity_beg_max = c->grav.ti_beg_max;
+
+ /* Loop over the gas particles in this cell. */
+ for (int k = 0; k < count; k++) {
+
+ /* Get a handle on the part. */
+ struct part *restrict p = &parts[k];
+ struct xpart *restrict xp = &xparts[k];
+
+ /* If the particle will be active no need to wake it up */
+ if (part_is_active(p, e) && p->wakeup != time_bin_not_awake)
+ p->wakeup = time_bin_not_awake;
+
+ /* Bip, bip, bip... wake-up time */
+ if (p->wakeup == time_bin_awake) {
+
+ /* Apply the limiter and get the new time-step size */
+ const integertime_t ti_new_step = timestep_limit_part(p, xp, e);
+
+ /* What is the next sync-point ? */
+ ti_hydro_end_min = min(ti_current + ti_new_step, ti_hydro_end_min);
+ ti_hydro_end_max = max(ti_current + ti_new_step, ti_hydro_end_max);
+
+ /* What is the next starting point for this cell ? */
+ ti_hydro_beg_max = max(ti_current, ti_hydro_beg_max);
+
+ /* Also limit the gpart counter-part */
+ if (p->gpart != NULL) {
+
+ /* Register the time-bin */
+ p->gpart->time_bin = p->time_bin;
+
+ /* What is the next sync-point ? */
+ ti_gravity_end_min =
+ min(ti_current + ti_new_step, ti_gravity_end_min);
+ ti_gravity_end_max =
+ max(ti_current + ti_new_step, ti_gravity_end_max);
+
+ /* What is the next starting point for this cell ? */
+ ti_gravity_beg_max = max(ti_current, ti_gravity_beg_max);
+ }
+ }
+ }
+
+ /* Store the updated values */
+ c->hydro.ti_end_min = min(c->hydro.ti_end_min, ti_hydro_end_min);
+ c->hydro.ti_end_max = max(c->hydro.ti_end_max, ti_hydro_end_max);
+ c->hydro.ti_beg_max = max(c->hydro.ti_beg_max, ti_hydro_beg_max);
+ c->grav.ti_end_min = min(c->grav.ti_end_min, ti_gravity_end_min);
+ c->grav.ti_end_max = max(c->grav.ti_end_max, ti_gravity_end_max);
+ c->grav.ti_beg_max = max(c->grav.ti_beg_max, ti_gravity_beg_max);
+ }
+
+ /* Clear the limiter flags. */
+ c->hydro.do_limiter = 0;
+ c->hydro.do_sub_limiter = 0;
+
+ if (timer) TIMER_TOC(timer_do_limiter);
+}
+
/**
* @brief End the force calculation of all active particles in a cell
* by multiplying the acccelerations by the relevant constants
@@ -2733,6 +2895,8 @@ void *runner_main(void *data) {
#endif
else if (t->subtype == task_subtype_force)
runner_doself2_branch_force(r, ci);
+ else if (t->subtype == task_subtype_limiter)
+ runner_doself2_branch_limiter(r, ci);
else if (t->subtype == task_subtype_grav)
runner_doself_recursive_grav(r, ci, 1);
else if (t->subtype == task_subtype_external_grav)
@@ -2754,6 +2918,8 @@ void *runner_main(void *data) {
#endif
else if (t->subtype == task_subtype_force)
runner_dopair2_branch_force(r, ci, cj);
+ else if (t->subtype == task_subtype_limiter)
+ runner_dopair2_branch_limiter(r, ci, cj);
else if (t->subtype == task_subtype_grav)
runner_dopair_recursive_grav(r, ci, cj, 1);
else if (t->subtype == task_subtype_stars_density)
@@ -2773,6 +2939,8 @@ void *runner_main(void *data) {
#endif
else if (t->subtype == task_subtype_force)
runner_dosub_self2_force(r, ci, 1);
+ else if (t->subtype == task_subtype_limiter)
+ runner_dosub_self2_limiter(r, ci, 1);
else if (t->subtype == task_subtype_stars_density)
runner_dosub_self_stars_density(r, ci, 1);
else if (t->subtype == task_subtype_stars_feedback)
@@ -2790,6 +2958,8 @@ void *runner_main(void *data) {
#endif
else if (t->subtype == task_subtype_force)
runner_dosub_pair2_force(r, ci, cj, t->flags, 1);
+ else if (t->subtype == task_subtype_limiter)
+ runner_dosub_pair2_limiter(r, ci, cj, t->flags, 1);
else if (t->subtype == task_subtype_stars_density)
runner_dosub_pair_stars_density(r, ci, cj, t->flags, 1);
else if (t->subtype == task_subtype_stars_feedback)
@@ -2849,6 +3019,9 @@ void *runner_main(void *data) {
case task_type_timestep:
runner_do_timestep(r, ci, 1);
break;
+ case task_type_timestep_limiter:
+ runner_do_limiter(r, ci, 0, 1);
+ break;
#ifdef WITH_MPI
case task_type_send:
if (t->subtype == task_subtype_tend) {
@@ -2865,6 +3038,8 @@ void *runner_main(void *data) {
runner_do_recv_part(r, ci, 0, 1);
} else if (t->subtype == task_subtype_gradient) {
runner_do_recv_part(r, ci, 0, 1);
+ } else if (t->subtype == task_subtype_limiter) {
+ runner_do_recv_part(r, ci, 0, 1);
} else if (t->subtype == task_subtype_gpart) {
runner_do_recv_gpart(r, ci, 1);
} else if (t->subtype == task_subtype_spart) {
diff --git a/src/scheduler.c b/src/scheduler.c
index 2ae8f6785434af021b52dd2d6586b4e2dc5d68bb..70705e3ee48fc3a038bd3ecc835b77017588f24f 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -263,6 +263,7 @@ void scheduler_write_dependencies(struct scheduler *s, int verbose) {
int density_cluster[4] = {0};
int gradient_cluster[4] = {0};
int force_cluster[4] = {0};
+ int limiter_cluster[4] = {0};
int gravity_cluster[5] = {0};
int stars_density_cluster[4] = {0};
@@ -283,6 +284,8 @@ void scheduler_write_dependencies(struct scheduler *s, int verbose) {
gradient_cluster[k] = 1;
if (type == task_type_self + k && subtype == task_subtype_force)
force_cluster[k] = 1;
+ if (type == task_type_self + k && subtype == task_subtype_limiter)
+ limiter_cluster[k] = 1;
if (type == task_type_self + k && subtype == task_subtype_grav)
gravity_cluster[k] = 1;
if (type == task_type_self + k &&
@@ -323,8 +326,17 @@ void scheduler_write_dependencies(struct scheduler *s, int verbose) {
subtaskID_names[task_subtype_gradient]);
fprintf(f, "\t};\n");
+ /* Make a cluster for the limiter tasks */
+ fprintf(f, "\t subgraph cluster2{\n");
+ fprintf(f, "\t\t label=\"\";\n");
+ for (int k = 0; k < 4; ++k)
+ if (limiter_cluster[k])
+ fprintf(f, "\t\t \"%s %s\";\n", taskID_names[task_type_self + k],
+ subtaskID_names[task_subtype_limiter]);
+ fprintf(f, "\t};\n");
+
/* Make a cluster for the gravity tasks */
- fprintf(f, "\t subgraph cluster3{\n");
+ fprintf(f, "\t subgraph cluster4{\n");
fprintf(f, "\t\t label=\"\";\n");
for (int k = 0; k < 2; ++k)
if (gravity_cluster[k])
@@ -1948,6 +1960,9 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
case task_type_timestep:
cost = wscale * count_i + wscale * gcount_i;
break;
+ case task_type_timestep_limiter:
+ cost = wscale * count_i;
+ break;
case task_type_send:
#if defined(WITH_MPI) && (defined(HAVE_PARMETIS) || defined(HAVE_METIS))
partcost = 0;
@@ -2157,7 +2172,8 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
subtaskMPI_comms[t->subtype], &t->req);
} else if (t->subtype == task_subtype_xv ||
t->subtype == task_subtype_rho ||
- t->subtype == task_subtype_gradient) {
+ t->subtype == task_subtype_gradient ||
+ t->subtype == task_subtype_limiter) {
err = MPI_Irecv(t->ci->hydro.parts, t->ci->hydro.count, part_mpi_type,
t->ci->nodeID, t->flags, subtaskMPI_comms[t->subtype],
&t->req);
@@ -2208,7 +2224,8 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
subtaskMPI_comms[t->subtype], &t->req);
} else if (t->subtype == task_subtype_xv ||
t->subtype == task_subtype_rho ||
- t->subtype == task_subtype_gradient) {
+ t->subtype == task_subtype_gradient ||
+ t->subtype == task_subtype_limiter) {
if ((t->ci->hydro.count * sizeof(struct part)) > s->mpi_message_limit)
err = MPI_Isend(t->ci->hydro.parts, t->ci->hydro.count,
part_mpi_type, t->cj->nodeID, t->flags,
diff --git a/src/space.c b/src/space.c
index 6401be586b784bd92684c2f70bc60f073b16148b..40cbd10688c65a4c93685894732d18c6ce178e2a 100644
--- a/src/space.c
+++ b/src/space.c
@@ -189,6 +189,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->hydro.density = NULL;
c->hydro.gradient = NULL;
c->hydro.force = NULL;
+ c->hydro.limiter = NULL;
c->grav.grav = NULL;
c->grav.mm = NULL;
c->hydro.dx_max_part = 0.0f;
@@ -223,12 +224,12 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->kick1 = NULL;
c->kick2 = NULL;
c->timestep = NULL;
+ c->timestep_limiter = NULL;
c->end_force = NULL;
c->hydro.drift = NULL;
c->grav.drift = NULL;
c->grav.drift_out = NULL;
c->hydro.cooling = NULL;
- c->sourceterms = NULL;
c->grav.long_range = NULL;
c->grav.down_in = NULL;
c->grav.down = NULL;
@@ -272,12 +273,14 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->mpi.hydro.recv_gradient = NULL;
c->mpi.grav.recv = NULL;
c->mpi.recv_ti = NULL;
+ c->mpi.limiter.recv = NULL;
c->mpi.hydro.send_xv = NULL;
c->mpi.hydro.send_rho = NULL;
c->mpi.hydro.send_gradient = NULL;
c->mpi.grav.send = NULL;
c->mpi.send_ti = NULL;
+ c->mpi.limiter.send = NULL;
#endif
}
}
@@ -4268,6 +4271,44 @@ void space_check_timesteps(struct space *s) {
#endif
}
+/**
+ * @brief #threadpool mapper function for the limiter debugging check
+ */
+void space_check_limiter_mapper(void *map_data, int nr_parts,
+ void *extra_data) {
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Unpack the data */
+ struct part *restrict parts = (struct part *)map_data;
+
+ /* Verify that all limited particles have been treated */
+ for (int k = 0; k < nr_parts; k++) {
+ if (parts[k].wakeup == time_bin_awake) error("Particle still woken up!");
+ if (parts[k].gpart != NULL)
+ if (parts[k].time_bin != parts[k].gpart->time_bin)
+ error("Gpart not on the same time-bin as part");
+ }
+#else
+ error("Calling debugging code without debugging flag activated.");
+#endif
+}
+
+/**
+ * @brief Checks that all particles have their wakeup flag in a correct state.
+ *
+ * Should only be used for debugging purposes.
+ *
+ * @param s The #space to check.
+ */
+void space_check_limiter(struct space *s) {
+#ifdef SWIFT_DEBUG_CHECKS
+
+ threadpool_map(&s->e->threadpool, space_check_limiter_mapper, s->parts,
+ s->nr_parts, sizeof(struct part), 1000, NULL);
+#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 a1280945d2aa232cbb5e5b519266bc7058e5dc57..84e71c9b55810d1075191df6132db09e2a0c8796 100644
--- a/src/space.h
+++ b/src/space.h
@@ -317,6 +317,7 @@ void space_check_drift_point(struct space *s, integertime_t ti_drift,
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_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 83db23645006f4d451c14ba0b9c2f997c9753655..f27b5fd413c1c216afb7d777244b0cee88159e6e 100644
--- a/src/task.c
+++ b/src/task.c
@@ -66,6 +66,7 @@ const char *taskID_names[task_type_count] = {"none",
"kick1",
"kick2",
"timestep",
+ "timestep_limiter",
"send",
"recv",
"grav_long_range",
@@ -83,10 +84,10 @@ const char *taskID_names[task_type_count] = {"none",
/* Sub-task type names. */
const char *subtaskID_names[task_subtype_count] = {
- "none", "density", "gradient", "force",
- "grav", "external_grav", "tend", "xv",
- "rho", "gpart", "multipole", "spart",
- "stars_density", "stars_feedback"};
+ "none", "density", "gradient", "force",
+ "limiter", "grav", "external_grav", "tend",
+ "xv", "rho", "gpart", "multipole",
+ "spart", "stars_density", "stars_feedback"};
#ifdef WITH_MPI
/* MPI communicators for the subtypes. */
@@ -140,6 +141,7 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
case task_type_sort:
case task_type_ghost:
case task_type_extra_ghost:
+ case task_type_timestep_limiter:
case task_type_cooling:
return task_action_part;
break;
@@ -161,6 +163,7 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
case task_subtype_density:
case task_subtype_gradient:
case task_subtype_force:
+ case task_subtype_limiter:
return task_action_part;
break;
@@ -337,6 +340,7 @@ void task_unlock(struct task *t) {
case task_type_drift_part:
case task_type_sort:
+ case task_type_timestep_limiter:
cell_unlocktree(ci);
break;
@@ -462,6 +466,7 @@ int task_lock(struct task *t) {
case task_type_drift_part:
case task_type_sort:
+ case task_type_timestep_limiter:
if (ci->hydro.hold) return 0;
if (cell_locktree(ci) != 0) return 0;
break;
diff --git a/src/task.h b/src/task.h
index e432a051f5f4aff6bddd5b3a91df50160ae4912f..83c71a73e23da8d25c99bf038045cde12aa64c89 100644
--- a/src/task.h
+++ b/src/task.h
@@ -58,6 +58,7 @@ enum task_types {
task_type_kick1,
task_type_kick2,
task_type_timestep,
+ task_type_timestep_limiter,
task_type_send,
task_type_recv,
task_type_grav_long_range,
@@ -83,6 +84,7 @@ enum task_subtypes {
task_subtype_density,
task_subtype_gradient,
task_subtype_force,
+ task_subtype_limiter,
task_subtype_grav,
task_subtype_external_grav,
task_subtype_tend,
diff --git a/src/timestep_limiter.h b/src/timestep_limiter.h
new file mode 100644
index 0000000000000000000000000000000000000000..0e1f68f00907c326654146d4d8f2fc4367f0a29f
--- /dev/null
+++ b/src/timestep_limiter.h
@@ -0,0 +1,128 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_TIMESTEP_LIMITER_H
+#define SWIFT_TIMESTEP_LIMITER_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/**
+ * @brief Wakes up a particle by rewinding it's kick1 back in time and applying
+ * a new one such that the particle becomes active again in the next time-step.
+ *
+ * @param p The #part to update.
+ * @param xp Its #xpart companion.
+ * @param e The #engine (to extract time-line information).
+ */
+__attribute__((always_inline)) INLINE static integertime_t timestep_limit_part(
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct engine *e) {
+
+ const int with_cosmology = e->policy & engine_policy_cosmology;
+ const double time_base = e->time_base;
+
+ integertime_t old_ti_beg, old_ti_end;
+ timebin_t old_time_bin;
+
+ /* Let's see when this particle started and used to end */
+ if (p->wakeup == time_bin_awake) {
+
+ /* Normal case */
+ old_ti_beg = get_integer_time_begin(e->ti_current, p->time_bin);
+ old_ti_end = get_integer_time_end(e->ti_current, p->time_bin);
+ old_time_bin = p->time_bin;
+ } else {
+
+ /* Particle that was limited in the previous step already */
+ old_ti_beg = get_integer_time_begin(e->ti_current, -p->wakeup);
+ old_ti_end = get_integer_time_end(e->ti_current, p->time_bin);
+ old_time_bin = -p->wakeup;
+ }
+
+ const integertime_t old_dti = old_ti_end - old_ti_beg;
+
+ /* The new fake time-step the particle will be on */
+ const integertime_t new_fake_ti_step =
+ get_integer_timestep(e->min_active_bin);
+
+ /* The actual time-step size this particle will use */
+ const integertime_t new_ti_beg = old_ti_beg;
+ const integertime_t new_ti_end = e->ti_current + new_fake_ti_step;
+ const integertime_t new_dti = new_ti_end - new_ti_beg;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Some basic safety checks */
+ if (old_ti_beg >= e->ti_current)
+ error(
+ "Incorrect value for old time-step beginning ti_current=%lld, "
+ "old_ti_beg=%lld",
+ e->ti_current, old_ti_beg);
+
+ if (old_ti_end <= e->ti_current)
+ error(
+ "Incorrect value for old time-step end ti_current=%lld, "
+ "old_ti_end=%lld",
+ e->ti_current, old_ti_end);
+
+ if (new_ti_end > old_ti_end) error("New end of time-step after the old one");
+
+ if (new_dti > old_dti) error("New time-step larger than old one");
+
+ if (new_fake_ti_step == 0) error("Wakeup call too early");
+#endif
+
+ double dt_kick_grav = 0., dt_kick_hydro = 0., dt_kick_therm = 0.,
+ dt_kick_corr = 0.;
+
+ /* Now we need to reverse the kick1... (the dt are negative here) */
+ if (with_cosmology) {
+ error("aa");
+ } else {
+ dt_kick_hydro = -(old_dti / 2) * time_base;
+ dt_kick_grav = -(old_dti / 2) * time_base;
+ dt_kick_therm = -(old_dti / 2) * time_base;
+ dt_kick_corr = -(old_dti / 2) * time_base;
+ }
+ kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm, dt_kick_corr,
+ e->cosmology, e->hydro_properties, old_ti_beg + old_dti / 2,
+ old_ti_beg);
+
+ /* ...and apply the new one (dt is positiive) */
+ if (with_cosmology) {
+ error("aa");
+ } else {
+ dt_kick_hydro = (new_dti / 2) * time_base;
+ dt_kick_grav = (new_dti / 2) * time_base;
+ dt_kick_therm = (new_dti / 2) * time_base;
+ dt_kick_corr = (new_dti / 2) * time_base;
+ }
+ kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm, dt_kick_corr,
+ e->cosmology, e->hydro_properties, new_ti_beg,
+ new_ti_beg + new_dti / 2);
+
+ /* Remember the old time-bin */
+ p->wakeup = old_time_bin;
+
+ /* Update the time bin of this particle */
+ p->time_bin = e->min_active_bin;
+
+ return new_fake_ti_step;
+}
+
+#endif /* SWIFT_TIMESTEP_LIMITER_H */
diff --git a/tools/task_plots/analyse_tasks.py b/tools/task_plots/analyse_tasks.py
index 5738ca068c215a78c6fb4ef2524ce3d73565633e..e897424a95be8937073bd16adf108fa4fa1456ad 100755
--- a/tools/task_plots/analyse_tasks.py
+++ b/tools/task_plots/analyse_tasks.py
@@ -82,6 +82,7 @@ TASKTYPES = [
"kick1",
"kick2",
"timestep",
+ "timestep_limiter",
"send",
"recv",
"grav_long_range",
@@ -104,6 +105,7 @@ SUBTYPES = [
"density",
"gradient",
"force",
+ "limiter",
"grav",
"external_grav",
"tend",
diff --git a/tools/task_plots/plot_tasks.py b/tools/task_plots/plot_tasks.py
index 82dc882becfc2a7a8a537b822aceb8d9d226792d..12fd4d241a268c9d45fd72f5cdda2727221ba94d 100755
--- a/tools/task_plots/plot_tasks.py
+++ b/tools/task_plots/plot_tasks.py
@@ -167,6 +167,7 @@ TASKTYPES = [
"kick1",
"kick2",
"timestep",
+ "timestep_limiter",
"send",
"recv",
"grav_long_range",
@@ -189,6 +190,7 @@ SUBTYPES = [
"density",
"gradient",
"force",
+ "limiter",
"grav",
"external_grav",
"tend",
@@ -204,15 +206,23 @@ SUBTYPES = [
# Task/subtypes of interest.
FULLTYPES = [
+ "self/limiter",
"self/force",
+ "self/gradient",
"self/density",
"self/grav",
+ "sub_self/limiter",
"sub_self/force",
+ "sub_self/gradient",
"sub_self/density",
+ "pair/limiter",
"pair/force",
+ "pair/gradient",
"pair/density",
"pair/grav",
+ "sub_pair/limiter",
"sub_pair/force",
+ "sub_pair/gradient",
"sub_pair/density",
"recv/xv",
"send/xv",