diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 01ddb6882f6e6a1ab607be314c2bf8c7a9cc9686..1fbec63fcae2d851c9acca04fb32f3e6b056d0ec 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -116,7 +116,6 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   /* Some smoothing length multiples. */
   const float h = p->h;
   const float ih = 1.0f / h;
-  const float ih2 = ih * ih;
   const float ihdim = pow_dimension(ih);
 
   float volume;
@@ -144,7 +143,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
 
   invert_dimension_by_dimension_matrix(p->geometry.matrix_E);
 
-  hydro_gradients_prepare_force_loop(p, ih2, volume);
+  hydro_gradients_prepare_force_loop(p, ih, volume);
 
   /* compute primitive variables */
   /* eqns (3)-(5) */
@@ -390,14 +389,18 @@ __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->a_hydro[1] = p->conserved.momentum[1] / p->conserved.mass - p->v[1];
-  p->a_hydro[2] = p->conserved.momentum[2] / p->conserved.mass - p->v[2];
+  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;
 }
 
 __attribute__((always_inline)) INLINE static void hydro_kick_extra(
diff --git a/src/hydro/Gizmo/hydro_gradients.h b/src/hydro/Gizmo/hydro_gradients.h
index f55bde5134307c4e07bf22a3f79aeb06d8693b77..11f79ebff2f2f0e9b0d6c342e8ab20bb67ab98eb 100644
--- a/src/hydro/Gizmo/hydro_gradients.h
+++ b/src/hydro/Gizmo/hydro_gradients.h
@@ -20,7 +20,7 @@
 #ifndef SWIFT_HYDRO_GRADIENTS_H
 #define SWIFT_HYDRO_GRADIENTS_H
 
-#define SPH_GRADIENTS
+//#define SPH_GRADIENTS
 
 #include "hydro_slope_limiters.h"
 
@@ -168,12 +168,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
                                             pj->primitives.gradients.v[1][1] +
                                             pj->primitives.gradients.v[2][2]));
 
-  //    printf("WL: %g %g %g %g %g\n", Wi[0], Wi[1], Wi[2], Wi[3], Wi[4]);
-  //    printf("WR: %g %g %g %g %g\n", Wj[0], Wj[1], Wj[2], Wj[3], Wj[4]);
-
-  //    printf("dWL: %g %g %g %g %g\n", dWi[0], dWi[1], dWi[2], dWi[3], dWi[4]);
-  //    printf("dWR: %g %g %g %g %g\n", dWj[0], dWj[1], dWj[2], dWj[3], dWj[4]);
-
   Wi[0] += dWi[0];
   Wi[1] += dWi[1];
   Wi[2] += dWi[2];
diff --git a/src/hydro/Gizmo/hydro_gradients_sph.h b/src/hydro/Gizmo/hydro_gradients_sph.h
index 29e36ec1903dbfa67b012b321c4a2d070dae91da..1671d9e43f9c0367495373b430ed11b71be2ef48 100644
--- a/src/hydro/Gizmo/hydro_gradients_sph.h
+++ b/src/hydro/Gizmo/hydro_gradients_sph.h
@@ -157,28 +157,30 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_density_loop(
  * @brief Calculations done before the force loop
  */
 __attribute__((always_inline)) INLINE static void
-hydro_gradients_prepare_force_loop(struct part *p, float ih2, float volume) {
+hydro_gradients_prepare_force_loop(struct part *p, float ih, float volume) {
+
+  const float ihdimp1 = pow_dimension_plus_one(ih);
 
   /* finalize gradients by multiplying with volume */
-  p->primitives.gradients.rho[0] *= ih2 * ih2 * volume;
-  p->primitives.gradients.rho[1] *= ih2 * ih2 * volume;
-  p->primitives.gradients.rho[2] *= ih2 * ih2 * volume;
+  p->primitives.gradients.rho[0] *= ihdimp1 * volume;
+  p->primitives.gradients.rho[1] *= ihdimp1 * volume;
+  p->primitives.gradients.rho[2] *= ihdimp1 * volume;
 
-  p->primitives.gradients.v[0][0] *= ih2 * ih2 * volume;
-  p->primitives.gradients.v[0][1] *= ih2 * ih2 * volume;
-  p->primitives.gradients.v[0][2] *= ih2 * ih2 * volume;
+  p->primitives.gradients.v[0][0] *= ihdimp1 * volume;
+  p->primitives.gradients.v[0][1] *= ihdimp1 * volume;
+  p->primitives.gradients.v[0][2] *= ihdimp1 * volume;
 
-  p->primitives.gradients.v[1][0] *= ih2 * ih2 * volume;
-  p->primitives.gradients.v[1][1] *= ih2 * ih2 * volume;
-  p->primitives.gradients.v[1][2] *= ih2 * ih2 * volume;
+  p->primitives.gradients.v[1][0] *= ihdimp1 * volume;
+  p->primitives.gradients.v[1][1] *= ihdimp1 * volume;
+  p->primitives.gradients.v[1][2] *= ihdimp1 * volume;
 
-  p->primitives.gradients.v[2][0] *= ih2 * ih2 * volume;
-  p->primitives.gradients.v[2][1] *= ih2 * ih2 * volume;
-  p->primitives.gradients.v[2][2] *= ih2 * ih2 * volume;
+  p->primitives.gradients.v[2][0] *= ihdimp1 * volume;
+  p->primitives.gradients.v[2][1] *= ihdimp1 * volume;
+  p->primitives.gradients.v[2][2] *= ihdimp1 * volume;
 
-  p->primitives.gradients.P[0] *= ih2 * ih2 * volume;
-  p->primitives.gradients.P[1] *= ih2 * ih2 * volume;
-  p->primitives.gradients.P[2] *= ih2 * ih2 * volume;
+  p->primitives.gradients.P[0] *= ihdimp1 * volume;
+  p->primitives.gradients.P[1] *= ihdimp1 * volume;
+  p->primitives.gradients.P[2] *= ihdimp1 * volume;
 
   hydro_slope_limit_cell(p);
 }
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index 82b26f3c78a9835f7e1359921e641e978550c8d0..d039ccf2077b3336491f055f58403015e2472577 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -113,8 +113,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
 
   float r = sqrtf(r2);
   float xi, xj;
-  float hi_inv, hi_inv2;
-  float hj_inv, hj_inv2;
+  float hi_inv, hi_inv_dim;
+  float hj_inv, hj_inv_dim;
   float wi, wj, wi_dx, wj_dx;
   int k, l;
   float A[3];
@@ -177,13 +177,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
 
   /* Compute kernel of pi. */
   hi_inv = 1.0 / hi;
-  hi_inv2 = hi_inv * hi_inv;
+  hi_inv_dim = pow_dimension(hi_inv);
   xi = r * hi_inv;
   kernel_deval(xi, &wi, &wi_dx);
 
   /* Compute kernel of pj. */
   hj_inv = 1.0 / hj;
-  hj_inv2 = hj_inv * hj_inv;
+  hj_inv_dim = pow_dimension(hj_inv);
   xj = r * hj_inv;
   kernel_deval(xj, &wj, &wj_dx);
 
@@ -193,9 +193,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   for (k = 0; k < 3; k++) {
     /* we add a minus sign since dx is pi->x - pj->x */
     A[k] = -Vi * (Bi[k][0] * dx[0] + Bi[k][1] * dx[1] + Bi[k][2] * dx[2]) * wi *
-               hi_inv * hi_inv2 -
+               hi_inv_dim -
            Vj * (Bj[k][0] * dx[0] + Bj[k][1] * dx[1] + Bj[k][2] * dx[2]) * wj *
-               hj_inv * hj_inv2;
+               hj_inv_dim;
     Anorm += A[k] * A[k];
   }
 
@@ -209,13 +209,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   Anorm = sqrtf(Anorm);
   for (k = 0; k < 3; k++) n_unit[k] = A[k] / Anorm;
 
-#ifdef PRINT_ID
-  if (pi->id == PRINT_ID || pj->id == PRINT_ID) {
-    printf("pi: %g %g %g\npj: %g %g %g\nA = %g %g %g\n", pi->x[0], pi->x[1],
-           pi->x[2], pj->x[0], pj->x[1], pj->x[2], A[0], A[1], A[2]);
-  }
-#endif
-
   /* Compute interface position (relative to pi, since we don't need the actual
    * position) */
   /* eqn. (8) */
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index bba7b5e6b961773162bf065ba9254731d3bd1e93..dd7f46ceb0f9809f1eba013e85572e5057c9ccd0 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -56,6 +56,10 @@ float convert_u(struct engine* e, struct part* p) {
   return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
 }
 
+float convert_A(struct engine* e, struct part* p) {
+  return p->primitives.P / pow_gamma(p->primitives.rho);
+}
+
 /**
  * @brief Specifies which particle fields to write to a dataset
  *
@@ -66,7 +70,7 @@ float convert_u(struct engine* e, struct part* p) {
 void hydro_write_particles(struct part* parts, struct io_props* list,
                            int* num_fields) {
 
-  *num_fields = 10;
+  *num_fields = 12;
 
   /* List what we want to write */
   list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -90,6 +94,10 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
                                  geometry.volume);
   list[9] = io_make_output_field("GradDensity", FLOAT, 3, UNIT_CONV_DENSITY,
                                  parts, primitives.gradients.rho);
+  list[10] = io_make_output_field_convert_part(
+      "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);
 }
 
 /**
diff --git a/src/riemann.h b/src/riemann.h
index d647b021167317d14f4cd7316d09c247794f3d23..06acc2db88243bbda4ede24c3fb9cf70868b666c 100644
--- a/src/riemann.h
+++ b/src/riemann.h
@@ -27,7 +27,7 @@
 #include "stdio.h"
 #include "stdlib.h"
 
-#define HLLC_SOLVER
+#define EXACT_SOLVER
 
 #ifdef EXACT_SOLVER
 #include "riemann/riemann_exact.h"