Make the minimal SPH scheme collect the minimal time-bin of the neighbours in...

Make the minimal SPH scheme collect the minimal time-bin of the neighbours in the force and use this when setting the time-step size.

src/const.h

@@ -26,6 +26,9 @@
 /* Time-step limiter maximal difference in signal velocity */
 #define const_limiter_max_v_sig_ratio 4.1f
 
+/* Maximal difference in time-bins between neighbouring particles */
+#define const_timestep_limiter_max_delta_bin 2
+
 /* I/O Constant; this determines the relative tolerance between the value of
  * redshift read from the snapshot, and the value from the parameter file. This
  * current value asserts that they must match within 0.1%. */

src/hydro/Minimal/hydro.h

@@ -622,6 +622,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   p->force.pressure = pressure;
   p->force.soundspeed = soundspeed;
   p->force.balsara = balsara;
+  p->min_ngb_time_bin = num_time_bins + 1;
 }
 
 /**
@@ -760,6 +761,10 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
     struct part *restrict p, const struct cosmology *cosmo) {
 
   p->force.h_dt *= p->h * hydro_dimension_inv;
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (p->min_ngb_time_bin == 0) error("Minimal time-bin of neighbours is 0");
+#endif
 }
 
 /**
@@ -863,6 +868,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
     struct part *restrict p, struct xpart *restrict xp) {
 
   p->time_bin = 0;
+  p->min_ngb_time_bin = num_time_bins + 1;
   p->wakeup = time_bin_not_awake;
   xp->v_full[0] = p->v[0];
   xp->v_full[1] = p->v[1];

src/hydro/Minimal/hydro_iact.h

@@ -306,6 +306,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   /* Update the signal velocity. */
   pi->force.v_sig = max(pi->force.v_sig, v_sig);
   pj->force.v_sig = max(pj->force.v_sig, v_sig);
+
+  /* Update the minimal time-bin */
+  if (pj->time_bin > 0)
+    pi->min_ngb_time_bin = min(pi->min_ngb_time_bin, pj->time_bin);
+  if (pi->time_bin > 0)
+    pj->min_ngb_time_bin = min(pj->min_ngb_time_bin, pi->time_bin);
 }
 
 /**
@@ -424,6 +430,10 @@ __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);
+
+  /* Update the minimal time-bin */
+  if (pj->time_bin > 0)
+    pi->min_ngb_time_bin = min(pi->min_ngb_time_bin, pj->time_bin);
 }
 
 /**

src/hydro/Minimal/hydro_part.h

@@ -183,6 +183,9 @@ struct part {
   /* Need waking-up ? */
   timebin_t wakeup;
 
+  /*! Minimal time-bin across all neighbours */
+  timebin_t min_ngb_time_bin;
+
 #ifdef SWIFT_DEBUG_CHECKS
 
   /* Time of the last drift */

src/timestep.h

@@ -180,7 +180,12 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
   const integertime_t new_dti = make_integer_timestep(
       new_dt, p->time_bin, e->ti_current, e->time_base_inv);
 
-  return new_dti;
+  /* Are we allowed to use this bin given the neighbours? */
+  timebin_t new_bin = get_time_bin(new_dti);
+  new_bin =
+      min(new_bin, p->min_ngb_time_bin + const_timestep_limiter_max_delta_bin);
+
+  return get_integer_timestep(new_bin);
 }
 
 /**