diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 1fbec63fcae2d851c9acca04fb32f3e6b056d0ec..ca795ab00860d3559d755dc4c3c44453b29743f5 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/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(
diff --git a/src/hydro/Gizmo/hydro_gradients_sph.h b/src/hydro/Gizmo/hydro_gradients_sph.h
index 1671d9e43f9c0367495373b430ed11b71be2ef48..f8dc7ea00458dd11c45aa0a340e6139376a286ae 100644
--- a/src/hydro/Gizmo/hydro_gradients_sph.h
+++ b/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;
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index d039ccf2077b3336491f055f58403015e2472577..792d5c4ed068eb4439204f002939fb30202f3353 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/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;
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index dd7f46ceb0f9809f1eba013e85572e5057c9ccd0..c94451aad943e7603f3393a400fe5864e624d16b 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/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);
 }
 
 /**
diff --git a/src/hydro/Gizmo/hydro_part.h b/src/hydro/Gizmo/hydro_part.h
index 8980bd4a4913f5e17b6dd9770e018c95c5433fe1..ce1d7a5e00568c69acd1d53ac429dd35d1c56900 100644
--- a/src/hydro/Gizmo/hydro_part.h
+++ b/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...). */
diff --git a/src/hydro/Gizmo/hydro_slope_limiters.h b/src/hydro/Gizmo/hydro_slope_limiters.h
index b341d2f15ee9d57ca08ecfef6f738633cc2ad302..357ba2a84eb0a0c340906ddfa9d20e6a4153a405 100644
--- a/src/hydro/Gizmo/hydro_slope_limiters.h
+++ b/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"