Gizmo swift_fixdt works.

src/hydro/Gizmo/hydro.h

@@ -171,6 +171,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
 
   /* Set the physical time step */
   p->force.dt = (p->ti_end - p->ti_begin) * timeBase;
+  /* Set the actual velocity of the particle */
+  p->force.v_full[0] = xp->v_full[0];
+  p->force.v_full[1] = xp->v_full[1];
+  p->force.v_full[2] = xp->v_full[2];
 }
 
 /**
@@ -389,18 +393,26 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
 __attribute__((always_inline)) INLINE static void hydro_end_force(
     struct part* p) {
 
-  return;
   /* Set the hydro acceleration, based on the new momentum and mass */
   /* NOTE: the momentum and mass are only correct for active particles, since
            only active particles have received flux contributions from all their
            neighbours. Since this method is only called for active particles,
            this is indeed the case. */
-  p->a_hydro[0] =
-      (p->conserved.momentum[0] / p->conserved.mass - p->v[0]) / p->force.dt;
-  p->a_hydro[1] =
-      (p->conserved.momentum[1] / p->conserved.mass - p->v[1]) / p->force.dt;
-  p->a_hydro[2] =
-      (p->conserved.momentum[2] / p->conserved.mass - p->v[2]) / p->force.dt;
+  if (p->force.dt) {
+    p->a_hydro[0] =
+        (p->conserved.momentum[0] / p->conserved.mass - p->primitives.v[0]) /
+        p->force.dt;
+    p->a_hydro[1] =
+        (p->conserved.momentum[1] / p->conserved.mass - p->primitives.v[1]) /
+        p->force.dt;
+    p->a_hydro[2] =
+        (p->conserved.momentum[2] / p->conserved.mass - p->primitives.v[2]) /
+        p->force.dt;
+  } else {
+    p->a_hydro[0] = 0.0f;
+    p->a_hydro[1] = 0.0f;
+    p->a_hydro[2] = 0.0f;
+  }
 }
 
 __attribute__((always_inline)) INLINE static void hydro_kick_extra(

src/hydro/Gizmo/hydro_gradients_sph.h

@@ -23,29 +23,6 @@
 __attribute__((always_inline)) INLINE static void
 hydro_gradients_init_density_loop(struct part *p) {
 
-  /* use the old volumes to estimate new primitive variables to be used for the
-     gradient calculation */
-  if (p->conserved.mass) {
-    p->primitives.rho = p->conserved.mass / p->geometry.volume;
-    p->primitives.v[0] = p->conserved.momentum[0] / p->conserved.mass;
-    p->primitives.v[1] = p->conserved.momentum[1] / p->conserved.mass;
-    p->primitives.v[2] = p->conserved.momentum[2] / p->conserved.mass;
-    p->primitives.P =
-        hydro_gamma_minus_one *
-        (p->conserved.energy -
-         0.5 * (p->conserved.momentum[0] * p->conserved.momentum[0] +
-                p->conserved.momentum[1] * p->conserved.momentum[1] +
-                p->conserved.momentum[2] * p->conserved.momentum[2]) /
-             p->conserved.mass) /
-        p->geometry.volume;
-  } else {
-    p->primitives.rho = 0.0f;
-    p->primitives.v[0] = 0.0f;
-    p->primitives.v[1] = 0.0f;
-    p->primitives.v[2] = 0.0f;
-    p->primitives.P = 0.0f;
-  }
-
   p->primitives.gradients.rho[0] = 0.0f;
   p->primitives.gradients.rho[1] = 0.0f;
   p->primitives.gradients.rho[2] = 0.0f;

src/hydro/Gizmo/hydro_iact.h

@@ -135,8 +135,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
       Bi[k][l] = pi->geometry.matrix_E[k][l];
       Bj[k][l] = pj->geometry.matrix_E[k][l];
     }
-    vi[k] = pi->v[k]; /* particle velocities */
-    vj[k] = pj->v[k];
+    vi[k] = pi->force.v_full[k]; /* particle velocities */
+    vj[k] = pj->force.v_full[k];
   }
   Vi = pi->geometry.volume;
   Vj = pj->geometry.volume;

src/hydro/Gizmo/hydro_io.h

@@ -70,7 +70,7 @@ float convert_A(struct engine* e, struct part* p) {
 void hydro_write_particles(struct part* parts, struct io_props* list,
                            int* num_fields) {
 
-  *num_fields = 12;
+  *num_fields = 13;
 
   /* List what we want to write */
   list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -98,6 +98,8 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
       "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, primitives.P, convert_A);
   list[11] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
                                   parts, primitives.P);
+  list[12] = io_make_output_field("TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY,
+                                  parts, conserved.energy);
 }
 
 /**

src/hydro/Gizmo/hydro_part.h

@@ -159,6 +159,9 @@ struct part {
     /* Physical time step of the particle. */
     float dt;
 
+    /* Actual velocity of the particle. */
+    float v_full[3];
+
   } force;
 
   /* Particle mass (this field is also part of the conserved quantities...). */

src/hydro/Gizmo/hydro_slope_limiters.h

@@ -20,8 +20,8 @@
 #ifndef SWIFT_HYDRO_SLOPE_LIMITERS_H
 #define SWIFT_HYDRO_SLOPE_LIMITERS_H
 
-//#define PER_FACE_LIMITER
-//#define CELL_WIDE_LIMITER
+#define PER_FACE_LIMITER
+#define CELL_WIDE_LIMITER
 
 #ifdef PER_FACE_LIMITER
 #include "hydro_slope_limiters_face.h"