diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index a2fd75e9c806a5666f56891ea9d8e12c809d71c8..7aa50d96b645879b6b3dc3e450b1b9ff64b81cd1 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -62,6 +62,9 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
 __attribute__((always_inline)) INLINE static void hydro_init_part(
     struct part* p) {
 
+  /* We need to do this before we reset the volume */
+  hydro_gradients_init_density_loop(p);
+
   p->density.wcount = 0.0f;
   p->density.wcount_dh = 0.0f;
   p->geometry.volume = 0.0f;
@@ -74,8 +77,6 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
   p->geometry.matrix_E[2][0] = 0.0f;
   p->geometry.matrix_E[2][1] = 0.0f;
   p->geometry.matrix_E[2][2] = 0.0f;
-
-  hydro_gradients_init_density_loop(p);
 }
 
 /**
diff --git a/src/hydro/Gizmo/hydro_gradients_sph.h b/src/hydro/Gizmo/hydro_gradients_sph.h
index 8c9d7be052920c9fc86216c63b55179c4dfd2153..29e36ec1903dbfa67b012b321c4a2d070dae91da 100644
--- a/src/hydro/Gizmo/hydro_gradients_sph.h
+++ b/src/hydro/Gizmo/hydro_gradients_sph.h
@@ -38,6 +38,12 @@ hydro_gradients_init_density_loop(struct part *p) {
                 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;
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index 0f847e1b8593f607cd8c4ddb48d418bda581252f..82b26f3c78a9835f7e1359921e641e978550c8d0 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -155,8 +155,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   dtj = pj->force.dt;
 
   /* calculate the maximal signal velocity */
-  vmax =
-      sqrtf(hydro_gamma * Wi[4] / Wi[0]) + sqrtf(hydro_gamma * Wj[4] / Wj[0]);
+  if (Wi[0] && Wj[0]) {
+    vmax =
+        sqrtf(hydro_gamma * Wi[4] / Wi[0]) + sqrtf(hydro_gamma * Wj[4] / Wj[0]);
+  } else {
+    vmax = 0.0f;
+  }
   dvdotdx = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
             (Wi[3] - Wj[3]) * dx[2];
   if (dvdotdx > 0.) {
diff --git a/tests/test125cells.c b/tests/test125cells.c
index f2cbe6606fc0467339fcf363b624154b90d463bf..a839c84f734cb6907c7807e6ed12253168c811e3 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -253,6 +253,25 @@ struct cell *make_cell(size_t n, const double offset[3], double size, double h,
         xpart->v_full[1] = part->v[1];
         xpart->v_full[2] = part->v[2];
         hydro_first_init_part(part, xpart);
+
+#if defined(GIZMO_SPH)
+        part->geometry.volume = volume;
+        part->primitives.rho = density;
+        part->primitives.v[0] = part->v[0];
+        part->primitives.v[1] = part->v[1];
+        part->primitives.v[2] = part->v[2];
+        part->conserved.mass = part->mass;
+        part->conserved.momentum[0] = part->conserved.mass * part->v[0];
+        part->conserved.momentum[1] = part->conserved.mass * part->v[1];
+        part->conserved.momentum[2] = part->conserved.mass * part->v[2];
+        part->conserved.energy =
+            part->primitives.P / hydro_gamma_minus_one * volume +
+            0.5f * (part->conserved.momentum[0] * part->conserved.momentum[0] +
+                    part->conserved.momentum[1] * part->conserved.momentum[1] +
+                    part->conserved.momentum[2] * part->conserved.momentum[2]) /
+                part->conserved.mass;
+#endif
+
         ++part;
         ++xpart;
       }