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..b4252581d6eb3b2932a074e7545b2d308be51865 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-3 # Time between statistics output
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/examples/SedovBlast_2D/run.sh b/examples/SedovBlast_2D/run.sh
index b481d4555241c17015452a2139c04c541ccf1cdc..e2136f8f5e6ee9bde61d5189ed7955d53a3a9a6e 100755
--- a/examples/SedovBlast_2D/run.sh
+++ b/examples/SedovBlast_2D/run.sh
@@ -13,7 +13,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_2D/sedov.yml b/examples/SedovBlast_2D/sedov.yml
index 84177ece31ef98ec55c41513276c9c0158e69bcf..b4252581d6eb3b2932a074e7545b2d308be51865 100644
--- a/examples/SedovBlast_2D/sedov.yml
+++ b/examples/SedovBlast_2D/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-4 # 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:
diff --git a/examples/SedovBlast_3D/run.sh b/examples/SedovBlast_3D/run.sh
index 88aec36a7b96b5fd2a7fde41f0e0c9dc7185f1e8..7f0788cc822f1a6427fb6dbee4a921f79c942808 100755
--- a/examples/SedovBlast_3D/run.sh
+++ b/examples/SedovBlast_3D/run.sh
@@ -13,7 +13,7 @@ then
fi
# Run SWIFT
-../swift --hydro --threads=4 sedov.yml 2>&1 | tee output.log
+../swift --hydro --limiter --threads=4 sedov.yml 2>&1 | tee output.log
# Plot the solution
python plotSolution.py 5
diff --git a/examples/SedovBlast_3D/sedov.yml b/examples/SedovBlast_3D/sedov.yml
index 6cf5b02427b8004787b646e6bcdd4bacaa25bc06..19e8c72538a748304ca4da076458c9ae27dc8f46 100644
--- a/examples/SedovBlast_3D/sedov.yml
+++ b/examples/SedovBlast_3D/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-4 # 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:
diff --git a/examples/main.c b/examples/main.c
index b5b41a2c48b1f9abffcdd2af32396c037bb3ce80..8327248f1bafac23981d59871cf5ce83987528a7 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)
@@ -896,6 +900,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 a7914a8a5a20d596a5516d61959e5c826c737b15..0f4103673e9a6ba17eafa3c2ed859dd80e4b5425 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.
*/
@@ -1622,7 +1633,10 @@ void cell_activate_drift_part(struct cell *c, struct scheduler *s) {
for (struct cell *parent = c->parent;
parent != NULL && !parent->hydro.do_sub_drift;
parent = parent->parent) {
+
+ /* Mark this cell for drifting */
parent->hydro.do_sub_drift = 1;
+
if (parent == c->hydro.super) {
#ifdef SWIFT_DEBUG_CHECKS
if (parent->hydro.drift == NULL)
@@ -1686,6 +1700,45 @@ 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 limiting. */
+ 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) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->timestep_limiter == NULL)
+ error("Trying to activate un-existing c->timestep_limiter");
+#endif
+ scheduler_activate(s, c->timestep_limiter);
+ } else {
+ for (struct cell *parent = c->parent;
+ parent != NULL && !parent->hydro.do_sub_limiter;
+ parent = parent->parent) {
+
+ /* Mark this cell for limiting */
+ parent->hydro.do_sub_limiter = 1;
+
+ if (parent == c->super) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (parent->timestep_limiter == NULL)
+ error("Trying to activate un-existing parent->timestep_limiter");
+#endif
+ scheduler_activate(s, parent->timestep_limiter);
+ break;
+ }
+ }
+ }
+}
+
/**
* @brief Activate the sorts up a cell hierarchy.
*/
@@ -1816,6 +1869,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 +1903,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 +2209,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 +2779,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 +2805,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 +2820,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 +2858,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 +2872,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 +2885,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 +2903,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 +2917,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 +2931,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 +2945,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 +2960,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 4c9fce2f43323f91c4f58fb3440b794a972cfec2..1cc5527751707ce497788d7798fe9144a6c8ee1a 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;
@@ -1908,6 +1909,10 @@ int engine_estimate_nr_tasks(const 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;
@@ -2586,8 +2591,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;
@@ -2595,6 +2603,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);
}
/**
@@ -2804,6 +2813,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);
@@ -3061,6 +3075,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 db4a67d8209295f9bbd230cd989a49ad35fc5ca4..757f2b00add656ed2ae7b40ebe431f0487cb9ad2 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 9aea036b7f2cc0c31f6e41ea18789450f3e5e689..adcb3e69707691c5218ebc9834f1e8dd0226c80a 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,16 @@ 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);
+ }
+
#if defined(WITH_LOGGER)
scheduler_addunlock(s, c->kick1, c->logger);
#endif
@@ -1281,7 +1332,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 */
@@ -1299,14 +1351,15 @@ 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);
@@ -1347,6 +1400,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];
@@ -1364,31 +1423,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 = 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);
+ if (with_limiter)
+ scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
#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);
+ if (with_limiter)
+ scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
#endif
}
@@ -1407,54 +1488,103 @@ 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 (with_limiter)
+ scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter, t->ci->super->timestep_limiter);
}
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);
- if (t->ci->super != t->cj->super)
- scheduler_addunlock(sched, t3, t->cj->super->end_force);
+ if (t->ci->hydro.super != t->cj->hydro.super) {
+ 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, t_force, t->cj->super->end_force);
+ if (with_limiter)
+ scheduler_addunlock(sched, t->cj->super->kick2, t_limiter);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter, t->cj->super->timestep);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter,
+ t->cj->super->timestep_limiter);
+ }
}
#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 (with_limiter)
+ scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter, t->ci->super->timestep_limiter);
}
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);
- if (t->ci->super != t->cj->super)
- scheduler_addunlock(sched, t2, t->cj->super->end_force);
+ if (t->ci->hydro.super != t->cj->hydro.super) {
+ 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, t_force, t->cj->super->end_force);
+ if (with_limiter)
+ scheduler_addunlock(sched, t->cj->super->kick2, t_limiter);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter, t->cj->super->timestep);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter,
+ t->cj->super->timestep_limiter);
+ }
}
#endif
@@ -1472,39 +1602,65 @@ 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);
+ if (with_limiter)
+ scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter, t->ci->super->timestep_limiter);
}
#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);
+ if (with_limiter)
+ scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter, t->ci->super->timestep_limiter);
}
#endif
}
@@ -1526,56 +1682,106 @@ 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 (with_limiter)
+ scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter, t->ci->super->timestep_limiter);
}
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);
- if (t->ci->super != t->cj->super)
- scheduler_addunlock(sched, t3, t->cj->super->end_force);
+ if (t->ci->hydro.super != t->cj->hydro.super) {
+ 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, t_force, t->cj->super->end_force);
+ if (with_limiter)
+ scheduler_addunlock(sched, t->cj->super->kick2, t_limiter);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter, t->cj->super->timestep);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter,
+ t->cj->super->timestep_limiter);
+ }
}
#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 (with_limiter)
+ scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter, t->ci->super->timestep_limiter);
}
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);
- if (t->ci->super != t->cj->super)
- scheduler_addunlock(sched, t2, t->cj->super->end_force);
+ if (t->ci->hydro.super != t->cj->hydro.super) {
+ 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, t_force, t->cj->super->end_force);
+ if (with_limiter)
+ scheduler_addunlock(sched, t->cj->super->kick2, t_limiter);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter, t->cj->super->timestep);
+ if (with_limiter)
+ scheduler_addunlock(sched, t_limiter,
+ t->cj->super->timestep_limiter);
+ }
}
#endif
}
@@ -1961,6 +2167,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++) {
@@ -1969,7 +2176,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. */
@@ -1988,6 +2195,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++) {
@@ -1995,7 +2203,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. */
diff --git a/src/engine_marktasks.c b/src/engine_marktasks.c
index 9c7a783c2547899816842cf9a05163e75d329aa8..3a26dbb2f47f9503aa0b93fa28d679f5eebaeede 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,10 @@ 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);
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index b4dc25495ab5be3c2e9c5ba0153e748a344f050f..2b1d19bc916889a5cfdc40b1357f1e3dfe9388af 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -645,6 +645,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/Default/hydro_debug.h b/src/hydro/Default/hydro_debug.h
index 3be9c9e1760591423edbd218d19b46ddf9aad01e..68367beaee97c285057cb055c1fbdbba5c370085 100644
--- a/src/hydro/Default/hydro_debug.h
+++ b/src/hydro/Default/hydro_debug.h
@@ -25,10 +25,11 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
"h=%.3e, "
- "wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, time_bin=%d\n",
+ "wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.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, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->time_bin);
+ p->h, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->time_bin,
+ p->wakeup);
}
#endif /* SWIFT_DEFAULT_HYDRO_DEBUG_H */
diff --git a/src/hydro/Default/hydro_iact.h b/src/hydro/Default/hydro_iact.h
index 72808874c3fc6b58005d0e3ad450eafea8aa4b4d..85c586a4e921e38296453b71a2a2b9637971c28c 100644
--- a/src/hydro/Default/hydro_iact.h
+++ b/src/hydro/Default/hydro_iact.h
@@ -378,4 +378,28 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->force.v_sig = max(pi->force.v_sig, v_sig);
}
+/**
+ * @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;
+ }
+}
+
#endif /* SWIFT_DEFAULT_HYDRO_IACT_H */
diff --git a/src/hydro/Default/hydro_part.h b/src/hydro/Default/hydro_part.h
index 2a18e03cb533ca860f227a31152ef2058e0dd37d..7230826dc3c7c2a3486001ca9060dd07d55d0931 100644
--- a/src/hydro/Default/hydro_part.h
+++ b/src/hydro/Default/hydro_part.h
@@ -21,6 +21,7 @@
#include "chemistry_struct.h"
#include "cooling_struct.h"
+#include "tracers_struct.h"
/* Extra particle data not needed during the SPH loops over neighbours. */
struct xpart {
@@ -40,6 +41,9 @@ struct xpart {
/* Additional data used to record cooling information */
struct cooling_xpart_data cooling_data;
+ /* Additional data used by the tracers */
+ struct tracers_xpart_data tracers_data;
+
float u_full;
/* Old density. */
@@ -132,6 +136,9 @@ struct part {
/* Particle time-bin */
timebin_t time_bin;
+ /* Need waking-up ? */
+ char wakeup;
+
#ifdef SWIFT_DEBUG_CHECKS
/* Time of the last drift */
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 77c4b5dd61b7ff9b84242df61f503fb42b472952..e28091815426e01d36ac0902340cd637969a91d9 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/GizmoMFM/hydro.h b/src/hydro/GizmoMFM/hydro.h
index b00a3578d02f492050c328af49a6108d566e9204..1ab1c1404f54450ddff8d95b51fdf3970daf7377 100644
--- a/src/hydro/GizmoMFM/hydro.h
+++ b/src/hydro/GizmoMFM/hydro.h
@@ -137,6 +137,9 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
p->conserved.momentum[2] * p->v[2]);
#endif
+ p->time_bin = 0;
+ p->wakeup = time_bin_not_awake;
+
/* initialize the particle velocity based on the primitive fluid velocity */
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
diff --git a/src/hydro/GizmoMFM/hydro_debug.h b/src/hydro/GizmoMFM/hydro_debug.h
index e8b0914bd3cf6a99210399c6fc654e526319009f..e3c9f793aec92c7bfa2527143e6ad771c3897a09 100644
--- a/src/hydro/GizmoMFM/hydro_debug.h
+++ b/src/hydro/GizmoMFM/hydro_debug.h
@@ -27,6 +27,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"a=[%.3e,%.3e,%.3e], "
"h=%.3e, "
"time_bin=%d, "
+ "wakeup=%d, "
"rho=%.3e, "
"P=%.3e, "
"gradients={"
@@ -51,7 +52,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"wcount_dh=%.3e, "
"wcount=%.3e}\n",
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->a_hydro[0],
- p->a_hydro[1], p->a_hydro[2], p->h, p->time_bin, p->rho, p->P,
+ p->a_hydro[1], p->a_hydro[2], p->h, p->time_bin, p->wakeup, p->rho, p->P,
p->gradients.rho[0], p->gradients.rho[1], p->gradients.rho[2],
p->gradients.v[0][0], p->gradients.v[0][1], p->gradients.v[0][2],
p->gradients.v[1][0], p->gradients.v[1][1], p->gradients.v[1][2],
diff --git a/src/hydro/GizmoMFM/hydro_iact.h b/src/hydro/GizmoMFM/hydro_iact.h
index 38a97cbea39c1ed5c6926c911941e655e52362aa..09d4c7c70ee2bae8a31d10cb4a568c4627c7b3cd 100644
--- a/src/hydro/GizmoMFM/hydro_iact.h
+++ b/src/hydro/GizmoMFM/hydro_iact.h
@@ -486,4 +486,29 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0, a, H);
}
+/**
+ * @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->timestepvars.vmax >
+ const_limiter_max_v_sig_ratio * pj->timestepvars.vmax) {
+
+ pj->wakeup = time_bin_awake;
+ }
+}
+
#endif /* SWIFT_GIZMO_MFM_HYDRO_IACT_H */
diff --git a/src/hydro/GizmoMFM/hydro_part.h b/src/hydro/GizmoMFM/hydro_part.h
index 0055d7d86a35746a8ba90015b3a6986f8ddb5f9f..a05cae18aaf18feb80f7a4ec383434eadece8a41 100644
--- a/src/hydro/GizmoMFM/hydro_part.h
+++ b/src/hydro/GizmoMFM/hydro_part.h
@@ -21,6 +21,7 @@
#include "chemistry_struct.h"
#include "cooling_struct.h"
+#include "tracers_struct.h"
/* Extra particle data not needed during the computation. */
struct xpart {
@@ -40,6 +41,9 @@ struct xpart {
/* Additional data used to record cooling information */
struct cooling_xpart_data cooling_data;
+ /* Additional data used by the tracers */
+ struct tracers_xpart_data tracers_data;
+
} SWIFT_STRUCT_ALIGN;
/* Data of a single particle. */
@@ -187,6 +191,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/GizmoMFV/hydro.h b/src/hydro/GizmoMFV/hydro.h
index 284b67b3b62cd7c6b75de192b299c0c48d170a05..f4e2b829769a58a4896516907317d02c936f2d65 100644
--- a/src/hydro/GizmoMFV/hydro.h
+++ b/src/hydro/GizmoMFV/hydro.h
@@ -121,6 +121,9 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
const float mass = p->conserved.mass;
+ p->time_bin = 0;
+ p->wakeup = time_bin_not_awake;
+
p->primitives.v[0] = p->v[0];
p->primitives.v[1] = p->v[1];
p->primitives.v[2] = p->v[2];
diff --git a/src/hydro/GizmoMFV/hydro_debug.h b/src/hydro/GizmoMFV/hydro_debug.h
index 8af3f824666529efad833c3bd520ace779718449..181bd6f82d547803c7303bd19be11cf66dc3a8a8 100644
--- a/src/hydro/GizmoMFV/hydro_debug.h
+++ b/src/hydro/GizmoMFV/hydro_debug.h
@@ -27,6 +27,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"a=[%.3e,%.3e,%.3e], "
"h=%.3e, "
"time_bin=%d, "
+ "wakeup=%d, "
"primitives={"
"v=[%.3e,%.3e,%.3e], "
"rho=%.3e, "
@@ -53,9 +54,9 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"wcount_dh=%.3e, "
"wcount=%.3e}\n",
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->a_hydro[0],
- p->a_hydro[1], p->a_hydro[2], p->h, p->time_bin, p->primitives.v[0],
- p->primitives.v[1], p->primitives.v[2], p->primitives.rho,
- p->primitives.P, p->primitives.gradients.rho[0],
+ p->a_hydro[1], p->a_hydro[2], p->h, p->time_bin, p->wakeup,
+ p->primitives.v[0], p->primitives.v[1], p->primitives.v[2],
+ p->primitives.rho, p->primitives.P, p->primitives.gradients.rho[0],
p->primitives.gradients.rho[1], p->primitives.gradients.rho[2],
p->primitives.gradients.v[0][0], p->primitives.gradients.v[0][1],
p->primitives.gradients.v[0][2], p->primitives.gradients.v[1][0],
diff --git a/src/hydro/GizmoMFV/hydro_iact.h b/src/hydro/GizmoMFV/hydro_iact.h
index 2f73e67ea2fdcecc527de8b1af0d15731f967b9b..d882549f8c55018419a2e1730d2ac099bbe1f5ee 100644
--- a/src/hydro/GizmoMFV/hydro_iact.h
+++ b/src/hydro/GizmoMFV/hydro_iact.h
@@ -501,4 +501,29 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0, a, H);
}
+/**
+ * @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->timestepvars.vmax >
+ const_limiter_max_v_sig_ratio * pj->timestepvars.vmax) {
+
+ pj->wakeup = time_bin_awake;
+ }
+}
+
#endif /* SWIFT_GIZMO_MFV_HYDRO_IACT_H */
diff --git a/src/hydro/GizmoMFV/hydro_part.h b/src/hydro/GizmoMFV/hydro_part.h
index 6248ddb11daf39a65be9a57fe51e40386ecda50b..8794b597712963e962cc23c796e9769efd4ea620 100644
--- a/src/hydro/GizmoMFV/hydro_part.h
+++ b/src/hydro/GizmoMFV/hydro_part.h
@@ -21,6 +21,7 @@
#include "chemistry_struct.h"
#include "cooling_struct.h"
+#include "tracers_struct.h"
/* Extra particle data not needed during the computation. */
struct xpart {
@@ -40,6 +41,9 @@ struct xpart {
/* Additional data used to record cooling information */
struct cooling_xpart_data cooling_data;
+ /* Additional data used by the tracers */
+ struct tracers_xpart_data tracers_data;
+
} SWIFT_STRUCT_ALIGN;
/* Data of a single particle. */
@@ -198,6 +202,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/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index d638c168f23c95dc3010f838846f4dfc0522bee5..524774435d03a6d808c4535a6c54b68ad16bcb66 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -740,6 +740,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/Minimal/hydro_debug.h b/src/hydro/Minimal/hydro_debug.h
index 73ffc26b8acf687a5445591ddccd72ea8e8fa8ae..3fadd05f9b93e53f1855c5daa7727d272ffe0fa5 100644
--- a/src/hydro/Minimal/hydro_debug.h
+++ b/src/hydro/Minimal/hydro_debug.h
@@ -41,12 +41,12 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"v_full=[%.3g, %.3g, %.3g], a=[%.3g, %.3g, %.3g], \n "
"m=%.3g, u=%.3g, du/dt=%.3g, P=%.3g, c_s=%.3g, \n "
"v_sig=%.3g, h=%.3g, dh/dt=%.3g, wcount=%.3g, rho=%.3g, \n "
- "dh_drho=%.3g, time_bin=%d \n",
+ "dh_drho=%.3g, 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->mass, p->u, p->u_dt, hydro_get_comoving_pressure(p),
p->force.soundspeed, p->force.v_sig, p->h, p->force.h_dt,
- p->density.wcount, p->rho, p->density.rho_dh, p->time_bin);
+ p->density.wcount, p->rho, p->density.rho_dh, p->time_bin, p->wakeup);
}
#endif /* SWIFT_MINIMAL_HYDRO_DEBUG_H */
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index b29f44588c2e13bb5b7c5c9cd5297205557c3fc9..7fc7a3c67f6c832d70109319ad964e25df30ff4e 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -424,4 +424,28 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->force.v_sig = max(pi->force.v_sig, v_sig);
}
+/**
+ * @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;
+ }
+}
+
#endif /* SWIFT_MINIMAL_HYDRO_IACT_H */
diff --git a/src/hydro/Minimal/hydro_part.h b/src/hydro/Minimal/hydro_part.h
index 1d14a94f2d91bf259df54c875a32bf3072ad33b6..80e472194e6a008859fa7e7fde9c79df6611142b 100644
--- a/src/hydro/Minimal/hydro_part.h
+++ b/src/hydro/Minimal/hydro_part.h
@@ -34,6 +34,7 @@
#include "chemistry_struct.h"
#include "cooling_struct.h"
+#include "tracers_struct.h"
/**
* @brief Particle fields not needed during the SPH loops over neighbours.
@@ -62,6 +63,9 @@ struct xpart {
/*! Additional data used to record cooling information */
struct cooling_xpart_data cooling_data;
+ /* Additional data used by the tracers */
+ struct tracers_xpart_data tracers_data;
+
} SWIFT_STRUCT_ALIGN;
/**
@@ -168,6 +172,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/Planetary/hydro.h b/src/hydro/Planetary/hydro.h
index 957e96dcf391b9027016926a969b28366590664f..ed7aa6b89d2b50ab2e00cedb0b3ef6779689feb1 100644
--- a/src/hydro/Planetary/hydro.h
+++ b/src/hydro/Planetary/hydro.h
@@ -735,6 +735,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/Planetary/hydro_debug.h b/src/hydro/Planetary/hydro_debug.h
index 74261f3b49e2881af1c403013005560efa53a7f1..306f7526404599a051f83dc1b61886ed2aa5b69e 100644
--- a/src/hydro/Planetary/hydro_debug.h
+++ b/src/hydro/Planetary/hydro_debug.h
@@ -42,12 +42,13 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"v_full=[%.3g, %.3g, %.3g], a=[%.3g, %.3g, %.3g], \n "
"m=%.3g, u=%.3g, du/dt=%.3g, P=%.3g, c_s=%.3g, \n "
"v_sig=%.3g, h=%.3g, dh/dt=%.3g, wcount=%.3g, rho=%.3g, \n "
- "dh_drho=%.3g, time_bin=%d, mat_id=%d \n",
+ "dh_drho=%.3g, time_bin=%d, wakeup=%d mat_id=%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->mass, p->u, p->u_dt, hydro_get_comoving_pressure(p),
p->force.soundspeed, p->force.v_sig, p->h, p->force.h_dt,
- p->density.wcount, p->rho, p->density.rho_dh, p->time_bin, p->mat_id);
+ p->density.wcount, p->rho, p->density.rho_dh, p->time_bin, p->wakeup,
+ p->mat_id);
}
#endif /* SWIFT_PLANETARY_HYDRO_DEBUG_H */
diff --git a/src/hydro/Planetary/hydro_iact.h b/src/hydro/Planetary/hydro_iact.h
index 19ee002b85c1b0bc8ed621a029059cd02c5e670f..afebb6a406bd310f38d51dcb32fc25da6b2674b5 100644
--- a/src/hydro/Planetary/hydro_iact.h
+++ b/src/hydro/Planetary/hydro_iact.h
@@ -346,4 +346,28 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->force.v_sig = max(pi->force.v_sig, v_sig);
}
+/**
+ * @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;
+ }
+}
+
#endif /* SWIFT_PLANETARY_HYDRO_IACT_H */
diff --git a/src/hydro/Planetary/hydro_part.h b/src/hydro/Planetary/hydro_part.h
index 4087cef62e873231a556f82869a7f6d848c8d72c..1955366da7265c4c40922d1e7290bc9128641600 100644
--- a/src/hydro/Planetary/hydro_part.h
+++ b/src/hydro/Planetary/hydro_part.h
@@ -36,6 +36,7 @@
#include "chemistry_struct.h"
#include "cooling_struct.h"
#include "equation_of_state.h" // For enum material_id
+#include "tracers_struct.h"
/**
* @brief Particle fields not needed during the SPH loops over neighbours.
@@ -64,6 +65,9 @@ struct xpart {
/*! Additional data used to record cooling information */
struct cooling_xpart_data cooling_data;
+ /* Additional data used by the tracers */
+ struct tracers_xpart_data tracers_data;
+
} SWIFT_STRUCT_ALIGN;
/**
@@ -173,6 +177,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/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h
index 8dd43cd72968f89cfc818342d618688f2f39cbd3..400a84915b700464b9b86f74400ba578b4efa446 100644
--- a/src/hydro/PressureEnergy/hydro.h
+++ b/src/hydro/PressureEnergy/hydro.h
@@ -763,6 +763,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];
@@ -802,4 +803,4 @@ hydro_set_init_internal_energy(struct part *p, float u_init) {
__attribute__((always_inline)) INLINE static void hydro_remove_part(
const struct part *p, const struct xpart *xp) {}
-#endif /* SWIFT_MINIMAL_HYDRO_H */
+#endif /* SWIFT_PRESSURE_ENERGY_HYDRO_H */
diff --git a/src/hydro/PressureEnergy/hydro_debug.h b/src/hydro/PressureEnergy/hydro_debug.h
index 6324167f12726e155eeaa3359be9741aca3a1e42..7ffc370ed4d6abd273fc3d8d5b887f5ccf8e001c 100644
--- a/src/hydro/PressureEnergy/hydro_debug.h
+++ b/src/hydro/PressureEnergy/hydro_debug.h
@@ -32,12 +32,12 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
"h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, \n"
"p_dh=%.3e, p_bar=%.3e \n"
- "time_bin=%d\n",
+ "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->u, p->u_dt, p->force.v_sig, hydro_get_comoving_pressure(p), p->h,
p->force.h_dt, (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho,
- p->density.pressure_bar_dh, p->pressure_bar, p->time_bin);
+ p->density.pressure_bar_dh, p->pressure_bar, p->time_bin, p->wakeup);
}
#endif /* SWIFT_MINIMAL_HYDRO_DEBUG_H */
diff --git a/src/hydro/PressureEnergy/hydro_iact.h b/src/hydro/PressureEnergy/hydro_iact.h
index 4146e61a53dd7ece57e263cb90308e2579aa3930..ae154ea549a52cb24ed7c69453533b7d59b39a85 100644
--- a/src/hydro/PressureEnergy/hydro_iact.h
+++ b/src/hydro/PressureEnergy/hydro_iact.h
@@ -17,8 +17,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
-#ifndef SWIFT_MINIMAL_HYDRO_IACT_H
-#define SWIFT_MINIMAL_HYDRO_IACT_H
+#ifndef SWIFT_PRESSURE_ENERGY_HYDRO_IACT_H
+#define SWIFT_PRESSURE_ENERGY_HYDRO_IACT_H
/**
* @file PressureEnergy/hydro_iact.h
@@ -418,5 +418,28 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Update the signal velocity. */
pi->force.v_sig = max(pi->force.v_sig, v_sig);
}
+/**
+ * @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;
+ }
+}
-#endif /* SWIFT_MINIMAL_HYDRO_IACT_H */
+#endif /* SWIFT_PRESSURE_ENERGY_HYDRO_IACT_H */
diff --git a/src/hydro/PressureEnergy/hydro_io.h b/src/hydro/PressureEnergy/hydro_io.h
index 06762c6124c2c726c4e687980455ab956a5fa79e..701c12283bf77acef4af77598f57705a2b364fa1 100644
--- a/src/hydro/PressureEnergy/hydro_io.h
+++ b/src/hydro/PressureEnergy/hydro_io.h
@@ -17,8 +17,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
-#ifndef SWIFT_PRESSURE_ENERGY_HYDRO_IACT_H
-#define SWIFT_PRESSURE_ENERGY_HYDRO_IACT_H
+#ifndef SWIFT_PRESSURE_ENERGY_HYDRO_IO_H
+#define SWIFT_PRESSURE_ENERGY_HYDRO_IO_H
/**
* @file PressureEnergy/hydro_io.h
* @brief P-U implementation of SPH (i/o routines)
diff --git a/src/hydro/PressureEnergy/hydro_part.h b/src/hydro/PressureEnergy/hydro_part.h
index bc7d14b612556dc722ecca67dd6ce823192e00f0..218fbf5dc17559b07974b68e42f69f4e7a0e8e3b 100644
--- a/src/hydro/PressureEnergy/hydro_part.h
+++ b/src/hydro/PressureEnergy/hydro_part.h
@@ -33,6 +33,7 @@
#include "chemistry_struct.h"
#include "cooling_struct.h"
+#include "tracers_struct.h"
/**
* @brief Particle fields not needed during the SPH loops over neighbours.
@@ -61,6 +62,9 @@ struct xpart {
/*! Additional data used to record cooling information */
struct cooling_xpart_data cooling_data;
+ /* Additional data used by the tracers */
+ struct tracers_xpart_data tracers_data;
+
} SWIFT_STRUCT_ALIGN;
/**
@@ -168,6 +172,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/PressureEnergyMorrisMonaghanAV/hydro.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
index deb013579fd33340236d3dd5817021fd100c0fcb..7ef55b86c24972f8f287273441da99f26285c531 100644
--- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
@@ -50,22 +50,26 @@
#include <float.h>
/**
- * @brief Returns the comoving internal energy of a particle
+ * @brief Returns the comoving internal energy of a particle at the last
+ * time the particle was kicked.
*
* For implementations where the main thermodynamic variable
* is not internal energy, this function computes the internal
* energy from the thermodynamic variable.
*
* @param p The particle of interest
+ * @param xp The extended data of the particle of interest.
*/
__attribute__((always_inline)) INLINE static float
-hydro_get_comoving_internal_energy(const struct part *restrict p) {
+hydro_get_comoving_internal_energy(const struct part *restrict p,
+ const struct xpart *restrict xp) {
- return p->u;
+ return xp->u_full;
}
/**
- * @brief Returns the physical internal energy of a particle
+ * @brief Returns the physical internal energy of a particle at the last
+ * time the particle was kicked.
*
* For implementations where the main thermodynamic variable
* is not internal energy, this function computes the internal
@@ -73,13 +77,15 @@ hydro_get_comoving_internal_energy(const struct part *restrict p) {
* physical coordinates.
*
* @param p The particle of interest.
+ * @param xp The extended data of the particle of interest.
* @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static float
hydro_get_physical_internal_energy(const struct part *restrict p,
+ const struct xpart *restrict xp,
const struct cosmology *cosmo) {
- return p->u * cosmo->a_factor_internal_energy;
+ return xp->u_full * cosmo->a_factor_internal_energy;
}
/**
@@ -734,6 +740,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/PressureEnergyMorrisMonaghanAV/hydro_debug.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h
index ead5fcc0c842d8018f784a1084941bdb9ebcb6ca..d0cd5367f94cd90f36cc2b738a63c7963adbd445 100644
--- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h
@@ -36,12 +36,13 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
"h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, \n"
"p_dh=%.3e, p_bar=%.3e \n"
- "time_bin=%d, alpha=%.3e\n",
+ "time_bin=%d, wakeup=%d alpha=%.3e\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->u, p->u_dt, p->force.v_sig, hydro_get_comoving_pressure(p), p->h,
p->force.h_dt, (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho,
- p->density.pressure_bar_dh, p->pressure_bar, p->time_bin, p->alpha);
+ p->density.pressure_bar_dh, p->pressure_bar, p->time_bin, p->wakeup,
+ p->alpha);
}
#endif /* SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_DEBUG_H */
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h
index 747fca714ce20d9c2b018e14ac24a6492c51a75f..69da511c7544a71ef381a0889c8b56c80d5211f1 100644
--- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h
@@ -424,4 +424,28 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->force.v_sig = max(pi->force.v_sig, v_sig);
}
+/**
+ * @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;
+ }
+}
+
#endif /* SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_IACT_H */
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h
index 1600679bc2e840d0b3b958531c279f5f29293b48..71662f14c61c92d65bcf493b6f5a43b8172e3697 100644
--- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h
@@ -69,12 +69,6 @@ INLINE static void hydro_read_particles(struct part* parts,
UNIT_CONV_DENSITY, parts, rho);
}
-INLINE static void convert_u(const struct engine* e, const struct part* p,
- const struct xpart* xp, float* ret) {
-
- ret[0] = hydro_get_comoving_internal_energy(p);
-}
-
INLINE static void convert_S(const struct engine* e, const struct part* p,
const struct xpart* xp, float* ret) {
@@ -170,9 +164,8 @@ INLINE static void hydro_write_particles(const struct part* parts,
io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
parts, h);
- list[4] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
- UNIT_CONV_ENERGY_PER_UNIT_MASS,
- parts, xparts, convert_u);
+ list[4] = io_make_output_field("InternalEnergy", FLOAT, 1,
+ UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
UNIT_CONV_NO_UNITS, parts, id);
list[6] =
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h
index da6391236811e2a907281c3db05462bb57602fe0..d66249ea179a830cedbd3c3f165ca5012fd18862 100644
--- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h
@@ -34,6 +34,7 @@
#include "chemistry_struct.h"
#include "cooling_struct.h"
+#include "tracers_struct.h"
/**
* @brief Particle fields not needed during the SPH loops over neighbours.
@@ -62,6 +63,9 @@ struct xpart {
/*! Additional data used to record cooling information */
struct cooling_xpart_data cooling_data;
+ /* Additional data used by the tracers */
+ struct tracers_xpart_data tracers_data;
+
} SWIFT_STRUCT_ALIGN;
/**
@@ -172,6 +176,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/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index 38e0f66fe7ecc1b6497717c9754bc36cd10a66f7..2e8d2d5db615f239bf5c3567e7beb155eab5cb38 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -730,6 +730,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;
p->rho_bar = 0.f;
p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
xp->v_full[0] = p->v[0];
diff --git a/src/hydro/PressureEntropy/hydro_debug.h b/src/hydro/PressureEntropy/hydro_debug.h
index 14d69bb650ff1bbd49394c0ca2f6256ad0cb188d..2163b70b94dde4e88f010d962358dccbde7960a3 100644
--- a/src/hydro/PressureEntropy/hydro_debug.h
+++ b/src/hydro/PressureEntropy/hydro_debug.h
@@ -36,14 +36,14 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
"h=%.3e, wcount=%.3f, wcount_dh=%.3e, m=%.3e, dh_drho=%.3e, rho=%.3e, "
"rho_bar=%.3e, P=%.3e, dP_dh=%.3e, P_over_rho2=%.3e, S=%.3e, S^1/g=%.3e, "
- "dS/dt=%.3e,\nc=%.3e v_sig=%e dh/dt=%.3e time_bin=%d\n",
+ "dS/dt=%.3e,\nc=%.3e 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, p->rho_bar, hydro_get_comoving_pressure(p),
p->density.pressure_dh, p->force.P_over_rho2, p->entropy,
p->entropy_one_over_gamma, p->entropy_dt, p->force.soundspeed,
- 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_PRESSURE_ENTROPY_HYDRO_DEBUG_H */
diff --git a/src/hydro/PressureEntropy/hydro_iact.h b/src/hydro/PressureEntropy/hydro_iact.h
index a018b39a99be5ed691485d93bd8dfd1735378bda..19279adec1f37117cf985e63a18a681ceee4f973 100644
--- a/src/hydro/PressureEntropy/hydro_iact.h
+++ b/src/hydro/PressureEntropy/hydro_iact.h
@@ -402,4 +402,28 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->entropy_dt += mj * visc_term * r_inv * dvdr;
}
+/**
+ * @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;
+ }
+}
+
#endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_IACT_H */
diff --git a/src/hydro/PressureEntropy/hydro_part.h b/src/hydro/PressureEntropy/hydro_part.h
index fb8424d66196b7013866acef6bec6ec9889a3353..a404a897b06ddc0777a493e2ecfd28b68e15defe 100644
--- a/src/hydro/PressureEntropy/hydro_part.h
+++ b/src/hydro/PressureEntropy/hydro_part.h
@@ -32,6 +32,7 @@
#include "chemistry_struct.h"
#include "cooling_struct.h"
+#include "tracers_struct.h"
/* Extra particle data not needed during the SPH loops over neighbours. */
struct xpart {
@@ -54,6 +55,9 @@ struct xpart {
/*! Additional data used to record cooling information */
struct cooling_xpart_data cooling_data;
+ /* Additional data used by the tracers */
+ struct tracers_xpart_data tracers_data;
+
} SWIFT_STRUCT_ALIGN;
/* Data of a single particle. */
@@ -148,6 +152,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/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
index 446219104dffb2939877ae2a7c782e66af153213..b0f3207dfce69ca79899b1134740d035d47251d1 100644
--- a/src/hydro/Shadowswift/hydro.h
+++ b/src/hydro/Shadowswift/hydro.h
@@ -103,6 +103,9 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
const float mass = p->conserved.mass;
+ p->time_bin = 0;
+ p->wakeup = time_bin_not_awake;
+
p->primitives.v[0] = p->v[0];
p->primitives.v[1] = p->v[1];
p->primitives.v[2] = p->v[2];
diff --git a/src/hydro/Shadowswift/hydro_debug.h b/src/hydro/Shadowswift/hydro_debug.h
index 7cd7f89c8112ebcf1930c5ca52cb389139191975..8ff85d62fc7d58d53220b1f77a7afb44c00c33b0 100644
--- a/src/hydro/Shadowswift/hydro_debug.h
+++ b/src/hydro/Shadowswift/hydro_debug.h
@@ -23,6 +23,8 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"x=[%.16e,%.16e,%.16e], "
"v=[%.3e,%.3e,%.3e], "
"a=[%.3e,%.3e,%.3e], "
+ "time_bin=%d, "
+ "wakeup=%d, "
"h=%.3e, "
"primitives={"
"v=[%.3e,%.3e,%.3e], "
@@ -47,9 +49,9 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"wcount_dh=%.3e, "
"wcount=%.3e}",
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->a_hydro[0],
- p->a_hydro[1], p->a_hydro[2], p->h, p->primitives.v[0],
- p->primitives.v[1], p->primitives.v[2], p->primitives.rho,
- p->primitives.P, p->primitives.gradients.rho[0],
+ p->a_hydro[1], p->a_hydro[2], p->time_bin, p->wakeup, p->h,
+ p->primitives.v[0], p->primitives.v[1], p->primitives.v[2],
+ p->primitives.rho, p->primitives.P, p->primitives.gradients.rho[0],
p->primitives.gradients.rho[1], p->primitives.gradients.rho[2],
p->primitives.gradients.v[0][0], p->primitives.gradients.v[0][1],
p->primitives.gradients.v[0][2], p->primitives.gradients.v[1][0],
diff --git a/src/hydro/Shadowswift/hydro_iact.h b/src/hydro/Shadowswift/hydro_iact.h
index eda8e3759d9e08dac8073ebed9fb36dd0c5b99f6..791e4c7924df9806fa9150d03c08a543771a7049 100644
--- a/src/hydro/Shadowswift/hydro_iact.h
+++ b/src/hydro/Shadowswift/hydro_iact.h
@@ -342,3 +342,28 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0, a, H);
}
+
+/**
+ * @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->timestepvars.vmax >
+ const_limiter_max_v_sig_ratio * pj->timestepvars.vmax) {
+
+ pj->wakeup = time_bin_awake;
+ }
+}
diff --git a/src/hydro/Shadowswift/hydro_part.h b/src/hydro/Shadowswift/hydro_part.h
index a7cc9daf0839216f098ac05c2267adc60ea11fb0..91ffaa85e5e6e80e7db577ce09363265f73e7f4c 100644
--- a/src/hydro/Shadowswift/hydro_part.h
+++ b/src/hydro/Shadowswift/hydro_part.h
@@ -21,6 +21,7 @@
#include "chemistry_struct.h"
#include "cooling_struct.h"
+#include "tracers_struct.h"
#include "voronoi_cell.h"
/* Extra particle data not needed during the computation. */
@@ -41,6 +42,9 @@ struct xpart {
/* Additional data used to record cooling information */
struct cooling_xpart_data cooling_data;
+ /* Additional data used by the tracers */
+ struct tracers_xpart_data tracers_data;
+
} SWIFT_STRUCT_ALIGN;
/* Data of a single particle. */
@@ -179,6 +183,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 d864eb00367a76c2287f53dbded0fb6feb60ef26..ac226c9856448c82a8d0b8332adad3b1c06ab411 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -268,6 +268,8 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
io_write_attribute_f(h_grpsph, "Viscosity decay length [internal units]",
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..1fb616bc8f434132f745f4e411f99c21b2cc941f 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,144 @@ 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;
+
+ /* Early abort? */
+ if (c->hydro.count == 0) {
+
+ /* Clear the limiter flags. */
+ c->hydro.do_limiter = 0;
+ c->hydro.do_sub_limiter = 0;
+ return;
+ }
+
+ /* 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];
+
+ /* Avoid inhibited particles */
+ if (part_is_inhibited(p, e)) continue;
+
+ /* 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 +2907,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 +2930,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 +2951,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 +2970,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 +3031,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 +3050,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 4e3eb4e29e6cd4a2cd91032d4ee81d203977a59e..ad6af73aec209a19106794636c3b6599baca21e1 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -59,12 +59,43 @@
*/
void scheduler_clear_active(struct scheduler *s) { s->active_count = 0; }
+/**
+ * @brief Increase the space available for unlocks. Only call when
+ * current index == s->size_unlock;
+ */
+static void scheduler_extend_unlocks(struct scheduler *s) {
+
+ /* Allocate the new buffer. */
+ const int size_unlocks_new = s->size_unlocks * 2;
+ struct task **unlocks_new =
+ (struct task **)malloc(sizeof(struct task *) * size_unlocks_new);
+ int *unlock_ind_new = (int *)malloc(sizeof(int) * size_unlocks_new);
+ if (unlocks_new == NULL || unlock_ind_new == NULL)
+ error("Failed to re-allocate unlocks.");
+
+ /* Wait for all writes to the old buffer to complete. */
+ while (s->completed_unlock_writes < s->size_unlocks)
+ ;
+
+ /* Copy the buffers. */
+ memcpy(unlocks_new, s->unlocks, sizeof(struct task *) * s->size_unlocks);
+ memcpy(unlock_ind_new, s->unlock_ind, sizeof(int) * s->size_unlocks);
+ free(s->unlocks);
+ free(s->unlock_ind);
+ s->unlocks = unlocks_new;
+ s->unlock_ind = unlock_ind_new;
+
+ /* Publish the new buffer size. */
+ s->size_unlocks = size_unlocks_new;
+}
+
/**
* @brief Add an unlock_task to the given task.
*
* @param s The #scheduler.
* @param ta The unlocking #task.
* @param tb The #task that will be unlocked.
+
*/
void scheduler_addunlock(struct scheduler *s, struct task *ta,
struct task *tb) {
@@ -77,37 +108,21 @@ void scheduler_addunlock(struct scheduler *s, struct task *ta,
const int ind = atomic_inc(&s->nr_unlocks);
/* Does the buffer need to be grown? */
- if (ind == s->size_unlocks) {
- /* Allocate the new buffer. */
- struct task **unlocks_new;
- int *unlock_ind_new;
- const int size_unlocks_new = s->size_unlocks * 2;
- if ((unlocks_new = (struct task **)malloc(sizeof(struct task *) *
- size_unlocks_new)) == NULL ||
- (unlock_ind_new = (int *)malloc(sizeof(int) * size_unlocks_new)) ==
- NULL)
- error("Failed to re-allocate unlocks.");
-
- /* Wait for all writes to the old buffer to complete. */
- while (s->completed_unlock_writes < ind)
- ;
-
- /* Copy the buffers. */
- memcpy(unlocks_new, s->unlocks, sizeof(struct task *) * ind);
- memcpy(unlock_ind_new, s->unlock_ind, sizeof(int) * ind);
- free(s->unlocks);
- free(s->unlock_ind);
- s->unlocks = unlocks_new;
- s->unlock_ind = unlock_ind_new;
-
- /* Publish the new buffer size. */
- s->size_unlocks = size_unlocks_new;
- }
+ if (ind == s->size_unlocks) scheduler_extend_unlocks(s);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ind > s->size_unlocks * 2)
+ message("unlocks guard enabled: %d / %d", ind, s->size_unlocks);
+#endif
/* Wait for there to actually be space at my index. */
while (ind > s->size_unlocks)
;
+ /* Guard against case when more than (old) s->size_unlocks unlocks
+ * are now pending. */
+ if (ind == s->size_unlocks) scheduler_extend_unlocks(s);
+
/* Write the unlock to the scheduler. */
s->unlocks[ind] = tb;
s->unlock_ind[ind] = ta - s->tasks;
@@ -115,7 +130,7 @@ void scheduler_addunlock(struct scheduler *s, struct task *ta,
}
/**
- * @brief compute the number of same dependencies
+ * @brief compute the number of similar dependencies
*
* @param s The #scheduler
* @param ta The #task
@@ -513,7 +528,7 @@ void scheduler_write_dependencies(struct scheduler *s, int verbose) {
/* Be clean */
free(task_dep);
- if (verbose && s->nodeID == 0)
+ if (verbose)
message("Printing task graph took %.3f %s.",
clocks_from_ticks(getticks() - tic), clocks_getunit());
}
diff --git a/src/space.c b/src/space.c
index 35aaffa66d5b921ec687350588ac91e3a52bb59f..e0cb1501d3307bd56345bcbfc0d2d97426ec33c9 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;
@@ -244,6 +245,8 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->stars.do_sub_sort = 0;
c->grav.do_sub_drift = 0;
c->hydro.do_sub_drift = 0;
+ c->hydro.do_sub_limiter = 0;
+ c->hydro.do_limiter = 0;
c->hydro.ti_end_min = -1;
c->hydro.ti_end_max = -1;
c->grav.ti_end_min = -1;
@@ -272,12 +275,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
}
}
@@ -2707,6 +2712,8 @@ void space_split_recursive(struct space *s, struct cell *c,
cp->stars.do_sub_sort = 0;
cp->grav.do_sub_drift = 0;
cp->hydro.do_sub_drift = 0;
+ cp->hydro.do_sub_limiter = 0;
+ cp->hydro.do_limiter = 0;
#ifdef WITH_MPI
cp->mpi.tag = -1;
#endif // WITH_MPI
@@ -4301,6 +4308,49 @@ 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].time_bin == time_bin_inhibited) continue;
+
+ if (parts[k].wakeup == time_bin_awake)
+ error("Particle still woken up! id=%lld", parts[k].id);
+
+ 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 f16aadc8afb7a2f811c4790688fb849ba1601ce3..4d5695f64c81e710c39fcc460a642a0887856814 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,8 @@ void task_unlock(struct task *t) {
case task_type_drift_part:
case task_type_sort:
+ case task_type_ghost:
+ case task_type_timestep_limiter:
cell_unlocktree(ci);
break;
@@ -462,6 +467,8 @@ int task_lock(struct task *t) {
case task_type_drift_part:
case task_type_sort:
+ case task_type_ghost:
+ case task_type_timestep_limiter:
if (ci->hydro.hold) return 0;
if (cell_locktree(ci) != 0) return 0;
break;
@@ -655,6 +662,9 @@ void task_get_group_name(int type, int subtype, char *cluster) {
case task_subtype_grav:
strcpy(cluster, "Gravity");
break;
+ case task_subtype_limiter:
+ strcpy(cluster, "Timestep_limiter");
+ break;
case task_subtype_stars_density:
strcpy(cluster, "Stars");
break;
diff --git a/src/task.h b/src/task.h
index a6782a6302e2f234f02d2b4e3052a11cb388dc31..100ac225bd5956e8d59d6a197c1257cb3e796ebb 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..cfadc2e62a872a2d2a8a578fe6bb48fd24c5ba29
--- /dev/null
+++ b/src/timestep_limiter.h
@@ -0,0 +1,143 @@
+/*******************************************************************************
+ * 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 struct cosmology *cosmo = e->cosmology;
+ 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) {
+ dt_kick_hydro = -cosmology_get_hydro_kick_factor(cosmo, old_ti_beg,
+ old_ti_beg + old_dti / 2);
+ dt_kick_grav = -cosmology_get_grav_kick_factor(cosmo, old_ti_beg,
+ old_ti_beg + old_dti / 2);
+ dt_kick_therm = -cosmology_get_therm_kick_factor(cosmo, old_ti_beg,
+ old_ti_beg + old_dti / 2);
+ dt_kick_corr = -cosmology_get_corr_kick_factor(cosmo, old_ti_beg,
+ old_ti_beg + old_dti / 2);
+ } 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) {
+ dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, new_ti_beg,
+ new_ti_beg + new_dti / 2);
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, new_ti_beg,
+ new_ti_beg + new_dti / 2);
+ dt_kick_therm = cosmology_get_therm_kick_factor(cosmo, new_ti_beg,
+ new_ti_beg + new_dti / 2);
+ dt_kick_corr = cosmology_get_corr_kick_factor(cosmo, new_ti_beg,
+ new_ti_beg + new_dti / 2);
+ } 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",