Gizmo time step criterion now properly takes into account the cell size.

src/hydro/Gizmo/hydro.h

@@ -49,25 +49,15 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
   return CFL_condition;
 #endif
 
-  if (p->timestepvars.vmax == 0.) {
-    /* vmax can be zero in vacuum cells that only have vacuum neighbours */
-    /* in this case, the time step should be limited by the maximally
-       allowed time step. Since we do not know what that value is here, we set
-       the time step to a very large value */
-    return FLT_MAX;
-  } else {
-    const float psize = powf(p->geometry.volume / hydro_dimension_unit_sphere,
-                             hydro_dimension_inv);
-    float vrel[3];
-    vrel[0] = p->primitives.v[0] - xp->v_full[0];
-    vrel[1] = p->primitives.v[1] - xp->v_full[1];
-    vrel[2] = p->primitives.v[2] - xp->v_full[2];
-    float vmax =
-        sqrtf(vrel[0] * vrel[0] + vrel[1] * vrel[1] + vrel[2] * vrel[2]) +
-        sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
-    vmax = max(vmax, p->timestepvars.vmax);
-    return CFL_condition * psize / vmax;
-  }
+  float vrel[3];
+  vrel[0] = p->primitives.v[0] - xp->v_full[0];
+  vrel[1] = p->primitives.v[1] - xp->v_full[1];
+  vrel[2] = p->primitives.v[2] - xp->v_full[2];
+  const float vmax =
+      sqrtf(vrel[0] * vrel[0] + vrel[1] * vrel[1] + vrel[2] * vrel[2]) +
+      sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
+  return CFL_condition *
+         min(0.5 * p->timestepvars.rmin / vmax, p->timestepvars.dt_min);
 }
 
 /**
@@ -429,7 +419,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
     struct part* restrict p, struct xpart* restrict xp) {
 
   /* Initialize time step criterion variables */
-  p->timestepvars.vmax = 0.0f;
+  p->timestepvars.dt_min = FLT_MAX;
+  p->timestepvars.rmin = FLT_MAX;
 
   /* Set the actual velocity of the particle */
   hydro_velocities_prepare_force(p, xp);

src/hydro/Gizmo/hydro_debug.h

@@ -46,7 +46,8 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       "volume=%.3e, "
       "matrix_E=[[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e]]}, "
       "timestepvars={"
-      "vmax=%.3e}, "
+      "dt_min=%.3e, "
+      "rmin=%.3e},"
       "density={"
       "div_v=%.3e, "
       "wcount_dh=%.3e, "
@@ -75,7 +76,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       p->geometry.matrix_E[1][0], p->geometry.matrix_E[1][1],
       p->geometry.matrix_E[1][2], p->geometry.matrix_E[2][0],
       p->geometry.matrix_E[2][1], p->geometry.matrix_E[2][2],
-      p->timestepvars.vmax, p->density.div_v, p->density.wcount_dh,
+      p->timestepvars.dt_min, p->timestepvars.rmin, p->density.div_v, p->density.wcount_dh,
       p->density.curl_v[0], p->density.curl_v[1], p->density.curl_v[2],
       p->density.wcount);
 }

src/hydro/Gizmo/hydro_iact.h

@@ -306,7 +306,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   float Bj[3][3];
   float Vi, Vj;
   float xij_i[3], xfac, xijdotdx;
-  float vmax, dvdotdx;
+  float vmax, dvdotdx, dt_min;
   float vi[3], vj[3], vij[3];
   float Wi[5], Wj[5];
   float dti, dtj, mindt;
@@ -349,9 +349,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   if (dvdotdx < 0.) {
     vmax -= dvdotdx / r;
   }
-  pi->timestepvars.vmax = max(pi->timestepvars.vmax, vmax);
+  dt_min = r / vmax;
+  pi->timestepvars.dt_min = min(pi->timestepvars.dt_min, dt_min);
+  pi->timestepvars.rmin = min(pi->timestepvars.rmin, r);
   if (mode == 1) {
-    pj->timestepvars.vmax = max(pj->timestepvars.vmax, vmax);
+    pj->timestepvars.dt_min = min(pj->timestepvars.dt_min, dt_min);
+    pj->timestepvars.rmin = min(pj->timestepvars.rmin, r);
   }
 
   /* The flux will be exchanged using the smallest time step of the two

src/hydro/Gizmo/hydro_part.h

@@ -153,11 +153,14 @@ struct part {
 
   } geometry;
 
-  /* Variables used for timestep calculation (currently not used). */
+  /* Variables used for timestep calculation. */
   struct {
 
-    /* Maximum fluid velocity among all neighbours. */
-    float vmax;
+    /* Minimum non-normalized timestep based on the neighbours. */
+    float dt_min;
+
+    /* Minimum distance between the particle and any of its neighbours. */
+    float rmin;
 
   } timestepvars;