Created gizmo-mfm, a MFM copy of gizmo-mfv.

configure.ac

@@ -855,7 +855,7 @@ fi
 # Hydro scheme.
 AC_ARG_WITH([hydro],
    [AS_HELP_STRING([--with-hydro=<scheme>],
-      [Hydro dynamics to use @<:@gadget2, minimal, hopkins, default, gizmo-mfv, shadowfax, minimal-multi-mat, debug default: gadget2@:>@]
+      [Hydro dynamics to use @<:@gadget2, minimal, hopkins, default, gizmo-mfv, gizmo-mfm, shadowfax, minimal-multi-mat, debug default: gadget2@:>@]
    )],
    [with_hydro="$withval"],
    [with_hydro="gadget2"]
@@ -885,6 +885,9 @@ case "$with_hydro" in
    gizmo-mfv)
       AC_DEFINE([GIZMO_MFV_SPH], [1], [GIZMO MFV SPH])
    ;;
+   gizmo-mfm)
+      AC_DEFINE([GIZMO_MFM_SPH], [1], [GIZMO MFM SPH])
+   ;;
    shadowfax)
       AC_DEFINE([SHADOWFAX_SPH], [1], [Shadowfax SPH])
    ;;

src/Makefile.am

@@ -92,6 +92,17 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
                  hydro/GizmoMFV/hydro_slope_limiters.h \
                  hydro/GizmoMFV/hydro_unphysical.h \
                  hydro/GizmoMFV/hydro_velocities.h \
+		 hydro/GizmoMFM/hydro.h hydro/GizmoMFM/hydro_iact.h \
+                 hydro/GizmoMFM/hydro_io.h hydro/GizmoMFM/hydro_debug.h \
+                 hydro/GizmoMFM/hydro_part.h \
+                 hydro/GizmoMFM/hydro_gradients_gizmo.h \
+                 hydro/GizmoMFM/hydro_gradients.h \
+                 hydro/GizmoMFM/hydro_gradients_sph.h \
+                 hydro/GizmoMFM/hydro_slope_limiters_cell.h \
+                 hydro/GizmoMFM/hydro_slope_limiters_face.h \
+                 hydro/GizmoMFM/hydro_slope_limiters.h \
+                 hydro/GizmoMFM/hydro_unphysical.h \
+                 hydro/GizmoMFM/hydro_velocities.h \
                  hydro/Shadowswift/hydro_debug.h \
                  hydro/Shadowswift/hydro_gradients.h hydro/Shadowswift/hydro.h \
                  hydro/Shadowswift/hydro_iact.h \

src/const.h

@@ -54,8 +54,6 @@
 //#define GIZMO_FIX_PARTICLES
 /* Try to keep cells regular by adding a correction velocity. */
 //#define GIZMO_STEER_MOTION
-/* Use a meshless finite mass method. */
-//#define GIZMO_MFM
 /* Use the total energy instead of the thermal energy as conserved variable. */
 //#define GIZMO_TOTAL_ENERGY
 

src/debug.c

@@ -50,6 +50,8 @@
 #include "./hydro/Default/hydro_debug.h"
 #elif defined(GIZMO_MFV_SPH)
 #include "./hydro/GizmoMFV/hydro_debug.h"
+#elif defined(GIZMO_MFM_SPH)
+#include "./hydro/GizmoMFM/hydro_debug.h"
 #elif defined(SHADOWFAX_SPH)
 #include "./hydro/Shadowswift/hydro_debug.h"
 #elif defined(MINIMAL_MULTI_MAT_SPH)

src/hydro.h

@@ -49,6 +49,10 @@
 #include "./hydro/GizmoMFV/hydro.h"
 #include "./hydro/GizmoMFV/hydro_iact.h"
 #define SPH_IMPLEMENTATION "GIZMO MFV (Hopkins 2015)"
+#elif defined(GIZMO_MFM_SPH)
+#include "./hydro/GizmoMFM/hydro.h"
+#include "./hydro/GizmoMFM/hydro_iact.h"
+#define SPH_IMPLEMENTATION "GIZMO MFM (Hopkins 2015)"
 #elif defined(SHADOWFAX_SPH)
 #include "./hydro/Shadowswift/hydro.h"
 #include "./hydro/Shadowswift/hydro_iact.h"

src/hydro/GizmoMFM/hydro_debug.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GIZMO_MFM_HYDRO_DEBUG_H
+#define SWIFT_GIZMO_MFM_HYDRO_DEBUG_H
+
+__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], "
+      "h=%.3e, "
+      "time_bin=%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={"
+      "wcount_dh=%.3e, "
+      "wcount=%.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->h, p->time_bin, 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.wcount_dh, p->density.wcount);
+}
+
+#endif /* SWIFT_GIZMO_MFM_HYDRO_DEBUG_H */

src/hydro/GizmoMFM/hydro_gradients.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_HYDRO_GIZMO_MFM_GRADIENTS_H
+#define SWIFT_HYDRO_GIZMO_MFM_GRADIENTS_H
+
+#include "hydro_slope_limiters.h"
+#include "hydro_unphysical.h"
+#include "riemann.h"
+
+#if defined(GRADIENTS_SPH)
+
+#define HYDRO_GRADIENT_IMPLEMENTATION "SPH gradients (Price 2012)"
+#include "hydro_gradients_sph.h"
+
+#elif defined(GRADIENTS_GIZMO)
+
+#define HYDRO_GRADIENT_IMPLEMENTATION "GIZMO gradients (Hopkins 2015)"
+#include "hydro_gradients_gizmo.h"
+
+#else
+
+/* No gradients. Perfectly acceptable, but we have to provide empty functions */
+#define HYDRO_GRADIENT_IMPLEMENTATION "No gradients (first order scheme)"
+
+/**
+ * @brief Initialize gradient variables
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_init(
+    struct part *p) {}
+
+/**
+ * @brief Gradient calculations done during the neighbour loop
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_collect(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj) {}
+
+/**
+ * @brief Gradient calculations done during the neighbour loop: non-symmetric
+ * version
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
+                               struct part *restrict pi,
+                               struct part *restrict pj) {}
+
+/**
+ * @brief Finalize the gradient variables after all data have been collected
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
+    struct part *p) {}
+
+#endif
+
+/**
+ * @brief Gradients reconstruction. Is the same for all gradient types (although
+ * gradients_none does nothing, since all gradients are zero -- are they?).
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_predict(
+    struct part* restrict pi, struct part* restrict pj, float hi, float hj,
+    const float* dx, float r, const float* xij_i, float* Wi, float* Wj) {
+
+  /* perform gradient reconstruction in space and time */
+  /* Compute interface position (relative to pj, since we don't need the actual
+   * position) eqn. (8) */
+  const float xfac = hj / (hi + hj);
+  const float xij_j[3] = {xfac * dx[0], xfac * dx[1], xfac * dx[2]};
+
+  float dWi[5];
+  dWi[0] = pi->primitives.gradients.rho[0] * xij_i[0] +
+           pi->primitives.gradients.rho[1] * xij_i[1] +
+           pi->primitives.gradients.rho[2] * xij_i[2];
+  dWi[1] = pi->primitives.gradients.v[0][0] * xij_i[0] +
+           pi->primitives.gradients.v[0][1] * xij_i[1] +
+           pi->primitives.gradients.v[0][2] * xij_i[2];
+  dWi[2] = pi->primitives.gradients.v[1][0] * xij_i[0] +
+           pi->primitives.gradients.v[1][1] * xij_i[1] +
+           pi->primitives.gradients.v[1][2] * xij_i[2];
+  dWi[3] = pi->primitives.gradients.v[2][0] * xij_i[0] +
+           pi->primitives.gradients.v[2][1] * xij_i[1] +
+           pi->primitives.gradients.v[2][2] * xij_i[2];
+  dWi[4] = pi->primitives.gradients.P[0] * xij_i[0] +
+           pi->primitives.gradients.P[1] * xij_i[1] +
+           pi->primitives.gradients.P[2] * xij_i[2];
+
+  float dWj[5];
+  dWj[0] = pj->primitives.gradients.rho[0] * xij_j[0] +
+           pj->primitives.gradients.rho[1] * xij_j[1] +
+           pj->primitives.gradients.rho[2] * xij_j[2];
+  dWj[1] = pj->primitives.gradients.v[0][0] * xij_j[0] +
+           pj->primitives.gradients.v[0][1] * xij_j[1] +
+           pj->primitives.gradients.v[0][2] * xij_j[2];
+  dWj[2] = pj->primitives.gradients.v[1][0] * xij_j[0] +
+           pj->primitives.gradients.v[1][1] * xij_j[1] +
+           pj->primitives.gradients.v[1][2] * xij_j[2];
+  dWj[3] = pj->primitives.gradients.v[2][0] * xij_j[0] +
+           pj->primitives.gradients.v[2][1] * xij_j[1] +
+           pj->primitives.gradients.v[2][2] * xij_j[2];
+  dWj[4] = pj->primitives.gradients.P[0] * xij_j[0] +
+           pj->primitives.gradients.P[1] * xij_j[1] +
+           pj->primitives.gradients.P[2] * xij_j[2];
+
+  /* Apply the slope limiter at this interface */
+  hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);
+
+  Wi[0] += dWi[0];
+  Wi[1] += dWi[1];
+  Wi[2] += dWi[2];
+  Wi[3] += dWi[3];
+  Wi[4] += dWi[4];
+
+  Wj[0] += dWj[0];
+  Wj[1] += dWj[1];
+  Wj[2] += dWj[2];
+  Wj[3] += dWj[3];
+  Wj[4] += dWj[4];
+
+  gizmo_check_physical_quantities("density", "pressure", Wi[0], Wi[1], Wi[2],
+                                  Wi[3], Wi[4]);
+  gizmo_check_physical_quantities("density", "pressure", Wj[0], Wj[1], Wj[2],
+                                  Wj[3], Wj[4]);
+}
+
+#endif /* SWIFT_HYDRO_GIZMO_MFM_GRADIENTS_H */

src/hydro/GizmoMFM/hydro_gradients_gizmo.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/**
+ * @brief Initialize gradient variables
+ *
+ * @param p Particle.
+ */
+#ifndef SWIFT_GIZMO_MFM_HYDRO_GRADIENTS_H
+#define SWIFT_GIZMO_MFM_HYDRO_GRADIENTS_H
+
+__attribute__((always_inline)) INLINE static void hydro_gradients_init(
+    struct part *p) {
+
+  p->primitives.gradients.rho[0] = 0.0f;
+  p->primitives.gradients.rho[1] = 0.0f;
+  p->primitives.gradients.rho[2] = 0.0f;
+
+  p->primitives.gradients.v[0][0] = 0.0f;
+  p->primitives.gradients.v[0][1] = 0.0f;
+  p->primitives.gradients.v[0][2] = 0.0f;
+
+  p->primitives.gradients.v[1][0] = 0.0f;
+  p->primitives.gradients.v[1][1] = 0.0f;
+  p->primitives.gradients.v[1][2] = 0.0f;
+
+  p->primitives.gradients.v[2][0] = 0.0f;
+  p->primitives.gradients.v[2][1] = 0.0f;
+  p->primitives.gradients.v[2][2] = 0.0f;
+
+  p->primitives.gradients.P[0] = 0.0f;
+  p->primitives.gradients.P[1] = 0.0f;
+  p->primitives.gradients.P[2] = 0.0f;
+
+  hydro_slope_limit_cell_init(p);
+}
+
+/**
+ * @brief Gradient calculations done during the neighbour loop
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_collect(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj) {
+
+  const float r_inv = 1.f / sqrtf(r2);
+  const float r = r2 * r_inv;
+
+  float wi, wj, wi_dx, wj_dx;
+  float Bi[3][3];
+  float Bj[3][3];
+  float Wi[5], Wj[5];
+
+  /* Initialize local variables */
+  for (int k = 0; k < 3; k++) {
+    for (int l = 0; l < 3; l++) {
+      Bi[k][l] = pi->geometry.matrix_E[k][l];
+      Bj[k][l] = pj->geometry.matrix_E[k][l];
+    }
+  }
+  Wi[0] = pi->primitives.rho;
+  Wi[1] = pi->primitives.v[0];
+  Wi[2] = pi->primitives.v[1];
+  Wi[3] = pi->primitives.v[2];
+  Wi[4] = pi->primitives.P;
+  Wj[0] = pj->primitives.rho;
+  Wj[1] = pj->primitives.v[0];
+  Wj[2] = pj->primitives.v[1];
+  Wj[3] = pj->primitives.v[2];
+  Wj[4] = pj->primitives.P;
+
+  /* Compute kernel of pi. */
+  const float hi_inv = 1.f / hi;
+  const float xi = r * hi_inv;
+  kernel_deval(xi, &wi, &wi_dx);
+
+  if (pi->density.wcorr > const_gizmo_min_wcorr) {
+    /* Compute gradients for pi */
+    /* there is a sign difference w.r.t. eqn. (6) because of the inverse
+     * definition of dx */
+    pi->primitives.gradients.rho[0] +=
+        (Wi[0] - Wj[0]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->primitives.gradients.rho[1] +=
+        (Wi[0] - Wj[0]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->primitives.gradients.rho[2] +=
+        (Wi[0] - Wj[0]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+    pi->primitives.gradients.v[0][0] +=
+        (Wi[1] - Wj[1]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->primitives.gradients.v[0][1] +=
+        (Wi[1] - Wj[1]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->primitives.gradients.v[0][2] +=
+        (Wi[1] - Wj[1]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+    pi->primitives.gradients.v[1][0] +=
+        (Wi[2] - Wj[2]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->primitives.gradients.v[1][1] +=
+        (Wi[2] - Wj[2]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->primitives.gradients.v[1][2] +=
+        (Wi[2] - Wj[2]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+    pi->primitives.gradients.v[2][0] +=
+        (Wi[3] - Wj[3]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->primitives.gradients.v[2][1] +=
+        (Wi[3] - Wj[3]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->primitives.gradients.v[2][2] +=
+        (Wi[3] - Wj[3]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+    pi->primitives.gradients.P[0] +=
+        (Wi[4] - Wj[4]) * wi *
+        (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+    pi->primitives.gradients.P[1] +=
+        (Wi[4] - Wj[4]) * wi *
+        (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+    pi->primitives.gradients.P[2] +=
+        (Wi[4] - Wj[4]) * wi *
+        (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+  } else {
+    /* The gradient matrix was not well-behaved, switch to SPH gradients */
+
+    pi->primitives.gradients.rho[0] -=
+        wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+    pi->primitives.gradients.rho[1] -=
+        wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+    pi->primitives.gradients.rho[2] -=
+        wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+
+    pi->primitives.gradients.v[0][0] -=
+        wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+    pi->primitives.gradients.v[0][1] -=
+        wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+    pi->primitives.gradients.v[0][2] -=
+        wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+
+    pi->primitives.gradients.v[1][0] -=
+        wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+    pi->primitives.gradients.v[1][1] -=
+        wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+    pi->primitives.gradients.v[1][2] -=
+        wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+
+    pi->primitives.gradients.v[2][0] -=
+        wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+    pi->primitives.gradients.v[2][1] -=
+        wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+    pi->primitives.gradients.v[2][2] -=
+        wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+
+    pi->primitives.gradients.P[0] -=
+        wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
+    pi->primitives.gradients.P[1] -=
+        wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
+    pi->primitives.gradients.P[2] -=
+        wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
+  }
+
+  hydro_slope_limit_cell_collect(pi, pj, r);
+
+  /* Compute kernel of pj. */
+  const float hj_inv = 1.f / hj;
+  const float xj = r * hj_inv;
+  kernel_deval(xj, &wj, &wj_dx);
+
+  if (pj->density.wcorr > const_gizmo_min_wcorr) {
+    /* Compute gradients for pj */
+    /* there is no sign difference w.r.t. eqn. (6) because dx is now what we
+     * want
+     * it to be */
+    pj->primitives.gradients.rho[0] +=
+        (Wi[0] - Wj[0]) * wj *
+        (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+    pj->primitives.gradients.rho[1] +=
+        (Wi[0] - Wj[0]) * wj *
+        (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+    pj->primitives.gradients.rho[2] +=
+        (Wi[0] - Wj[0]) * wj *
+        (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+
+    pj->primitives.gradients.v[0][0] +=
+        (Wi[1] - Wj[1]) * wj *
+        (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+    pj->primitives.gradients.v[0][1] +=
+        (Wi[1] - Wj[1]) * wj *
+        (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+    pj->primitives.gradients.v[0][2] +=
+        (Wi[1] - Wj[1]) * wj *
+        (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+    pj->primitives.gradients.v[1][0] +=
+        (Wi[2] - Wj[2]) * wj *
+        (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+    pj->primitives.gradients.v[1][1] +=
+        (Wi[2] - Wj[2]) * wj *
+        (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+    pj->primitives.gradients.v[1][2] +=
+        (Wi[2] - Wj[2]) * wj *
+        (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+    pj->primitives.gradients.v[2][0] +=
+        (Wi[3] - Wj[3]) * wj *
+        (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+    pj->primitives.gradients.v[2][1] +=
+        (Wi[3] - Wj[3]) * wj *
+        (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+    pj->primitives.gradients.v[2][2] +=
+        (Wi[3] - Wj[3]) * wj *
+        (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+
+    pj->primitives.gradients.P[0] +=
+        (Wi[4] - Wj[4]) * wj *
+        (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+    pj->primitives.gradients.P[1] +=
+        (Wi[4] - Wj[4]) * wj *
+        (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+    pj->primitives.gradients.P[2] +=
+        (Wi[4] - Wj[4]) * wj *
+        (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+
+  } else {
+    /* SPH gradients */
+
+    pj->primitives.gradients.rho[0] -=
+        wj_dx * dx[0] * (pi->primitives