Merge branch 'gizmo_mfm' into 'master'

GIZMO MFM

See merge request !549

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, 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"]
@@ -882,8 +882,11 @@ case "$with_hydro" in
    default)
       AC_DEFINE([DEFAULT_SPH], [1], [Default SPH])
    ;;
-   gizmo)
-      AC_DEFINE([GIZMO_SPH], [1], [GIZMO SPH])
+   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

@@ -81,17 +81,28 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
                  hydro/Gadget2/hydro_debug.h hydro/Gadget2/hydro_part.h \
 		 hydro/PressureEntropy/hydro.h hydro/PressureEntropy/hydro_iact.h hydro/PressureEntropy/hydro_io.h \
                  hydro/PressureEntropy/hydro_debug.h hydro/PressureEntropy/hydro_part.h \
-		 hydro/Gizmo/hydro.h hydro/Gizmo/hydro_iact.h \
-                 hydro/Gizmo/hydro_io.h hydro/Gizmo/hydro_debug.h \
-                 hydro/Gizmo/hydro_part.h \
-                 hydro/Gizmo/hydro_gradients_gizmo.h \
-                 hydro/Gizmo/hydro_gradients.h \
-                 hydro/Gizmo/hydro_gradients_sph.h \
-                 hydro/Gizmo/hydro_slope_limiters_cell.h \
-                 hydro/Gizmo/hydro_slope_limiters_face.h \
-                 hydro/Gizmo/hydro_slope_limiters.h \
-                 hydro/Gizmo/hydro_unphysical.h \
-                 hydro/Gizmo/hydro_velocities.h \
+		 hydro/GizmoMFV/hydro.h hydro/GizmoMFV/hydro_iact.h \
+                 hydro/GizmoMFV/hydro_io.h hydro/GizmoMFV/hydro_debug.h \
+                 hydro/GizmoMFV/hydro_part.h \
+                 hydro/GizmoMFV/hydro_gradients_gizmo.h \
+                 hydro/GizmoMFV/hydro_gradients.h \
+                 hydro/GizmoMFV/hydro_gradients_sph.h \
+                 hydro/GizmoMFV/hydro_slope_limiters_cell.h \
+                 hydro/GizmoMFV/hydro_slope_limiters_face.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

@@ -53,7 +53,8 @@
 /* This option disables particle movement */
 //#define GIZMO_FIX_PARTICLES
 /* Try to keep cells regular by adding a correction velocity. */
-#define GIZMO_STEER_MOTION
+//#define GIZMO_STEER_MOTION
+/* Use the total energy instead of the thermal energy as conserved variable. */
 //#define GIZMO_TOTAL_ENERGY
 
 /* Options to control handling of unphysical values (GIZMO_SPH only). */

src/debug.c

@@ -48,8 +48,10 @@
 #include "./hydro/PressureEntropy/hydro_debug.h"
 #elif defined(DEFAULT_SPH)
 #include "./hydro/Default/hydro_debug.h"
-#elif defined(GIZMO_SPH)
-#include "./hydro/Gizmo/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

@@ -45,10 +45,14 @@
 #include "./hydro/Default/hydro.h"
 #include "./hydro/Default/hydro_iact.h"
 #define SPH_IMPLEMENTATION "Default version of SPH"
-#elif defined(GIZMO_SPH)
-#include "./hydro/Gizmo/hydro.h"
-#include "./hydro/Gizmo/hydro_iact.h"
-#define SPH_IMPLEMENTATION "GIZMO (Hopkins 2015)"
+#elif defined(GIZMO_MFV_SPH)
+#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/Gizmo/hydro.h → src/hydro/GizmoMFM/hydro.h

@@ -17,8 +17,8 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
-#ifndef SWIFT_GIZMO_HYDRO_H
-#define SWIFT_GIZMO_HYDRO_H
+#ifndef SWIFT_GIZMO_MFM_HYDRO_H
+#define SWIFT_GIZMO_MFM_HYDRO_H
 
 #include "adiabatic_index.h"
 #include "approx_math.h"
@@ -972,4 +972,4 @@ hydro_set_init_internal_energy(struct part* p, float u_init) {
   p->primitives.P = hydro_gamma_minus_one * p->primitives.rho * u_init;
 }
 
-#endif /* SWIFT_GIZMO_HYDRO_H */
+#endif /* SWIFT_GIZMO_MFM_HYDRO_H */

src/hydro/Gizmo/hydro_debug.h → src/hydro/GizmoMFM/hydro_debug.h

@@ -16,8 +16,8 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
-#ifndef SWIFT_GIZMO_HYDRO_DEBUG_H
-#define SWIFT_GIZMO_HYDRO_DEBUG_H
+#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) {
@@ -78,4 +78,4 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       p->timestepvars.vmax, p->density.wcount_dh, p->density.wcount);
 }
 
-#endif /* SWIFT_GIZMO_HYDRO_DEBUG_H */
+#endif /* SWIFT_GIZMO_MFM_HYDRO_DEBUG_H */

src/hydro/Gizmo/hydro_gradients.h → src/hydro/GizmoMFM/hydro_gradients.h

@@ -17,8 +17,8 @@
  *
  ******************************************************************************/
 
-#ifndef SWIFT_HYDRO_GIZMO_GRADIENTS_H
-#define SWIFT_HYDRO_GIZMO_GRADIENTS_H
+#ifndef SWIFT_HYDRO_GIZMO_MFM_GRADIENTS_H
+#define SWIFT_HYDRO_GIZMO_MFM_GRADIENTS_H
 
 #include "hydro_slope_limiters.h"
 #include "hydro_unphysical.h"
@@ -156,4 +156,4 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
                                   Wj[3], Wj[4]);
 }
 
-#endif /* SWIFT_HYDRO_GIZMO_GRADIENTS_H */
+#endif /* SWIFT_HYDRO_GIZMO_MFM_GRADIENTS_H */

src/hydro/Gizmo/hydro_gradients_gizmo.h → src/hydro/GizmoMFM/hydro_gradients_gizmo.h

@@ -22,8 +22,8 @@
  *
  * @param p Particle.
  */
-#ifndef SWIFT_GIZMO_HYDRO_GRADIENTS_H
-#define SWIFT_GIZMO_HYDRO_GRADIENTS_H
+#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) {
@@ -485,4 +485,4 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
   hydro_slope_limit_cell(p);
 }
 
-#endif /* SWIFT_GIZMO_HYDRO_GRADIENTS_H */
+#endif /* SWIFT_GIZMO_MFM_HYDRO_GRADIENTS_H */

src/hydro/Gizmo/hydro_gradients_sph.h → src/hydro/GizmoMFM/hydro_gradients_sph.h

@@ -22,8 +22,8 @@
  *
  * @param p Particle.
  */
-#ifndef SWIFT_GIZMO_HYDRO_SPH_GRADIENTS_H
-#define SWIFT_GIZMO_HYDRO_SPH_GRADIENTS_H
+#ifndef SWIFT_GIZMO_MFM_HYDRO_SPH_GRADIENTS_H
+#define SWIFT_GIZMO_MFM_HYDRO_SPH_GRADIENTS_H
 
 __attribute__((always_inline)) INLINE static void hydro_gradients_init(
     struct part *p) {
@@ -253,4 +253,4 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
   hydro_slope_limit_cell(p);
 }
 
-#endif /* SWIFT_GIZMO_HYDRO_SPH_GRADIENTS_H */
+#endif /* SWIFT_GIZMO_MFM_HYDRO_SPH_GRADIENTS_H */

src/hydro/GizmoMFM/hydro_iact.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *                    Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ *
+ * 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_IACT_H
+#define SWIFT_GIZMO_MFM_HYDRO_IACT_H
+
+#include "adiabatic_index.h"
+#include "hydro_gradients.h"
+#include "riemann.h"
+
+#define GIZMO_VOLUME_CORRECTION
+
+/**
+ * @brief Calculate the volume interaction between particle i and particle j
+ *
+ * The volume is in essence the same as the weighted number of neighbours in a
+ * classical SPH density calculation.
+ *
+ * We also calculate the components of the matrix E, which is used for second
+ * order accurate gradient calculations and for the calculation of the interface
+ * surface areas.
+ *
+ * @param r2 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_density(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H) {
+
+  float wi, wj, wi_dx, wj_dx;
+
+  /* Get r and h inverse. */
+  const float r = sqrtf(r2);
+
+  /* Compute density of pi. */
+  const float hi_inv = 1.f / hi;
+  const float xi = r * hi_inv;
+  kernel_deval(xi, &wi, &wi_dx);
+
+  pi->density.wcount += wi;
+  pi->density.wcount_dh -= (hydro_dimension * wi + xi * wi_dx);
+
+  /* these are eqns. (1) and (2) in the summary */
+  pi->geometry.volume += wi;
+  for (int k = 0; k < 3; k++)
+    for (int l = 0; l < 3; l++)
+      pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
+
+  pi->geometry.centroid[0] -= dx[0] * wi;
+  pi->geometry.centroid[1] -= dx[1] * wi;
+  pi->geometry.centroid[2] -= dx[2] * wi;
+
+  /* Compute density of pj. */
+  const float hj_inv = 1.f / hj;
+  const float xj = r * hj_inv;
+  kernel_deval(xj, &wj, &wj_dx);
+
+  pj->density.wcount += wj;
+  pj->density.wcount_dh -= (hydro_dimension * wj + xj * wj_dx);
+
+  /* these are eqns. (1) and (2) in the summary */
+  pj->geometry.volume += wj;
+  for (int k = 0; k < 3; k++)
+    for (int l = 0; l < 3; l++)
+      pj->geometry.matrix_E[k][l] += dx[k] * dx[l] * wj;
+
+  pj->geometry.centroid[0] += dx[0] * wj;
+  pj->geometry.centroid[1] += dx[1] * wj;
+  pj->geometry.centroid[2] += dx[2] * wj;
+}
+
+/**
+ * @brief Calculate the volume interaction between particle i and particle j:
+ * non-symmetric version
+ *
+ * The volume is in essence the same as the weighted number of neighbours in a
+ * classical SPH density calculation.
+ *
+ * We also calculate the components of the matrix E, which is used for second
+ * order accurate gradient calculations and for the calculation of the interface
+ * surface areas.
+ *
+ * @param r2 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    const struct part *restrict pj, float a, float H) {
+
+  float wi, wi_dx;
+
+  /* Get r and h inverse. */
+  const float r = sqrtf(r2);
+
+  const float hi_inv = 1.f / hi;
+  const float xi = r * hi_inv;
+  kernel_deval(xi, &wi, &wi_dx);
+
+  pi->density.wcount += wi;
+  pi->density.wcount_dh -= (hydro_dimension * wi + xi * wi_dx);
+
+  /* these are eqns. (1) and (2) in the summary */
+  pi->geometry.volume += wi;
+  for (int k = 0; k < 3; k++)
+    for (int l = 0; l < 3; l++)
+      pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
+
+  pi->geometry.centroid[0] -= dx[0] * wi;
+  pi->geometry.centroid[1] -= dx[1] * wi;
+  pi->geometry.centroid[2] -= dx[2] * wi;
+}
+
+/**
+ * @brief Calculate the gradient interaction between particle i and particle j
+ *
+ * This method wraps around hydro_gradients_collect, which can be an empty
+ * method, in which case no gradients are used.
+ *
+ * @param r2 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_gradient(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H) {
+
+  hydro_gradients_collect(r2, dx, hi, hj, pi, pj);
+}
+
+/**
+ * @brief Calculate the gradient interaction between particle i and particle j:
+ * non-symmetric version
+ *
+ * This method wraps around hydro_gradients_nonsym_collect, which can be an
+ * empty method, in which case no gradients are used.
+ *
+ * @param r2 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H) {
+
+  hydro_gradients_nonsym_collect(r2, dx, hi, hj, pi, pj);
+}
+
+/**
+ * @brief Common part of the flux calculation between particle i and j
+ *
+ * Since the only difference between the symmetric and non-symmetric version
+ * of the flux calculation  is in the update of the conserved variables at the
+ * very end (which is not done for particle j if mode is 0), both
+ * runner_iact_force and runner_iact_nonsym_force call this method, with an
+ * appropriate mode.
+ *
+ * This method calculates the surface area of the interface between particle i
+ * and particle j, as well as the interface position and velocity. These are
+ * then used to reconstruct and predict the primitive variables, which are then
+ * fed to a Riemann solver that calculates a flux. This flux is used to update
+ * the conserved variables of particle i or both particles.
+ *
+ * This method also calculates the maximal velocity used to calculate the time
+ * step.
+ *
+ * @param r2 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, int mode, float a, float H) {
+
+  const float r_inv = 1.f / sqrtf(r2);
+  const float r = r2 * r_inv;
+
+  /* Initialize local variables */
+  float Bi[3][3];
+  float Bj[3][3];
+  float vi[3], vj[3];
+  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];
+    }
+    vi[k] = pi->v[k]; /* particle velocities */
+    vj[k] = pj->v[k];
+  }
+  const float Vi = pi->geometry.volume;
+  const float Vj = pj->geometry.volume;
+  float Wi[5], Wj[5];
+  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;
+
+  /* calculate the maximal signal velocity */
+  float vmax;
+  if (Wi[0] > 0.0f && Wj[0] > 0.0f) {
+    const float ci = gas_soundspeed_from_pressure(Wi[0], Wi[4]);
+    const float cj = gas_soundspeed_from_pressure(Wj[0], Wj[4]);
+    vmax = ci + cj;
+  } else
+    vmax = 0.0f;
+
+  float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] +
+               (pi->v[2] - pj->v[2]) * dx[2];
+
+  /* Velocity on the axis linking the particles */
+  float dvdotdx = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
+                  (Wi[3] - Wj[3]) * dx[2];
+  dvdotdx = min(dvdotdx, dvdr);
+
+  /* We only care about this velocity for particles moving towards each others
+   */
+  dvdotdx = min(dvdotdx, 0.f);
+
+  /* Get the signal velocity */
+  /* the magical factor 3 also appears in Gadget2 */
+  vmax -= 3.f * dvdotdx * r_inv;
+
+  /* Store the signal velocity */
+  pi->timestepvars.vmax = max(pi->timestepvars.vmax, vmax);
+  if (mode == 1) pj->timestepvars.vmax = max(pj->timestepvars.vmax, vmax);
+
+  /* Compute kernel of pi. */
+  float wi, wi_dx;
+  const float hi_inv = 1.f / hi;
+  const float hi_inv_dim = pow_dimension(hi_inv);
+  const float xi = r * hi_inv;
+  kernel_deval(xi, &wi, &wi_dx);
+
+  /* Compute kernel of pj. */
+  float wj, wj_dx;
+  const float hj_inv = 1.f / hj;
+  const float hj_inv_dim = pow_dimension(hj_inv);
+  const float xj = r * hj_inv;
+  kernel_deval(xj, &wj, &wj_dx);
+
+  /* Compute h_dt. We are going to use an SPH-like estimate of div_v for that */
+  const float hidp1 = pow_dimension_plus_one(hi_inv);
+  const float hjdp1 = pow_dimension_plus_one(hj_inv);
+  const float wi_dr = hidp1 * wi_dx;
+  const float wj_dr = hjdp1 * wj_dx;
+  dvdr *= r_inv;
+  if (pj->primitives.rho > 0.)
+    pi->force.h_dt -= pj->conserved.mass * dvdr / pj->primitives.rho * wi_dr;
+  if (mode == 1 && pi->primitives.rho > 0.)
+    pj->force.h_dt -= pi->conserved.mass * dvdr / pi->primitives.rho * wj_dr;
+
+  /* Compute (square of) area */
+  /* eqn. (7) */
+  float Anorm2 = 0.0f;
+  float A[3];
+  if (pi->density.wcorr > const_gizmo_min_wcorr &&
+      pj->density.wcorr > const_gizmo_min_wcorr) {
+    /* in principle, we use Vi and Vj as weights for the left and right
+       contributions to the generalized surface vector.
+       However, if Vi and Vj are very different (because they have very
+       different
+       smoothing lengths), then the expressions below are more stable. */
+    float Xi = Vi;
+    float Xj = Vj;
+#ifdef GIZMO_VOLUME_CORRECTION
+    if (fabsf(Vi - Vj) / min(Vi, Vj) > 1.5f * hydro_dimension) {
+      Xi = (Vi * hj + Vj * hi) / (hi + hj);
+      Xj = Xi;
+    }
+#endif
+    for (int k = 0; k < 3; k++) {
+      /* we add a minus sign since dx is pi->x - pj->x */
+      A[k] = -Xi * (Bi[k][0] * dx[0] + Bi[k][1] * dx[1] + Bi[k][2] * dx[2]) *
+                 wi * hi_inv_dim -
+             Xj * (Bj[k][0] * dx[0] + Bj[k][1] * dx[1] + Bj[k][2] * dx[2]) *
+                 wj * hj_inv_dim;
+      Anorm2 += A[k] * A[k];
+    }
+  } else {
+    /* ill condition gradient matrix: revert to SPH face area */
+    const float Anorm =
+        -(hidp1 * Vi * Vi * wi_dx + hjdp1 * Vj * Vj * wj_dx) * r_inv;
+    A[0] = -Anorm * dx[0];
+    A[1] = -Anorm * dx[1];
+    A[2] = -Anorm * dx[2];
+    Anorm2 = Anorm * Anorm * r2;
+  }
+
+  /* if the interface has no area, nothing happens and we return */
+  /* continuing results in dividing by zero and NaN's... */
+  if (Anorm2 == 0.f) return;
+
+  /* Compute the area */
+  const float Anorm_inv = 1. / sqrtf(Anorm2);
+  const float Anorm = Anorm2 * Anorm_inv;
+
+#ifdef SWIFT_DEBUG_CHECKS
+  /* For stability reasons, we do require A and dx to have opposite
+     directions (basically meaning that the surface normal for the surface
+     always points from particle i to particle j, as it would in a real
+     moving-mesh code). If not, our scheme is no longer upwind and hence can
+     become unstable. */
+  const float dA_dot_dx = A[0] * dx[0] + A[1] * dx[1] + A[2] * dx[2];
+  /* In GIZMO, Phil Hopkins reverts to an SPH integration scheme if this
+     happens. We curently just ignore this case and display a message. */
+  const float rdim = pow_dimension(r);
+  if (dA_dot_dx > 1.e-6f * rdim) {
+    message("Ill conditioned gradient matrix (%g %g %g %g %g)!", dA_dot_dx,
+            Anorm, Vi, Vj, r);
+  }
+#endif
+
+  /* compute the normal vector of the interface */
+  const float n_unit[3] = {A[0] * Anorm_inv, A[1] * Anorm_inv,
+                           A[2] * Anorm_inv};
+
+  /* Compute interface position (relative to pi, since we don't need the actual
+   * position) eqn. (8) */
+  const float xfac = -hi / (hi + hj);
+  const float xij_i[3] = {xfac * dx[0], xfac * dx[1], xfac * dx[2]};
+
+  /* Compute interface velocity */
+  /* eqn. (9) */
+  float xijdotdx = xij_i[0] * dx[0] + xij_i[1] * dx[1] + xij_i[2] * dx[2];
+  xijdotdx *= r_inv * r_inv;
+  const float vij[3] = {vi[0] + (vi[0] - vj[0]) * xijdotdx,
+                        vi[1] + (vi[1] - vj[1]) * xijdotdx,
+                        vi[2] + (vi[2] - vj[2]) * xijdotdx};
+
+  /* complete calculation of position of interface */
+  /* NOTE: dx is not necessarily just pi->x - pj->x but can also contain
+           correction terms for periodicity. If we do the interpolation,
+           we have to use xij w.r.t. the actual particle.
+           => we need a separate xij for pi and pj... */
+  /* tldr: we do not need the code below, but we do need the same code as above
+     but then with i and j swapped */
+  //    for ( k = 0 ; k < 3 ; k++ )
+  //      xij[k] += pi->x[k];
+
+  /* Boost the primitive variables to the frame of reference of the interface */
+  /* Note that velocities are indices 1-3 in W */
+  Wi[1] -= vij[0];
+  Wi[2] -= vij[1];
+  Wi[3] -= vij[2];
+  Wj[1] -= vij[0];
+  Wj[2] -= vij[1];
+  Wj[3] -= vij[2];
+
+  hydro_gradients_predict(pi, pj, hi, hj, dx, r, xij_i, Wi, Wj);
+
+  /* we don't need to rotate, we can use the unit vector in the Riemann problem
+   * itself (see GIZMO) */
+
+  float totflux[5];
+  riemann_solve_for_middle_state_flux(Wi, Wj, n_unit, vij, totflux);
+
+  /* Multiply with the interface surface area */
+  totflux[0] *= Anorm;
+  totflux[1] *= Anorm;
+  totflux[2