Re-generated the code from timestep_limiter2. Code runs until the first rebuild.

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.

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:
   

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:
 

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

examples/SedovBlast_1D/sedov.yml

@@ -11,7 +11,7 @@ TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
   time_end:   5e-2  # The end time of the simulation (in internal units).
   dt_min:     1e-7  # The minimal time-step size of the simulation (in internal units).
-  dt_max:     1e-5  # The maximal time-step size of the simulation (in internal units).
+  dt_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
 
 # Parameters governing the snapshots
 Snapshots:
@@ -21,7 +21,7 @@ Snapshots:
 
 # Parameters governing the conserved quantities statistics
 Statistics:
-  delta_time:          1e-5 # Time between statistics output
+  delta_time:          1e-2 # Time between statistics output
 
 # Parameters for the hydrodynamics scheme
 SPH:

examples/main.c

@@ -153,6 +153,7 @@ int main(int argc, char *argv[]) {
   int with_stars = 0;
   int with_star_formation = 0;
   int with_feedback = 0;
+  int with_limiter = 0;
   int with_fp_exceptions = 0;
   int with_drift_all = 0;
   int with_mpole_reconstruction = 0;
@@ -202,6 +203,8 @@ int main(int argc, char *argv[]) {
       OPT_BOOLEAN('S', "stars", &with_stars, "Run with stars.", NULL, 0, 0),
       OPT_BOOLEAN('x', "velociraptor", &with_structure_finding,
                   "Run with structure finding.", NULL, 0, 0),
+      OPT_BOOLEAN(0, "limiter", &with_limiter, "Run with time-step limiter.",
+                  NULL, 0, 0),
 
       OPT_GROUP("  Control options:\n"),
       OPT_BOOLEAN('a', "pin", &with_aff,
@@ -456,6 +459,7 @@ int main(int argc, char *argv[]) {
   if (with_feedback) error("Can't run with feedback over MPI (yet).");
   if (with_star_formation)
     error("Can't run with star formation over MPI (yet)");
+  if (with_limiter) error("Can't run with time-step limiter over MPI (yet)");
 #endif
 
 #if defined(WITH_MPI) && defined(HAVE_VELOCIRAPTOR)
@@ -878,6 +882,7 @@ int main(int argc, char *argv[]) {
       engine_policies |= engine_policy_external_gravity;
     if (with_cosmology) engine_policies |= engine_policy_cosmology;
     if (with_temperature) engine_policies |= engine_policy_temperature;
+    if (with_limiter) engine_policies |= engine_policy_limiter;
     if (with_cooling) engine_policies |= engine_policy_cooling;
     if (with_stars) engine_policies |= engine_policy_stars;
     if (with_star_formation) engine_policies |= engine_policy_star_formation;

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  \

src/cell.c

@@ -1232,8 +1232,11 @@ void cell_clean_links(struct cell *c, void *data) {
   c->hydro.density = NULL;
   c->hydro.gradient = NULL;
   c->hydro.force = NULL;
+  c->hydro.limiter = NULL;
   c->grav.grav = NULL;
   c->grav.mm = NULL;
+  c->stars.density = NULL;
+  c->stars.feedback = NULL;
 }
 
 /**
@@ -1599,6 +1602,14 @@ void cell_clear_drift_flags(struct cell *c, void *data) {
   c->grav.do_sub_drift = 0;
 }
 
+/**
+ * @brief Clear the limiter flags on the given cell.
+ */
+void cell_clear_limiter_flags(struct cell *c, void *data) {
+  c->hydro.do_limiter = 0;
+  c->hydro.do_sub_limiter = 0;
+}
+
 /**
  * @brief Activate the #part drifts on the given cell.
  */
@@ -1686,6 +1697,34 @@ void cell_activate_drift_spart(struct cell *c, struct scheduler *s) {
   cell_activate_drift_gpart(c, s);
 }
 
+/**
+ * @brief Activate the drifts on the given cell.
+ */
+void cell_activate_limiter(struct cell *c, struct scheduler *s) {
+
+  /* If this cell is already marked for drift, quit early. */
+  if (c->hydro.do_limiter) return;
+
+  /* Mark this cell for drifting. */
+  c->hydro.do_limiter = 1;
+
+  /* Set the do_sub_limiter all the way up and activate the super limiter
+     if this has not yet been done. */
+  if (c == c->super) {
+    scheduler_activate(s, c->timestep_limiter);
+  } else {
+    for (struct cell *parent = c->parent;
+         parent != NULL && !parent->hydro.do_sub_limiter;
+         parent = parent->parent) {
+      parent->hydro.do_sub_limiter = 1;
+      if (parent == c->super) {
+        scheduler_activate(s, parent->timestep_limiter);
+        break;
+      }
+    }
+  }
+}
+
 /**
  * @brief Activate the sorts up a cell hierarchy.
  */
@@ -1816,6 +1855,7 @@ void cell_activate_stars_sorts(struct cell *c, int sid, struct scheduler *s) {
 void cell_activate_subcell_hydro_tasks(struct cell *ci, struct cell *cj,
                                        struct scheduler *s) {
   const struct engine *e = s->space->e;
+  const int with_limiter = (e->policy & engine_policy_limiter);
 
   /* Store the current dx_max and h_max values. */
   ci->hydro.dx_max_part_old = ci->hydro.dx_max_part;
@@ -1849,6 +1889,7 @@ void cell_activate_subcell_hydro_tasks(struct cell *ci, struct cell *cj,
 
       /* We have reached the bottom of the tree: activate drift */
       cell_activate_drift_part(ci, s);
+      if (with_limiter) cell_activate_limiter(ci, s);
     }
   }
 
@@ -2154,6 +2195,12 @@ void cell_activate_subcell_hydro_tasks(struct cell *ci, struct cell *cj,
       if (ci->nodeID == engine_rank) cell_activate_drift_part(ci, s);
       if (cj->nodeID == engine_rank) cell_activate_drift_part(cj, s);
 
+      /* Also activate the time-step limiter */
+      if (ci->nodeID == engine_rank && with_limiter)
+        cell_activate_limiter(ci, s);
+      if (cj->nodeID == engine_rank && with_limiter)
+        cell_activate_limiter(cj, s);
+
       /* Do we need to sort the cells? */
       cell_activate_hydro_sorts(ci, sid, s);
       cell_activate_hydro_sorts(cj, sid, s);
@@ -2718,6 +2765,7 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
 
   struct engine *e = s->space->e;
   const int nodeID = e->nodeID;
+  const int with_limiter = (e->policy & engine_policy_limiter);
   int rebuild = 0;
 
   /* Un-skip the density tasks involved with this cell. */
@@ -2743,6 +2791,7 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
       /* Activate hydro drift */
       if (t->type == task_type_self) {
         if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
+        if (ci->nodeID == nodeID && with_limiter) cell_activate_limiter(ci, s);
       }
 
       /* Set the correct sorting flags and activate hydro drifts */
@@ -2757,6 +2806,10 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
         if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
         if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
 
+        /* Activate the limiter tasks. */
+        if (ci->nodeID == nodeID && with_limiter) cell_activate_limiter(ci, s);
+        if (cj->nodeID == nodeID && with_limiter) cell_activate_limiter(cj, s);
+
         /* Check the sorts and activate them if needed. */
         cell_activate_hydro_sorts(ci, t->flags, s);
         cell_activate_hydro_sorts(cj, t->flags, s);
@@ -2791,7 +2844,11 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
         }
 
         /* If the foreign cell is active, we want its ti_end values. */
-        if (ci_active) scheduler_activate(s, ci->mpi.recv_ti);
+        if (ci_active || with_limiter) scheduler_activate(s, ci->mpi.recv_ti);
+
+        if (with_limiter) scheduler_activate(s, ci->mpi.limiter.recv);
+        if (with_limiter)
+          scheduler_activate_send(s, cj->mpi.limiter.send, ci->nodeID);
 
         /* Is the foreign cell active and will need stuff from us? */
         if (ci_active) {
@@ -2801,6 +2858,7 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
           /* Drift the cell which will be sent; note that not all sent
              particles will be drifted, only those that are needed. */
           cell_activate_drift_part(cj, s);
+          if (with_limiter) cell_activate_limiter(cj, s);
 
           /* If the local cell is also active, more stuff will be needed. */
           if (cj_active) {
@@ -2813,7 +2871,8 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
         }
 
         /* If the local cell is active, send its ti_end values. */
-        if (cj_active) scheduler_activate_send(s, cj->mpi.send_ti, ci_nodeID);
+        if (cj_active || with_limiter)
+          scheduler_activate_send(s, cj->mpi.send_ti, ci_nodeID);
 
       } else if (cj_nodeID != nodeID) {
 
@@ -2830,7 +2889,11 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
         }
 
         /* If the foreign cell is active, we want its ti_end values. */
-        if (cj_active) scheduler_activate(s, cj->mpi.recv_ti);
+        if (cj_active || with_limiter) scheduler_activate(s, cj->mpi.recv_ti);
+
+        if (with_limiter) scheduler_activate(s, cj->mpi.limiter.recv);
+        if (with_limiter)
+          scheduler_activate_send(s, ci->mpi.limiter.send, cj->nodeID);
 
         /* Is the foreign cell active and will need stuff from us? */
         if (cj_active) {
@@ -2840,6 +2903,7 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
           /* Drift the cell which will be sent; note that not all sent
              particles will be drifted, only those that are needed. */
           cell_activate_drift_part(ci, s);
+          if (with_limiter) cell_activate_limiter(ci, s);
 
           /* If the local cell is also active, more stuff will be needed. */
           if (ci_active) {
@@ -2853,7 +2917,8 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
         }
 
         /* If the local cell is active, send its ti_end values. */
-        if (ci_active) scheduler_activate_send(s, ci->mpi.send_ti, cj_nodeID);
+        if (ci_active || with_limiter)
+          scheduler_activate_send(s, ci->mpi.send_ti, cj_nodeID);
       }
 #endif
     }
@@ -2866,6 +2931,8 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
       scheduler_activate(s, l->t);
     for (struct link *l = c->hydro.force; l != NULL; l = l->next)
       scheduler_activate(s, l->t);
+    for (struct link *l = c->hydro.limiter; l != NULL; l = l->next)
+      scheduler_activate(s, l->t);
 
     if (c->hydro.extra_ghost != NULL)
       scheduler_activate(s, c->hydro.extra_ghost);
@@ -2879,7 +2946,6 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
     if (c->hydro.cooling != NULL) scheduler_activate(s, c->hydro.cooling);
     if (c->hydro.star_formation != NULL)
       scheduler_activate(s, c->hydro.star_formation);
-    if (c->sourceterms != NULL) scheduler_activate(s, c->sourceterms);
     if (c->logger != NULL) scheduler_activate(s, c->logger);
   }
 

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);

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

src/engine.c

@@ -113,7 +113,8 @@ const char *engine_policy_names[] = {"none",
                                      "stars",
                                      "structure finding",
                                      "star formation",
-                                     "feedback"};
+                                     "feedback",
+                                     "time-step limiter"};
 
 /** The rank of the engine as a global variable (for messages). */
 int engine_rank;
@@ -1907,6 +1908,10 @@ int engine_estimate_nr_tasks(struct engine *e) {
 #endif
 #endif
   }
+  if (e->policy & engine_policy_limiter) {
+    n1 += 18;
+    n2 += 1;
+  }
   if (e->policy & engine_policy_self_gravity) {
     n1 += 125;
     n2 += 8;
@@ -2585,8 +2590,11 @@ void engine_skip_force_and_kick(struct engine *e) {
     /* Skip everything that updates the particles */
     if (t->type == task_type_drift_part || t->type == task_type_drift_gpart ||
         t->type == task_type_kick1 || t->type == task_type_kick2 ||
-        t->type == task_type_timestep || t->subtype == task_subtype_force ||
-        t->subtype == task_subtype_grav || t->type == task_type_end_force ||
+        t->type == task_type_timestep ||
+        t->type == task_type_timestep_limiter ||
+        t->subtype == task_subtype_force ||
+        t->subtype == task_subtype_limiter || t->subtype == task_subtype_grav ||
+        t->type == task_type_end_force ||
         t->type == task_type_grav_long_range || t->type == task_type_grav_mm ||
         t->type == task_type_grav_down || t->type == task_type_cooling)
       t->skip = 1;
@@ -2594,6 +2602,7 @@ void engine_skip_force_and_kick(struct engine *e) {
 
   /* Run through the cells and clear some flags. */
   space_map_cells_pre(e->s, 1, cell_clear_drift_flags, NULL);
+  space_map_cells_pre(e->s, 1, cell_clear_limiter_flags, NULL);
 }
 
 /**
@@ -2803,6 +2812,11 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
     gravity_exact_force_check(e->s, e, 1e-1);
 #endif
 
+#ifdef SWIFT_DEBUG_CHECKS
+  /* Make sure all woken-up particles have been processed */
+  space_check_limiter(e->s);
+#endif
+
   /* Recover the (integer) end of the next time-step */
   engine_collect_end_of_step(e, 1);
 
@@ -3060,6 +3074,11 @@ void engine_step(struct engine *e) {
     gravity_exact_force_check(e->s, e, 1e-1);
 #endif
 
+#ifdef SWIFT_DEBUG_CHECKS
+  /* Make sure all woken-up particles have been processed */
+  space_check_limiter(e->s);
+#endif
+
   /* Collect information about the next time-step */
   engine_collect_end_of_step(e, 1);
   e->forcerebuild = e->collect_group1.forcerebuild;

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];
 
 /**

src/engine_marktasks.c

@@ -69,6 +69,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
   struct scheduler *s = (struct scheduler *)(((size_t *)extra_data)[2]);
   struct engine *e = (struct engine *)((size_t *)extra_data)[0];
   const int nodeID = e->nodeID;
+  const int with_limiter = e->policy & engine_policy_limiter;
 
   for (int ind = 0; ind < num_elements; ind++) {
 
@@ -90,6 +91,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
         if (cell_is_active_hydro(ci, e)) {
           scheduler_activate(s, t);
           cell_activate_drift_part(ci, s);
+          if (with_limiter) cell_activate_limiter(ci, s);
         }
       }
 
@@ -99,6 +101,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
         if (cell_is_active_hydro(ci, e)) {
           scheduler_activate(s, t);
           cell_activate_subcell_hydro_tasks(ci, NULL, s);
+          if (with_limiter) cell_activate_limiter(ci, s);
         }
       }
 
@@ -111,6 +114,16 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
         if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
       }
 
+      else if (t->type == task_type_self &&
+               t->subtype == task_subtype_limiter) {
+        if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
+      }
+
+      else if (t->type == task_type_sub_self &&
+               t->subtype == task_subtype_limiter) {
+        if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
+      }
+
 #ifdef EXTRA_HYDRO_LOOP
       else if (t_type == task_type_self && t_subtype == task_subtype_gradient) {
         if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
@@ -207,6 +220,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
       /* Only activate tasks that involve a local active cell. */
       if ((t_subtype == task_subtype_density ||
            t_subtype == task_subtype_gradient ||
+           t_subtype == task_subtype_limiter ||
            t_subtype == task_subtype_force) &&
           ((ci_active_hydro && ci_nodeID == nodeID) ||
            (cj_active_hydro && cj_nodeID == nodeID))) {
@@ -226,6 +240,12 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
           if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
           if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
 
+          /* And the limiter */
+          if (ci->nodeID == nodeID && with_limiter)
+            cell_activate_limiter(ci, s);
+          if (cj->nodeID == nodeID && with_limiter)
+            cell_activate_limiter(cj, s);
+
           /* Check the sorts and activate them if needed. */
           cell_activate_hydro_sorts(ci, t->flags, s);
           cell_activate_hydro_sorts(cj, t->flags, s);
@@ -560,6 +580,11 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
         scheduler_activate(s, t);
     }
 
+    /* Time-step limiter? */
+    else if (t->type == task_type_timestep) {
+      if (cell_is_active_hydro(t->ci, e)) scheduler_activate(s, t);
+    }
+
     /* Subgrid tasks */
     else if (t_type == task_type_cooling) {
       if (cell_is_active_hydro(t->ci, e) || cell_is_active_gravity(t->ci, e))

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];

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 */

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 */

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 */

src/hydro_properties.c

@@ -262,6 +262,8 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
                        p->viscosity.alpha_min);
   io_write_attribute_f(h_grpsph, "Viscosity decay length", p->viscosity.length);
   io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
+  io_write_attribute_f(h_grpsph, "Max v_sig ratio (limiter)",
+                       const_limiter_max_v_sig_ratio);
 }
 #endif
 

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