Time-step limiter and time-step synchronization over MPI

examples/main.c

@@ -583,10 +583,6 @@ int main(int argc, char *argv[]) {
 #ifdef WITH_MPI
   if (with_mpole_reconstruction && nr_nodes > 1)
     error("Cannot reconstruct m-poles every step over MPI (yet).");
-  if (with_timestep_limiter)
-    error("Can't run with time-step limiter over MPI (yet)");
-  if (with_timestep_sync)
-    error("Can't run with time-step synchronization over MPI (yet)");
 #ifdef WITH_LOGGER
   error("Can't run with the particle logger over MPI (yet)");
 #endif

src/active.h

@@ -136,9 +136,9 @@ __attribute__((always_inline)) INLINE static int cell_is_active_hydro(
   if (c->hydro.ti_end_min < e->ti_current)
     error(
         "cell in an impossible time-zone! c->ti_end_min=%lld (t=%e) and "
-        "e->ti_current=%lld (t=%e, a=%e)",
+        "e->ti_current=%lld (t=%e, a=%e) c->nodeID=%d",
         c->hydro.ti_end_min, c->hydro.ti_end_min * e->time_base, e->ti_current,
-        e->ti_current * e->time_base, e->cosmology->a);
+        e->ti_current * e->time_base, e->cosmology->a, c->nodeID);
 #endif
 
   return (c->hydro.ti_end_min == e->ti_current);

src/cell.c

@@ -3476,16 +3476,15 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
           }
         }
 
+        /* If the foreign cell is active, we want its particles for the limiter
+         */
+        if (ci_active && with_timestep_limiter)
+          scheduler_activate_recv(s, ci->mpi.recv, task_subtype_limiter);
+
         /* If the foreign cell is active, we want its ti_end values. */
-        if (ci_active || with_timestep_limiter)
+        if (ci_active)
           scheduler_activate_recv(s, ci->mpi.recv, task_subtype_tend_part);
 
-        if (with_timestep_limiter)
-          scheduler_activate_recv(s, ci->mpi.recv, task_subtype_limiter);
-        if (with_timestep_limiter)
-          scheduler_activate_send(s, cj->mpi.send, task_subtype_limiter,
-                                  ci->nodeID);
-
         /* Is the foreign cell active and will need stuff from us? */
         if (ci_active) {
 
@@ -3508,8 +3507,13 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
           }
         }
 
+        /* If the local cell is active, send its particles for the limiting. */
+        if (cj_active && with_timestep_limiter)
+          scheduler_activate_send(s, cj->mpi.send, task_subtype_limiter,
+                                  ci_nodeID);
+
         /* If the local cell is active, send its ti_end values. */
-        if (cj_active || with_timestep_limiter)
+        if (cj_active)
           scheduler_activate_send(s, cj->mpi.send, task_subtype_tend_part,
                                   ci_nodeID);
 
@@ -3544,16 +3548,15 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
           }
         }
 
+        /* If the foreign cell is active, we want its particles for the limiter
+         */
+        if (cj_active && with_timestep_limiter)
+          scheduler_activate_recv(s, cj->mpi.recv, task_subtype_limiter);
+
         /* If the foreign cell is active, we want its ti_end values. */
-        if (cj_active || with_timestep_limiter)
+        if (cj_active)
           scheduler_activate_recv(s, cj->mpi.recv, task_subtype_tend_part);
 
-        if (with_timestep_limiter)
-          scheduler_activate_recv(s, cj->mpi.recv, task_subtype_limiter);
-        if (with_timestep_limiter)
-          scheduler_activate_send(s, ci->mpi.send, task_subtype_limiter,
-                                  cj->nodeID);
-
         /* Is the foreign cell active and will need stuff from us? */
         if (cj_active) {
 
@@ -3577,8 +3580,13 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
           }
         }
 
+        /* If the local cell is active, send its particles for the limiting. */
+        if (ci_active && with_timestep_limiter)
+          scheduler_activate_send(s, ci->mpi.send, task_subtype_limiter,
+                                  cj_nodeID);
+
         /* If the local cell is active, send its ti_end values. */
-        if (ci_active || with_timestep_limiter)
+        if (ci_active)
           scheduler_activate_send(s, ci->mpi.send, task_subtype_tend_part,
                                   cj_nodeID);
 
@@ -5281,12 +5289,16 @@ void cell_check_timesteps(const struct cell *c, const integertime_t ti_current,
   integertime_t ti_end_max = 0;
   integertime_t ti_beg_max = 0;
 
+  int count = 0;
+
   for (int i = 0; i < c->hydro.count; ++i) {
 
     const struct part *p = &c->hydro.parts[i];
     if (p->time_bin == time_bin_inhibited) continue;
     if (p->time_bin == time_bin_not_created) continue;
 
+    ++count;
+
     integertime_t ti_end, ti_beg;
 
     if (p->time_bin <= max_bin) {
@@ -5303,14 +5315,15 @@ void cell_check_timesteps(const struct cell *c, const integertime_t ti_current,
     ti_beg_max = max(ti_beg, ti_beg_max);
   }
 
-  if (c->hydro.count > 0) {
+  /* Only check cells that have at least one non-inhibited particle */
+  if (count > 0) {
 
     if (ti_end_min != c->hydro.ti_end_min)
       error(
           "Non-matching ti_end_min. Cell=%lld true=%lld ti_current=%lld "
           "depth=%d",
           c->hydro.ti_end_min, ti_end_min, ti_current, c->depth);
-    if (ti_end_max != c->hydro.ti_end_max)
+    if (ti_end_max > c->hydro.ti_end_max)
       error(
           "Non-matching ti_end_max. Cell=%lld true=%lld ti_current=%lld "
           "depth=%d",

src/engine.c

@@ -1855,9 +1855,9 @@ void engine_skip_drift(struct engine *e) {
  * @brief Launch the runners.
  *
  * @param e The #engine.
- * @param fof Are we launching the FOF tasks or the regular tasks?
+ * @param call What kind of tasks are we running? (For time analysis)
  */
-void engine_launch(struct engine *e, const int fof) {
+void engine_launch(struct engine *e, const char *call) {
   const ticks tic = getticks();
 
 #ifdef SWIFT_DEBUG_CHECKS
@@ -1883,14 +1883,9 @@ void engine_launch(struct engine *e, const int fof) {
   /* Sit back and wait for the runners to come home. */
   swift_barrier_wait(&e->wait_barrier);
 
-  if (e->verbose) {
-    if (fof)
-      message("(fof) took %.3f %s.", clocks_from_ticks(getticks() - tic),
-              clocks_getunit());
-    else
-      message("(tasks) took %.3f %s.", clocks_from_ticks(getticks() - tic),
-              clocks_getunit());
-  }
+  if (e->verbose)
+    message("(%s) took %.3f %s.", call, clocks_from_ticks(getticks() - tic),
+            clocks_getunit());
 }
 
 /**
@@ -1977,7 +1972,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
 
   /* Now, launch the calculation */
   TIMER_TIC;
-  engine_launch(e, /*fof=*/0);
+  engine_launch(e, "tasks");
   TIMER_TOC(timer_runners);
 
   /* Apply some conversions (e.g. internal energy -> entropy) */
@@ -1991,7 +1986,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
     if (hydro_need_extra_init_loop) {
       engine_marktasks(e);
       engine_skip_force_and_kick(e);
-      engine_launch(e, /*fof=*/0);
+      engine_launch(e, "tasks");
     }
   }
 
@@ -2039,9 +2034,19 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
 
   /* Run the 0th time-step */
   TIMER_TIC2;
-  engine_launch(e, /*fof=*/0);
+  engine_launch(e, "tasks");
   TIMER_TOC2(timer_runners);
 
+  /* Since the time-steps may have changed because of the limiter's
+   * action, we need to communicate the new time-step sizes */
+  if ((e->policy & engine_policy_timestep_sync) ||
+      (e->policy & engine_policy_timestep_limiter)) {
+#ifdef WITH_MPI
+    engine_unskip_timestep_communications(e);
+    engine_launch(e, "timesteps");
+#endif
+  }
+
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   /* Check the accuracy of the gravity calculation */
   if (e->policy & engine_policy_self_gravity)
@@ -2337,9 +2342,19 @@ void engine_step(struct engine *e) {
 
   /* Start all the tasks. */
   TIMER_TIC;
-  engine_launch(e, /*fof=*/0);
+  engine_launch(e, "tasks");
   TIMER_TOC(timer_runners);
 
+  /* Since the time-steps may have changed because of the limiter's
+   * action, we need to communicate the new time-step sizes */
+  if ((e->policy & engine_policy_timestep_sync) ||
+      (e->policy & engine_policy_timestep_limiter)) {
+#ifdef WITH_MPI
+    engine_unskip_timestep_communications(e);
+    engine_launch(e, "timesteps");
+#endif
+  }
+
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   /* Check the accuracy of the gravity calculation */
   if (e->policy & engine_policy_self_gravity)

src/engine.h

@@ -492,6 +492,7 @@ void engine_compute_next_fof_time(struct engine *e);
 void engine_compute_next_statistics_time(struct engine *e);
 void engine_recompute_displacement_constraint(struct engine *e);
 void engine_unskip(struct engine *e);
+void engine_unskip_timestep_communications(struct engine *e);
 void engine_drift_all(struct engine *e, const int drift_mpoles);
 void engine_drift_top_multipoles(struct engine *e);
 void engine_reconstruct_multipoles(struct engine *e);
@@ -523,7 +524,7 @@ void engine_config(int restart, int fof, struct engine *e,
                    int nr_threads, int with_aff, int verbose,
                    const char *restart_file);
 void engine_dump_index(struct engine *e);
-void engine_launch(struct engine *e, const int fof);
+void engine_launch(struct engine *e, const char *call);
 void engine_prepare(struct engine *e);
 void engine_init_particles(struct engine *e, int flag_entropy_ICs,
                            int clean_h_values);

src/engine_fof.c

@@ -127,7 +127,7 @@ void engine_fof(struct engine *e, const int dump_results,
   if (e->verbose) engine_print_task_counts(e);
 
   /* Perform local FOF tasks. */
-  engine_launch(e, /*fof=*/1);
+  engine_launch(e, "fof");
 
   /* Perform FOF search over foreign particles and
    * find groups which require black hole seeding.  */

src/engine_maketasks.c

@@ -133,12 +133,16 @@ void engine_addtasks_send_gravity(struct engine *e, struct cell *ci,
  * @param t_xv The send_xv #task, if it has already been created.
  * @param t_rho The send_rho #task, if it has already been created.
  * @param t_gradient The send_gradient #task, if already created.
- * @param t_ti The recv_ti_end #task, if it has already been created.
+ * @param t_ti The send_ti_end #task, if it has already been created.
+ * @param t_limiter The send_limiter #task, if it has already been created.
+ * @param with_limiter Are we running with the time-step limiter?
+ * @param with_sync Are we running with time-step synchronization?
  */
 void engine_addtasks_send_hydro(struct engine *e, struct cell *ci,
                                 struct cell *cj, struct task *t_xv,
                                 struct task *t_rho, struct task *t_gradient,
-                                struct task *t_ti) {
+                                struct task *t_ti, struct task *t_limiter,
+                                const int with_limiter, const int with_sync) {
 
 #ifdef WITH_MPI
   struct link *l = NULL;
@@ -176,6 +180,11 @@ void engine_addtasks_send_hydro(struct engine *e, struct cell *ci,
       t_ti = scheduler_addtask(s, task_type_send, task_subtype_tend_part,
                                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);
+      }
+
 #ifdef EXTRA_HYDRO_LOOP
 
       scheduler_addunlock(s, t_gradient, ci->hydro.super->hydro.end_force);
@@ -210,6 +219,7 @@ void engine_addtasks_send_hydro(struct engine *e, struct cell *ci,
       scheduler_addunlock(s, ci->hydro.super->hydro.drift, t_xv);
 
       scheduler_addunlock(s, ci->super->timestep, t_ti);
+      if (with_limiter) scheduler_addunlock(s, ci->super->timestep, t_limiter);
     }
 
     /* Add them to the local cell. */
@@ -219,6 +229,7 @@ void engine_addtasks_send_hydro(struct engine *e, struct cell *ci,
     engine_addlink(e, &ci->mpi.send, t_gradient);
 #endif
     engine_addlink(e, &ci->mpi.send, t_ti);
+    if (with_limiter) engine_addlink(e, &ci->mpi.send, t_limiter);
   }
 
   /* Recurse? */
@@ -226,7 +237,8 @@ void engine_addtasks_send_hydro(struct engine *e, struct cell *ci,
     for (int k = 0; k < 8; k++)
       if (ci->progeny[k] != NULL)
         engine_addtasks_send_hydro(e, ci->progeny[k], cj, t_xv, t_rho,
-                                   t_gradient, t_ti);
+                                   t_gradient, t_ti, t_limiter, with_limiter,
+                                   with_sync);
 
 #else
   error("SWIFT was not compiled with MPI support.");
@@ -441,10 +453,15 @@ void engine_addtasks_send_black_holes(struct engine *e, struct cell *ci,
  * @param t_rho The recv_rho #task, if it has already been created.
  * @param t_gradient The recv_gradient #task, if it has already been created.
  * @param t_ti The recv_ti_end #task, if it has already been created.
+ * @param t_limiter The recv_limiter #task, if it has already been created.
+ * @param with_limiter Are we running with the time-step limiter?
+ * @param with_sync Are we running with time-step synchronization?
  */
 void engine_addtasks_recv_hydro(struct engine *e, struct cell *c,
                                 struct task *t_xv, struct task *t_rho,
-                                struct task *t_gradient, struct task *t_ti) {
+                                struct task *t_gradient, struct task *t_ti,
+                                struct task *t_limiter, const int with_limiter,
+                                const int with_sync) {
 
 #ifdef WITH_MPI
   struct scheduler *s = &e->sched;
@@ -458,7 +475,7 @@ void engine_addtasks_recv_hydro(struct engine *e, struct cell *c,
 #ifdef SWIFT_DEBUG_CHECKS
     /* Make sure this cell has a valid tag. */
     if (c->mpi.tag < 0) error("Trying to receive from untagged cell.");
-#endif  // SWIFT_DEBUG_CHECKS
+#endif /* SWIFT_DEBUG_CHECKS */
 
     /* Create the tasks. */
     t_xv = scheduler_addtask(s, task_type_recv, task_subtype_xv, c->mpi.tag, 0,
@@ -472,6 +489,11 @@ void engine_addtasks_recv_hydro(struct engine *e, struct cell *c,
 
     t_ti = scheduler_addtask(s, task_type_recv, task_subtype_tend_part,
                              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);
+    }
   }
 
   if (t_xv != NULL) {
@@ -481,6 +503,7 @@ void engine_addtasks_recv_hydro(struct engine *e, struct cell *c,
     engine_addlink(e, &c->mpi.recv, t_gradient);
 #endif
     engine_addlink(e, &c->mpi.recv, t_ti);
+    if (with_limiter) engine_addlink(e, &c->mpi.recv, t_limiter);
 
     /* Add dependencies. */
     if (c->hydro.sorts != NULL) {
@@ -508,8 +531,15 @@ void engine_addtasks_recv_hydro(struct engine *e, struct cell *c,
     }
 #endif
 
-    /* Make sure the density has been computed before the stars compute theirs.
-     */
+    if (with_limiter) {
+      for (struct link *l = c->hydro.limiter; l != NULL; l = l->next) {
+        scheduler_addunlock(s, t_ti, l->t);
+        scheduler_addunlock(s, t_limiter, l->t);
+      }
+    }
+
+    /* Make sure the gas density has been computed before the
+     * stars compute theirs. */
     for (struct link *l = c->stars.density; l != NULL; l = l->next) {
       scheduler_addunlock(s, t_rho, l->t);
     }
@@ -526,7 +556,7 @@ void engine_addtasks_recv_hydro(struct engine *e, struct cell *c,
     for (int k = 0; k < 8; k++)
       if (c->progeny[k] != NULL)
         engine_addtasks_recv_hydro(e, c->progeny[k], t_xv, t_rho, t_gradient,
-                                   t_ti);
+                                   t_ti, t_limiter, with_limiter, with_sync);
 
 #else
   error("SWIFT was not compiled with MPI support.");
@@ -2914,6 +2944,8 @@ void engine_addtasks_send_mapper(void *map_data, int num_elements,
 
   struct engine *e = (struct engine *)extra_data;
   const int with_star_formation = (e->policy & engine_policy_star_formation);
+  const int with_limiter = (e->policy & engine_policy_timestep_limiter);
+  const int with_sync = (e->policy & engine_policy_timestep_sync);
   struct cell_type_pair *cell_type_pairs = (struct cell_type_pair *)map_data;
 
   for (int k = 0; k < num_elements; k++) {
@@ -2926,7 +2958,8 @@ void engine_addtasks_send_mapper(void *map_data, int num_elements,
     if ((e->policy & engine_policy_hydro) && (type & proxy_cell_type_hydro))
       engine_addtasks_send_hydro(e, ci, cj, /*t_xv=*/NULL,
                                  /*t_rho=*/NULL, /*t_gradient=*/NULL,
-                                 /*t_ti=*/NULL);
+                                 /*t_ti=*/NULL, /*t_limiter=*/NULL,
+                                 with_limiter, with_sync);
 
     /* Add the send tasks for the cells in the proxy that have a stars
      * connection. */
@@ -2958,6 +2991,8 @@ void engine_addtasks_recv_mapper(void *map_data, int num_elements,
 
   struct engine *e = (struct engine *)extra_data;
   const int with_star_formation = (e->policy & engine_policy_star_formation);
+  const int with_limiter = (e->policy & engine_policy_timestep_limiter);
+  const int with_sync = (e->policy & engine_policy_timestep_sync);
   struct cell_type_pair *cell_type_pairs = (struct cell_type_pair *)map_data;
 
   for (int k = 0; k < num_elements; k++) {
@@ -2968,7 +3003,8 @@ void engine_addtasks_recv_mapper(void *map_data, int num_elements,
      * connection. */
     if ((e->policy & engine_policy_hydro) && (type & proxy_cell_type_hydro))
       engine_addtasks_recv_hydro(e, ci, /*t_xv=*/NULL, /*t_rho=*/NULL,
-                                 /*t_gradient=*/NULL, /*t_ti=*/NULL);
+                                 /*t_gradient=*/NULL, /*t_ti=*/NULL,
+                                 /*t_limiter=*/NULL, with_limiter, with_sync);
 
     /* Add the recv tasks for the cells in the proxy that have a stars
      * connection. */

src/engine_marktasks.c

@@ -491,6 +491,11 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
             }
           }
 
+          /* If the foreign cell is active, we want its particles for the
+           * limiter */
+          if (ci_active_hydro && with_timestep_limiter)
+            scheduler_activate_recv(s, ci->mpi.recv, task_subtype_limiter);
+
           /* If the foreign cell is active, we want its ti_end values. */
           if (ci_active_hydro)
             scheduler_activate_recv(s, ci->mpi.recv, task_subtype_tend_part);
@@ -517,6 +522,12 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
             }
           }
 
+          /* If the local cell is active, send its particles for the limiting.
+           */
+          if (cj_active_hydro && with_timestep_limiter)
+            scheduler_activate_send(s, cj->mpi.send, task_subtype_limiter,