Partially restored old time step criterion. Isolated flux limiter into a...

Partially restored old time step criterion. Isolated flux limiter into a separate file and made it configurable. Removed some superfluous gravity variables.

src/Makefile.am

@@ -86,6 +86,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h
                  hydro/Gizmo/hydro_slope_limiters_cell.h \
                  hydro/Gizmo/hydro_slope_limiters_face.h \
                  hydro/Gizmo/hydro_slope_limiters.h \
+                 hydro/Gizmo/hydro_flux_limiters.h \
                  hydro/Gizmo/hydro_unphysical.h \
                  hydro/Gizmo/hydro_velocities.h \
                  hydro/Shadowswift/hydro_debug.h \

src/const.h

@@ -49,6 +49,9 @@
 #define SLOPE_LIMITER_PER_FACE
 #define SLOPE_LIMITER_CELL_WIDE
 
+/* Types of flux limiter to use (GIZMO_SPH only) */
+#define GIZMO_FLUX_LIMITER
+
 /* Options to control the movement of particles for GIZMO_SPH. */
 /* This option disables particle movement */
 //#define GIZMO_FIX_PARTICLES

src/hydro/Gizmo/hydro.h

@@ -53,12 +53,17 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
   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 =
+  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);
   const float psize = powf(p->geometry.volume / hydro_dimension_unit_sphere,
                            hydro_dimension_inv);
-  return CFL_condition * min(psize / vmax, p->timestepvars.dt_min);
+  float dt = FLT_MAX;
+  if (vmax > 0.) {
+    dt = psize / vmax;
+  }
+  return CFL_condition * dt;
 }
 
 /**
@@ -420,7 +425,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
     struct part* restrict p, struct xpart* restrict xp) {
 
   /* Initialize time step criterion variables */
-  p->timestepvars.dt_min = FLT_MAX;
+  p->timestepvars.vmax = 0.;
 
   /* Set the actual velocity of the particle */
   hydro_velocities_prepare_force(p, xp);
@@ -637,12 +642,6 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
     a_grav[1] = p->gpart->a_grav[1];
     a_grav[2] = p->gpart->a_grav[2];
 
-    /* Store the gravitational acceleration for later use. */
-    /* This is used for the prediction step. */
-    p->gravity.old_a[0] = a_grav[0];
-    p->gravity.old_a[1] = a_grav[1];
-    p->gravity.old_a[2] = a_grav[2];
-
     /* Make sure the gpart knows the mass has changed. */
     p->gpart->mass = p->conserved.mass;
 
@@ -664,7 +663,6 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
     p->conserved.momentum[0] += dt * p->conserved.mass * a_grav[0];
     p->conserved.momentum[1] += dt * p->conserved.mass * a_grav[1];
     p->conserved.momentum[2] += dt * p->conserved.mass * a_grav[2];
-
   }
 
   /* reset fluxes */

src/hydro/Gizmo/hydro_debug.h

@@ -46,7 +46,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       "volume=%.3e, "
       "matrix_E=[[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e]]}, "
       "timestepvars={"
-      "dt_min=%.3e},"
+      "vmax=%.3e},"
       "density={"
       "div_v=%.3e, "
       "wcount_dh=%.3e, "
@@ -75,7 +75,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.dt_min, p->density.div_v, p->density.wcount_dh,
+      p->timestepvars.vmax, 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_flux_limiters.h

+
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * 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_HYDRO_FLUX_LIMITERS_H
+#define SWIFT_HYDRO_FLUX_LIMITERS_H
+
+#ifdef GIZMO_FLUX_LIMITER
+
+#define HYDRO_FLUX_LIMITER_IMPLEMENTATION "GIZMO flux limiter"
+
+/**
+ * @brief Limit the flux between two particles.
+ *
+ * @param flux Unlimited flux between the particles.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void hydro_flux_limiters_apply(
+    float* flux, struct part* pi, struct part* pj) {
+
+  float flux_limit_factor = 1.;
+  const float timefac = max(pi->force.dt, pj->force.dt);
+  const float areafac = max(pi->geometry.Atot, pj->geometry.Atot);
+  const float totfac = timefac * areafac;
+  if (flux[0] * totfac > pi->conserved.mass) {
+    flux_limit_factor = pi->conserved.mass / (flux[0] * totfac);
+  }
+  if (flux[0] * totfac > pj->conserved.mass) {
+    flux_limit_factor =
+        min(pj->conserved.mass / (flux[0] * totfac), flux_limit_factor);
+  }
+  if (flux[4] * totfac > pi->conserved.energy) {
+    flux_limit_factor =
+        min(pi->conserved.energy / (flux[4] * totfac), flux_limit_factor);
+  }
+  if (flux[4] * totfac > pj->conserved.energy) {
+    flux_limit_factor =
+        min(pj->conserved.energy / (flux[4] * totfac), flux_limit_factor);
+  }
+
+  flux[0] *= flux_limit_factor;
+  flux[1] *= flux_limit_factor;
+  flux[2] *= flux_limit_factor;
+  flux[3] *= flux_limit_factor;
+  flux[4] *= flux_limit_factor;
+}
+
+#else
+
+#define HYDRO_FLUX_LIMITER_IMPLEMENTATION "No flux limiter"
+
+/**
+ * @brief Limit the flux between two particles.
+ *
+ * @param flux Unlimited flux between the particles.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void hydro_flux_limiters_apply(
+    float* flux, struct part* pi, struct part* pj) {}
+
+#endif
+
+#endif  // SWIFT_HYDRO_FLUX_LIMITERS_H

src/hydro/Gizmo/hydro_gradients.h

@@ -99,7 +99,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
   float xij_j[3];
   int k;
   float xfac;
-  float a_grav_i[3], a_grav_j[3];
 
   /* perform gradient reconstruction in space and time */
   /* space */
@@ -141,34 +140,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
            pj->primitives.gradients.P[1] * xij_j[1] +
            pj->primitives.gradients.P[2] * xij_j[2];
 
-  a_grav_i[0] = pi->gravity.old_a[0];
-  a_grav_i[1] = pi->gravity.old_a[1];
-  a_grav_i[2] = pi->gravity.old_a[2];
-
-  a_grav_i[0] += pi->gravity.grad_a[0][0] * xij_i[0] +
-                 pi->gravity.grad_a[0][1] * xij_i[1] +
-                 pi->gravity.grad_a[0][2] * xij_i[2];
-  a_grav_i[1] += pi->gravity.grad_a[1][0] * xij_i[0] +
-                 pi->gravity.grad_a[1][1] * xij_i[1] +
-                 pi->gravity.grad_a[1][2] * xij_i[2];
-  a_grav_i[2] += pi->gravity.grad_a[2][0] * xij_i[0] +
-                 pi->gravity.grad_a[2][1] * xij_i[1] +
-                 pi->gravity.grad_a[2][2] * xij_i[2];
-
-  a_grav_j[0] = pj->gravity.old_a[0];
-  a_grav_j[1] = pj->gravity.old_a[1];
-  a_grav_j[2] = pj->gravity.old_a[2];
-
-  a_grav_j[0] += pj->gravity.grad_a[0][0] * xij_j[0] +
-                 pj->gravity.grad_a[0][1] * xij_j[1] +
-                 pj->gravity.grad_a[0][2] * xij_j[2];
-  a_grav_j[1] += pj->gravity.grad_a[1][0] * xij_j[0] +
-                 pj->gravity.grad_a[1][1] * xij_j[1] +
-                 pj->gravity.grad_a[1][2] * xij_j[2];
-  a_grav_j[2] += pj->gravity.grad_a[2][0] * xij_j[0] +
-                 pj->gravity.grad_a[2][1] * xij_j[1] +
-                 pj->gravity.grad_a[2][2] * xij_j[2];
-
   hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);
 
   /* time */
@@ -198,10 +169,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
                hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
                                       pi->primitives.gradients.v[1][1] +
                                       pi->primitives.gradients.v[2][2]));
-
-    dWi[1] += 0.5 * mindt * a_grav_i[0];
-    dWi[2] += 0.5 * mindt * a_grav_i[1];
-    dWi[3] += 0.5 * mindt * a_grav_i[2];
   }
 
   if (Wj[0] > 0.0f) {
@@ -230,10 +197,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
                hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
                                       pj->primitives.gradients.v[1][1] +
                                       pj->primitives.gradients.v[2][2]));
-
-    dWj[1] += 0.5 * mindt * a_grav_j[0];
-    dWj[2] += 0.5 * mindt * a_grav_j[1];
-    dWj[3] += 0.5 * mindt * a_grav_j[2];
   }
 
   Wi[0] += dWi[0];

src/hydro/Gizmo/hydro_gradients_gizmo.h

@@ -45,18 +45,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_init(
   p->primitives.gradients.P[1] = 0.0f;
   p->primitives.gradients.P[2] = 0.0f;
 
-  p->gravity.grad_a[0][0] = 0.0f;
-  p->gravity.grad_a[0][1] = 0.0f;
-  p->gravity.grad_a[0][2] = 0.0f;
-
-  p->gravity.grad_a[1][0] = 0.0f;
-  p->gravity.grad_a[1][1] = 0.0f;
-  p->gravity.grad_a[1][2] = 0.0f;
-
-  p->gravity.grad_a[2][0] = 0.0f;
-  p->gravity.grad_a[2][1] = 0.0f;
-  p->gravity.grad_a[2][2] = 0.0f;
-
   hydro_slope_limit_cell_init(p);
 }
 
@@ -157,35 +145,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
         (Wi[4] - Wj[4]) * wi *
         (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
 
-    pi->gravity.grad_a[0][0] +=
-        (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
-        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
-    pi->gravity.grad_a[0][1] +=
-        (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
-        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
-    pi->gravity.grad_a[0][2] +=
-        (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
-        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
-    pi->gravity.grad_a[1][0] +=
-        (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
-        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
-    pi->gravity.grad_a[1][1] +=
-        (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
-        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
-    pi->gravity.grad_a[1][2] +=
-        (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
-        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
-    pi->gravity.grad_a[2][0] +=
-        (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
-        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
-    pi->gravity.grad_a[2][1] +=
-        (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
-        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
-    pi->gravity.grad_a[2][2] +=
-        (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
-        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
   } else {
     /* The gradient matrix was not well-behaved, switch to SPH gradients */
 
@@ -223,27 +182,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
         wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
     pi->primitives.gradients.P[2] -=
         wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
-
-    pi->gravity.grad_a[0][0] -=
-        wi_dx * dx[0] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
-    pi->gravity.grad_a[0][1] -=
-        wi_dx * dx[1] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
-    pi->gravity.grad_a[0][2] -=
-        wi_dx * dx[2] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
-
-    pi->gravity.grad_a[1][0] -=
-        wi_dx * dx[0] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
-    pi->gravity.grad_a[1][1] -=
-        wi_dx * dx[1] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
-    pi->gravity.grad_a[1][2] -=
-        wi_dx * dx[2] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
-
-    pi->gravity.grad_a[2][0] -=
-        wi_dx * dx[0] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
-    pi->gravity.grad_a[2][1] -=
-        wi_dx * dx[1] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
-    pi->gravity.grad_a[2][2] -=
-        wi_dx * dx[2] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
   }
 
   hydro_slope_limit_cell_collect(pi, pj, r);
@@ -306,35 +244,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
         (Wi[4] - Wj[4]) * wj *
         (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
 
-    pj->gravity.grad_a[0][0] +=
-        (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wj *
-        (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
-    pj->gravity.grad_a[0][1] +=
-        (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wj *
-        (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
-    pj->gravity.grad_a[0][2] +=
-        (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wj *
-        (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
-
-    pj->gravity.grad_a[1][0] +=
-        (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wj *
-        (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
-    pj->gravity.grad_a[1][1] +=
-        (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wj *
-        (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
-    pj->gravity.grad_a[1][2] +=
-        (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wj *
-        (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
-
-    pj->gravity.grad_a[2][0] +=
-        (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wj *
-        (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
-    pj->gravity.grad_a[2][1] +=
-        (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wj *
-        (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
-    pj->gravity.grad_a[2][2] +=
-        (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wj *
-        (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
   } else {
     /* SPH gradients */
 
@@ -371,27 +280,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
         wj_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
     pj->primitives.gradients.P[2] -=
         wj_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
-
-    pj->gravity.grad_a[0][0] -=
-        wj_dx * dx[0] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
-    pj->gravity.grad_a[0][1] -=
-        wj_dx * dx[1] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
-    pj->gravity.grad_a[0][2] -=
-        wj_dx * dx[2] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
-
-    pj->gravity.grad_a[1][0] -=
-        wj_dx * dx[0] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
-    pj->gravity.grad_a[1][1] -=
-        wj_dx * dx[1] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
-    pj->gravity.grad_a[1][2] -=
-        wj_dx * dx[2] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
-
-    pj->gravity.grad_a[2][0] -=
-        wj_dx * dx[0] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
-    pj->gravity.grad_a[2][1] -=
-        wj_dx * dx[1] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
-    pj->gravity.grad_a[2][2] -=
-        wj_dx * dx[2] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
   }
 
   hydro_slope_limit_cell_collect(pj, pi, r);
@@ -493,35 +381,6 @@ hydro_gradients_nonsym_collect(float r2, float *dx, float hi, float hj,
         (Wi[4] - Wj[4]) * wi *
         (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
 
-    pi->gravity.grad_a[0][0] +=
-        (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
-        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
-    pi->gravity.grad_a[0][1] +=
-        (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
-        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
-    pi->gravity.grad_a[0][2] +=
-        (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
-        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
-    pi->gravity.grad_a[1][0] +=
-        (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
-        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
-    pi->gravity.grad_a[1][1] +=
-        (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
-        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
-    pi->gravity.grad_a[1][2] +=
-        (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
-        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
-    pi->gravity.grad_a[2][0] +=
-        (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
-        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
-    pi->gravity.grad_a[2][1] +=
-        (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
-        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
-    pi->gravity.grad_a[2][2] +=
-        (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
-        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
   } else {
     /* Gradient matrix is not well-behaved, switch to SPH gradients */
 
@@ -558,27 +417,6 @@ hydro_gradients_nonsym_collect(float r2, float *dx, float hi, float hj,
         wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
     pi->primitives.gradients.P[2] -=
         wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
-
-    pi->gravity.grad_a[0][0] -=
-        wi_dx * dx[0] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
-    pi->gravity.grad_a[0][1] -=
-        wi_dx * dx[1] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
-    pi->gravity.grad_a[0][2] -=
-        wi_dx * dx[2] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
-
-    pi->gravity.grad_a[1][0] -=
-        wi_dx * dx[0] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
-    pi->gravity.grad_a[1][1] -=
-        wi_dx * dx[1] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
-    pi->gravity.grad_a[1][2] -=
-        wi_dx * dx[2] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
-
-    pi->gravity.grad_a[2][0] -=
-        wi_dx * dx[0] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
-    pi->gravity.grad_a[2][1] -=
-        wi_dx * dx[1] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
-    pi->gravity.grad_a[2][2] -=
-        wi_dx * dx[2] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
   }
 
   hydro_slope_limit_cell_collect(pi, pj, r);
@@ -618,17 +456,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
     p->primitives.gradients.P[1] *= ihdim;
     p->primitives.gradients.P[2] *= ihdim;
 
-    p->gravity.grad_a[0][0] *= ihdim;
-    p->gravity.grad_a[0][1] *= ihdim;
-    p->gravity.grad_a[0][2] *= ihdim;
-
-    p->gravity.grad_a[1][0] *= ihdim;
-    p->gravity.grad_a[1][1] *= ihdim;
-    p->gravity.grad_a[1][2] *= ihdim;
-
-    p->gravity.grad_a[2][0] *= ihdim;
-    p->gravity.grad_a[2][1] *= ihdim;
-    p->gravity.grad_a[2][2] *= ihdim;
   } else {
     const float ihdimp1 = pow_dimension_plus_one(ih);
 
@@ -653,18 +480,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
     p->primitives.gradients.P[0] *= ihdimp1 * volume;
     p->primitives.gradients.P[1] *= ihdimp1 * volume;
     p->primitives.gradients.P[2] *= ihdimp1 * volume;
-
-    p->gravity.grad_a[0][0] *= ihdimp1 * volume;
-    p->gravity.grad_a[0][1] *= ihdimp1 * volume;
-    p->gravity.grad_a[0][2] *= ihdimp1 * volume;
-
-    p->gravity.grad_a[1][0] *= ihdimp1 * volume;
-    p->gravity.grad_a[1][1] *= ihdimp1 * volume;
-    p->gravity.grad_a[1][2] *= ihdimp1 * volume;
-
-    p->gravity.grad_a[2][0] *= ihdimp1 * volume;
-    p->gravity.grad_a[2][1] *= ihdimp1 * volume;
-    p->gravity.grad_a[2][2] *= ihdimp1 * volume;
   }
 
   hydro_slope_limit_cell(p);

src/hydro/Gizmo/hydro_iact.h

@@ -20,6 +20,7 @@
  ******************************************************************************/
 
 #include "adiabatic_index.h"
+#include "hydro_flux_limiters.h"
 #include "hydro_gradients.h"
 #include "riemann.h"
 
@@ -306,7 +307,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, dt_min;
+  float vmax, dvdotdx;
   float vi[3], vj[3], vij[3];
   float Wi[5], Wj[5];
   float dti, dtj, mindt;
@@ -346,16 +347,15 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   }
   dvdotdx = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
             (Wi[3] - Wj[3]) * dx[2];
-  dvdotdx = min(dvdotdx,
-            (vi[0] - vj[0])*dx[0] + (vi[1]-vj[1])*dx[1] + (vi[2]-vj[2])*dx[2]);
+  dvdotdx = min(dvdotdx, (vi[0] - vj[0]) * dx[0] + (vi[1] - vj[1]) * dx[1] +
+                             (vi[2] - vj[2]) * dx[2]);
   if (dvdotdx < 0.) {
     /* the magical factor 3 also appears in Gadget2 */
-    vmax -= 3.*dvdotdx / r;
+    vmax -= 3. * dvdotdx / r;
   }
-  dt_min = r / vmax;
-  pi->timestepvars.dt_min = min(pi->timestepvars.dt_min, dt_min);
+  pi->timestepvars.vmax = max(pi->timestepvars.vmax, vmax);
   if (mode == 1) {
-    pj->timestepvars.dt_min = min(pj->timestepvars.dt_min, dt_min);
+    pj->timestepvars.vmax = max(pj->timestepvars.vmax, vmax);
   }
 
   /* The flux will be exchanged using the smallest time step of the two
@@ -491,33 +491,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   float totflux[5];
   riemann_solve_for_flux(Wi, Wj, n_unit, vij, totflux);
 
-  /* Flux limiter */
-  float flux_limit_factor = 1.;
-  float timefac = max(dti, dtj);
-  float areafac = max(pi->geometry.Atot, pj->geometry.Atot);
-  if (totflux[0] * areafac * timefac > pi->conserved.mass) {
-    flux_limit_factor = pi->conserved.mass / (totflux[0] * areafac * timefac);
-  }
-  if (totflux[0] * areafac * timefac > pj->conserved.mass) {
-    flux_limit_factor =
-        min(pj->conserved.mass / (totflux[0] * areafac * timefac),
-            flux_limit_factor);
-  }
-  if (totflux[4] * areafac * timefac > pi->conserved.energy) {
-    flux_limit_factor =
-        min(pi->conserved.energy / (totflux[4] * areafac * timefac),
-            flux_limit_factor);
-  }
-  if (totflux[4] * areafac * timefac > pj->conserved.energy) {
-    flux_limit_factor =
-        min(pj->conserved.energy / (totflux[4] * areafac * timefac),
-            flux_limit_factor);
-  }
-  totflux[0] *= flux_limit_factor;
-  totflux[1] *= flux_limit_factor;
-  totflux[2] *= flux_limit_factor;
-  totflux[3] *= flux_limit_factor;
-  totflux[4] *= flux_limit_factor;
+  hydro_flux_limiters_apply(totflux, pi, pj);
 
   /* Store mass flux */
   float mflux = Anorm * totflux[0];

src/hydro/Gizmo/hydro_io.h

@@ -18,6 +18,7 @@
  ******************************************************************************/
 
 #include "adiabatic_index.h"