diff --git a/src/hydro/Gizmo/hydro_debug.h b/src/hydro/Gizmo/hydro_debug.h
index e42ebb875f257001bf346c82a8c7559b47ef6b86..db672bab16662961a99558b563b9320d4e75ecc4 100644
--- a/src/hydro/Gizmo/hydro_debug.h
+++ b/src/hydro/Gizmo/hydro_debug.h
@@ -21,7 +21,63 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
     const struct part* p, const struct xpart* xp) {
   printf(
       "x=[%.16e,%.16e,%.16e], "
-      "v=[%.3e,%.3e,%.3e], a=[%.3e,%.3e,%.3e], volume=%.3e\n",
+      "v=[%.3e,%.3e,%.3e], "
+      "a=[%.3e,%.3e,%.3e], "
+      "h=%.3e, "
+      "ti_begin=%d, "
+      "ti_end=%d, "
+      "primitives={"
+      "v=[%.3e,%.3e,%.3e], "
+      "rho=%.3e, "
+      "P=%.3e, "
+      "gradients={"
+      "rho=[%.3e,%.3e,%.3e], "
+      "v=[[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e]], "
+      "P=[%.3e,%.3e,%.3e]}, "
+      "limiter={"
+      "rho=[%.3e,%.3e], "
+      "v=[[%.3e,%.3e],[%.3e,%.3e],[%.3e,%.3e]], "
+      "P=[%.3e,%.3e], "
+      "maxr=%.3e}}, "
+      "conserved={"
+      "momentum=[%.3e,%.3e,%.3e], "
+      "mass=%.3e, "
+      "energy=%.3e}, "
+      "geometry={"
+      "volume=%.3e, "
+      "matrix_E=[[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e]]}, "
+      "timestepvars={"
+      "vmax=%.3e}, "
+      "density={"
+      "div_v=%.3e, "
+      "wcount_dh=%.3e, "
+      "curl_v=[%.3e,%.3e,%.3e], "
+      "wcount=%.3e}, "
+      "mass=%.3e\n",
       p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->a_hydro[0],
-      p->a_hydro[1], p->a_hydro[2], p->geometry.volume);
+      p->a_hydro[1], p->a_hydro[2], p->h, p->ti_begin, p->ti_end,
+      p->primitives.v[0], p->primitives.v[1], p->primitives.v[2],
+      p->primitives.rho, p->primitives.P, p->primitives.gradients.rho[0],
+      p->primitives.gradients.rho[1], p->primitives.gradients.rho[2],
+      p->primitives.gradients.v[0][0], p->primitives.gradients.v[0][1],
+      p->primitives.gradients.v[0][2], p->primitives.gradients.v[1][0],
+      p->primitives.gradients.v[1][1], p->primitives.gradients.v[1][2],
+      p->primitives.gradients.v[2][0], p->primitives.gradients.v[2][1],
+      p->primitives.gradients.v[2][2], p->primitives.gradients.P[0],
+      p->primitives.gradients.P[1], p->primitives.gradients.P[2],
+      p->primitives.limiter.rho[0], p->primitives.limiter.rho[1],
+      p->primitives.limiter.v[0][0], p->primitives.limiter.v[0][1],
+      p->primitives.limiter.v[1][0], p->primitives.limiter.v[1][1],
+      p->primitives.limiter.v[2][0], p->primitives.limiter.v[2][1],
+      p->primitives.limiter.P[0], p->primitives.limiter.P[1],
+      p->primitives.limiter.maxr, p->conserved.momentum[0],
+      p->conserved.momentum[1], p->conserved.momentum[2], p->conserved.mass,
+      p->conserved.energy, p->geometry.volume, p->geometry.matrix_E[0][0],
+      p->geometry.matrix_E[0][1], p->geometry.matrix_E[0][2],
+      p->geometry.matrix_E[1][0], p->geometry.matrix_E[1][1],
+      p->geometry.matrix_E[1][2], p->geometry.matrix_E[2][0],
+      p->geometry.matrix_E[2][1], p->geometry.matrix_E[2][2],
+      p->timestepvars.vmax, p->density.div_v, p->density.wcount_dh,
+      p->density.curl_v[0], p->density.curl_v[1], p->density.curl_v[2],
+      p->density.wcount, p->mass);
 }
diff --git a/tests/testSymmetry.c b/tests/testSymmetry.c
index eb3fab6becca08e9ef87e7c60cc8c04bd2a0290c..f03bc68eac9d05e1e2c3b07b221bd0b5172e6838 100644
--- a/tests/testSymmetry.c
+++ b/tests/testSymmetry.c
@@ -46,6 +46,59 @@ int main(int argc, char *argv[]) {
   pi.id = 1;
   pj.id = 2;
 
+#if defined(GIZMO_SPH)
+  /* Give the primitive variables sensible values, since the Riemann solver does
+     not like negative densities and pressures */
+  pi.primitives.rho = random_uniform(0.1f, 1.0f);
+  pi.primitives.v[0] = random_uniform(-10.0f, 10.0f);
+  pi.primitives.v[1] = random_uniform(-10.0f, 10.0f);
+  pi.primitives.v[2] = random_uniform(-10.0f, 10.0f);
+  pi.primitives.P = random_uniform(0.1f, 1.0f);
+  /*  pj.primitives.rho = random_uniform(0.1f, 1.0f);
+    pj.primitives.v[0] = random_uniform(-10.0f, 10.0f);
+    pj.primitives.v[1] = random_uniform(-10.0f, 10.0f);
+    pj.primitives.v[2] = random_uniform(-10.0f, 10.0f);
+    pj.primitives.P = random_uniform(0.1f, 1.0f);*/
+  /* make the values for pj the same, since otherwise we suffer from the swap of
+     left and right states in the Riemann solver */
+  pj.primitives.rho = pi.primitives.rho;
+  pj.primitives.v[0] = pi.primitives.v[0];
+  pj.primitives.v[1] = pi.primitives.v[1];
+  pj.primitives.v[2] = pi.primitives.v[2];
+  pj.primitives.P = pi.primitives.P;
+  /* make gradients zero */
+  pi.primitives.gradients.rho[0] = 0.0f;
+  pi.primitives.gradients.rho[1] = 0.0f;
+  pi.primitives.gradients.rho[2] = 0.0f;
+  pi.primitives.gradients.v[0][0] = 0.0f;
+  pi.primitives.gradients.v[0][1] = 0.0f;
+  pi.primitives.gradients.v[0][2] = 0.0f;
+  pi.primitives.gradients.v[1][0] = 0.0f;
+  pi.primitives.gradients.v[1][1] = 0.0f;
+  pi.primitives.gradients.v[1][2] = 0.0f;
+  pi.primitives.gradients.v[2][0] = 0.0f;
+  pi.primitives.gradients.v[2][1] = 0.0f;
+  pi.primitives.gradients.v[2][2] = 0.0f;
+  pi.primitives.gradients.P[0] = 0.0f;
+  pi.primitives.gradients.P[1] = 0.0f;
+  pi.primitives.gradients.P[2] = 0.0f;
+  pj.primitives.gradients.rho[0] = 0.0f;
+  pj.primitives.gradients.rho[1] = 0.0f;
+  pj.primitives.gradients.rho[2] = 0.0f;
+  pj.primitives.gradients.v[0][0] = 0.0f;
+  pj.primitives.gradients.v[0][1] = 0.0f;
+  pj.primitives.gradients.v[0][2] = 0.0f;
+  pj.primitives.gradients.v[1][0] = 0.0f;
+  pj.primitives.gradients.v[1][1] = 0.0f;
+  pj.primitives.gradients.v[1][2] = 0.0f;
+  pj.primitives.gradients.v[2][0] = 0.0f;
+  pj.primitives.gradients.v[2][1] = 0.0f;
+  pj.primitives.gradients.v[2][2] = 0.0f;
+  pj.primitives.gradients.P[0] = 0.0f;
+  pj.primitives.gradients.P[1] = 0.0f;
+  pj.primitives.gradients.P[2] = 0.0f;
+#endif
+
   /* Make an xpart companion */
   struct xpart xpi, xpj;
   bzero(&xpi, sizeof(struct xpart));
@@ -104,8 +157,16 @@ int main(int argc, char *argv[]) {
   i_ok = memcmp(&pi, &pi2, sizeof(struct part));
   j_ok = memcmp(&pj, &pj2, sizeof(struct part));
 
-  if (i_ok) error("Particles 'pi' do not match after force");
-  if (j_ok) error("Particles 'pj' do not match after force");
+  if (i_ok) {
+    printParticle_single(&pi, &xpi);
+    printParticle_single(&pi2, &xpi);
+    error("Particles 'pi' do not match after force");
+  }
+  if (j_ok) {
+    printParticle_single(&pj, &xpj);
+    printParticle_single(&pj2, &xpj);
+    error("Particles 'pj' do not match after force");
+  }
 
   return 0;
 }