Tried to fix test125cells. Found a few bugs in the process.

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);
 }
 
 /**

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;

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.) {

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;
       }