Put slope limiters in separate files.

src/hydro/Gizmo/hydro_gradients.h

@@ -31,6 +31,8 @@
 #include "hydro_gradients_none.h"
 #endif
 
+#include "hydro_slope_limiters.h"
+
 /**
  * @brief Gradients reconstruction. Is the same for all gradient types (although
  * gradients_none does nothing, since all gradients are zero -- are they?).
@@ -84,69 +86,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
            pj->primitives.gradients.P[1] * xij_j[1] +
            pj->primitives.gradients.P[2] * xij_j[2];
 
-  float xij_i_norm;
-  float phi_i, phi_j;
-  float delta1, delta2;
-  float phiminus, phiplus;
-  float phimin, phimax;
-  float phibar;
-  /* free parameters, values from Hopkins */
-  float psi1 = 0.5, psi2 = 0.25;
-  float phi_mid0, phi_mid;
-
-  for (k = 0; k < 10; k++) {
-    if (k < 5) {
-      phi_i = Wi[k];
-      phi_j = Wj[k];
-      phi_mid0 = Wi[k] + dWi[k];
-      xij_i_norm = sqrtf(xij_i[0] * xij_i[0] + xij_i[1] * xij_i[1] +
-                         xij_i[2] * xij_i[2]);
-    } else {
-      phi_i = Wj[k - 5];
-      phi_j = Wi[k - 5];
-      phi_mid0 = Wj[k - 5] + dWj[k - 5];
-      xij_i_norm = sqrtf(xij_j[0] * xij_j[0] + xij_j[1] * xij_j[1] +
-                         xij_j[2] * xij_j[2]);
-    }
-
-    delta1 = psi1 * fabs(phi_i - phi_j);
-    delta2 = psi2 * fabs(phi_i - phi_j);
-
-    phimin = fmin(phi_i, phi_j);
-    phimax = fmax(phi_i, phi_j);
-
-    phibar = phi_i + xij_i_norm / r * (phi_j - phi_i);
-
-    /* if sign(phimax+delta1) == sign(phimax) */
-    if ((phimax + delta1) * phimax > 0.0f) {
-      phiplus = phimax + delta1;
-    } else {
-      phiplus = phimax / (1.0f + delta1 / fabs(phimax));
-    }
-
-    /* if sign(phimin-delta1) == sign(phimin) */
-    if ((phimin - delta1) * phimin > 0.0f) {
-      phiminus = phimin - delta1;
-    } else {
-      phiminus = phimin / (1.0f + delta1 / fabs(phimin));
-    }
-
-    if (phi_i == phi_j) {
-      phi_mid = phi_i;
-    } else {
-      if (phi_i < phi_j) {
-        phi_mid = fmax(phiminus, fmin(phibar + delta2, phi_mid0));
-      } else {
-        phi_mid = fmin(phiplus, fmax(phibar - delta2, phi_mid0));
-      }
-    }
-
-    if (k < 5) {
-      dWi[k] = phi_mid - phi_i;
-    } else {
-      dWj[k - 5] = phi_mid - phi_i;
-    }
-  }
+  hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);
 
   /* time */
   dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +

src/hydro/Gizmo/hydro_gradients_sph.h

@@ -59,18 +59,7 @@ hydro_gradients_init_density_loop(struct part *p) {
   p->primitives.gradients.P[1] = 0.0f;
   p->primitives.gradients.P[2] = 0.0f;
 
-  p->primitives.limiter.rho[0] = FLT_MAX;
-  p->primitives.limiter.rho[1] = -FLT_MAX;
-  p->primitives.limiter.v[0][0] = FLT_MAX;
-  p->primitives.limiter.v[0][1] = -FLT_MAX;
-  p->primitives.limiter.v[1][0] = FLT_MAX;
-  p->primitives.limiter.v[1][1] = -FLT_MAX;
-  p->primitives.limiter.v[2][0] = FLT_MAX;
-  p->primitives.limiter.v[2][1] = -FLT_MAX;
-  p->primitives.limiter.P[0] = FLT_MAX;
-  p->primitives.limiter.P[1] = -FLT_MAX;
-
-  p->primitives.limiter.maxr = -FLT_MAX;
+  hydro_slope_limit_cell_init(p);
 }
 
 /**
@@ -116,32 +105,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_density_loop(
   pi->primitives.gradients.P[2] -=
       wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
 
-  /* basic slope limiter: collect the maximal and the minimal value for the
-   * primitive variables among the ngbs */
-  pi->primitives.limiter.rho[0] =
-      fmin(pj->primitives.rho, pi->primitives.limiter.rho[0]);
-  pi->primitives.limiter.rho[1] =
-      fmax(pj->primitives.rho, pi->primitives.limiter.rho[1]);
-
-  pi->primitives.limiter.v[0][0] =
-      fmin(pj->primitives.v[0], pi->primitives.limiter.v[0][0]);
-  pi->primitives.limiter.v[0][1] =
-      fmax(pj->primitives.v[0], pi->primitives.limiter.v[0][1]);
-  pi->primitives.limiter.v[1][0] =
-      fmin(pj->primitives.v[1], pi->primitives.limiter.v[1][0]);
-  pi->primitives.limiter.v[1][1] =
-      fmax(pj->primitives.v[1], pi->primitives.limiter.v[1][1]);
-  pi->primitives.limiter.v[2][0] =
-      fmin(pj->primitives.v[2], pi->primitives.limiter.v[2][0]);
-  pi->primitives.limiter.v[2][1] =
-      fmax(pj->primitives.v[2], pi->primitives.limiter.v[2][1]);
-
-  pi->primitives.limiter.P[0] =
-      fmin(pj->primitives.P, pi->primitives.limiter.P[0]);
-  pi->primitives.limiter.P[1] =
-      fmax(pj->primitives.P, pi->primitives.limiter.P[1]);
-
-  pi->primitives.limiter.maxr = fmax(r, pi->primitives.limiter.maxr);
+  hydro_slope_limit_cell_collect(pi, pj, r);
 
   if (mode == 1) {
     /* signs are the same as before, since we swap i and j twice */
@@ -179,31 +143,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_density_loop(
     pj->primitives.gradients.P[2] -=
         wj_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
 
-    /* basic slope limiter: collect the maximal and the minimal value for the
-     * primitive variables among the ngbs */
-    pj->primitives.limiter.rho[0] =
-        fmin(pi->primitives.rho, pj->primitives.limiter.rho[0]);
-    pj->primitives.limiter.rho[1] =
-        fmax(pi->primitives.rho, pj->primitives.limiter.rho[1]);
-    pj->primitives.limiter.v[0][0] =
-        fmin(pi->primitives.v[0], pj->primitives.limiter.v[0][0]);
-    pj->primitives.limiter.v[0][1] =
-        fmax(pi->primitives.v[0], pj->primitives.limiter.v[0][1]);
-    pj->primitives.limiter.v[1][0] =
-        fmin(pi->primitives.v[1], pj->primitives.limiter.v[1][0]);
-    pj->primitives.limiter.v[1][1] =
-        fmax(pi->primitives.v[1], pj->primitives.limiter.v[1][1]);
-    pj->primitives.limiter.v[2][0] =
-        fmin(pi->primitives.v[2], pj->primitives.limiter.v[2][0]);
-    pj->primitives.limiter.v[2][1] =
-        fmax(pi->primitives.v[2], pj->primitives.limiter.v[2][1]);
-
-    pj->primitives.limiter.P[0] =
-        fmin(pi->primitives.P, pj->primitives.limiter.P[0]);
-    pj->primitives.limiter.P[1] =
-        fmax(pi->primitives.P, pj->primitives.limiter.P[1]);
-
-    pj->primitives.limiter.maxr = fmax(r, pj->primitives.limiter.maxr);
+    hydro_slope_limit_cell_collect(pj, pi, r);
   }
 }
 
@@ -213,9 +153,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_density_loop(
 __attribute__((always_inline)) INLINE static void
 hydro_gradients_prepare_force_loop(struct part *p, float ih2, float volume) {
 
-  float gradrho[3], gradv[3][3], gradP[3];
-  float gradtrue, gradmax, gradmin, alpha;
-
   /* finalize gradients by multiplying with volume */
   p->primitives.gradients.rho[0] *= ih2 * ih2 * volume;
   p->primitives.gradients.rho[1] *= ih2 * ih2 * volume;
@@ -237,88 +174,7 @@ hydro_gradients_prepare_force_loop(struct part *p, float ih2, float volume) {
   p->primitives.gradients.P[1] *= ih2 * ih2 * volume;
   p->primitives.gradients.P[2] *= ih2 * ih2 * volume;
 
-  /* slope limiter */
-  gradrho[0] = p->primitives.gradients.rho[0];
-  gradrho[1] = p->primitives.gradients.rho[1];
-  gradrho[2] = p->primitives.gradients.rho[2];
-
-  gradv[0][0] = p->primitives.gradients.v[0][0];
-  gradv[0][1] = p->primitives.gradients.v[0][1];
-  gradv[0][2] = p->primitives.gradients.v[0][2];
-
-  gradv[1][0] = p->primitives.gradients.v[1][0];
-  gradv[1][1] = p->primitives.gradients.v[1][1];
-  gradv[1][2] = p->primitives.gradients.v[1][2];
-
-  gradv[2][0] = p->primitives.gradients.v[2][0];
-  gradv[2][1] = p->primitives.gradients.v[2][1];
-  gradv[2][2] = p->primitives.gradients.v[2][2];
-
-  gradP[0] = p->primitives.gradients.P[0];
-  gradP[1] = p->primitives.gradients.P[1];
-  gradP[2] = p->primitives.gradients.P[2];
-
-  gradtrue = sqrtf(gradrho[0] * gradrho[0] + gradrho[1] * gradrho[1] +
-                   gradrho[2] * gradrho[2]); /* gradtrue might be zero. In this
-    case, there is no gradient and we don't
-    need to slope limit anything... */
-  if (gradtrue) {
-    gradtrue *= p->primitives.limiter.maxr;
-    gradmax = p->primitives.limiter.rho[1] - p->primitives.rho;
-    gradmin = p->primitives.rho - p->primitives.limiter.rho[0];
-    alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
-    p->primitives.gradients.rho[0] *= alpha;
-    p->primitives.gradients.rho[1] *= alpha;
-    p->primitives.gradients.rho[2] *= alpha;
-  }
-
-  gradtrue = sqrtf(gradv[0][0] * gradv[0][0] + gradv[0][1] * gradv[0][1] +
-                   gradv[0][2] * gradv[0][2]);
-  if (gradtrue) {
-    gradtrue *= p->primitives.limiter.maxr;
-    gradmax = p->primitives.limiter.v[0][1] - p->primitives.v[0];
-    gradmin = p->primitives.v[0] - p->primitives.limiter.v[0][0];
-    alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
-    p->primitives.gradients.v[0][0] *= alpha;
-    p->primitives.gradients.v[0][1] *= alpha;
-    p->primitives.gradients.v[0][2] *= alpha;
-  }
-
-  gradtrue = sqrtf(gradv[1][0] * gradv[1][0] + gradv[1][1] * gradv[1][1] +
-                   gradv[1][2] * gradv[1][2]);
-  if (gradtrue) {
-    gradtrue *= p->primitives.limiter.maxr;
-    gradmax = p->primitives.limiter.v[1][1] - p->primitives.v[1];
-    gradmin = p->primitives.v[1] - p->primitives.limiter.v[1][0];
-    alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
-    p->primitives.gradients.v[1][0] *= alpha;
-    p->primitives.gradients.v[1][1] *= alpha;
-    p->primitives.gradients.v[1][2] *= alpha;
-  }
-
-  gradtrue = sqrtf(gradv[2][0] * gradv[2][0] + gradv[2][1] * gradv[2][1] +
-                   gradv[2][2] * gradv[2][2]);
-  if (gradtrue) {
-    gradtrue *= p->primitives.limiter.maxr;
-    gradmax = p->primitives.limiter.v[2][1] - p->primitives.v[2];
-    gradmin = p->primitives.v[2] - p->primitives.limiter.v[2][0];
-    alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
-    p->primitives.gradients.v[2][0] *= alpha;
-    p->primitives.gradients.v[2][1] *= alpha;
-    p->primitives.gradients.v[2][2] *= alpha;
-  }
-
-  gradtrue =
-      sqrtf(gradP[0] * gradP[0] + gradP[1] * gradP[1] + gradP[2] * gradP[2]);
-  if (gradtrue) {
-    gradtrue *= p->primitives.limiter.maxr;
-    gradmax = p->primitives.limiter.P[1] - p->primitives.P;
-    gradmin = p->primitives.P - p->primitives.limiter.P[0];
-    alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
-    p->primitives.gradients.P[0] *= alpha;
-    p->primitives.gradients.P[1] *= alpha;
-    p->primitives.gradients.P[2] *= alpha;
-  }
+  hydro_slope_limit_cell(p);
 }
 
 /**

src/hydro/Gizmo/hydro_slope_limiters.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_SLOPE_LIMITERS_H
+#define SWIFT_HYDRO_SLOPE_LIMITERS_H
+
+//#define PER_FACE_LIMITER
+//#define CELL_WIDE_LIMITER
+
+#ifdef PER_FACE_LIMITER
+#include "hydro_slope_limiters_face.h"
+#else
+
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_face(
+    float *Wi, float *Wj, float *dWi, float *dWj, float *xij_i, float *xij_j,
+    float r) {}
+
+#endif
+
+#ifdef CELL_WIDE_LIMITER
+#include "hydro_slope_limiters_cell.h"
+#else
+
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell_init(
+    struct part *p) {}
+
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell_collect(
+    struct part* pi, struct part* pj, float r) {}
+
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
+    struct part *p) {}
+
+#endif
+
+#endif  // SWIFT_HYDRO_SLOPE_LIMITERS_H

src/hydro/Gizmo/hydro_slope_limiters_cell.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/>.
+ *
+ ******************************************************************************/
+
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell_init(
+    struct part* p) {
+
+  p->primitives.limiter.rho[0] = FLT_MAX;
+  p->primitives.limiter.rho[1] = -FLT_MAX;
+  p->primitives.limiter.v[0][0] = FLT_MAX;
+  p->primitives.limiter.v[0][1] = -FLT_MAX;
+  p->primitives.limiter.v[1][0] = FLT_MAX;
+  p->primitives.limiter.v[1][1] = -FLT_MAX;
+  p->primitives.limiter.v[2][0] = FLT_MAX;
+  p->primitives.limiter.v[2][1] = -FLT_MAX;
+  p->primitives.limiter.P[0] = FLT_MAX;
+  p->primitives.limiter.P[1] = -FLT_MAX;
+
+  p->primitives.limiter.maxr = -FLT_MAX;
+}
+
+__attribute__((always_inline)) INLINE static void
+hydro_slope_limit_cell_collect(struct part* pi, struct part* pj, float r) {
+
+  /* basic slope limiter: collect the maximal and the minimal value for the
+   * primitive variables among the ngbs */
+  pi->primitives.limiter.rho[0] =
+      fmin(pj->primitives.rho, pi->primitives.limiter.rho[0]);
+  pi->primitives.limiter.rho[1] =
+      fmax(pj->primitives.rho, pi->primitives.limiter.rho[1]);
+
+  pi->primitives.limiter.v[0][0] =
+      fmin(pj->primitives.v[0], pi->primitives.limiter.v[0][0]);
+  pi->primitives.limiter.v[0][1] =
+      fmax(pj->primitives.v[0], pi->primitives.limiter.v[0][1]);
+  pi->primitives.limiter.v[1][0] =
+      fmin(pj->primitives.v[1], pi->primitives.limiter.v[1][0]);
+  pi->primitives.limiter.v[1][1] =
+      fmax(pj->primitives.v[1], pi->primitives.limiter.v[1][1]);
+  pi->primitives.limiter.v[2][0] =
+      fmin(pj->primitives.v[2], pi->primitives.limiter.v[2][0]);
+  pi->primitives.limiter.v[2][1] =
+      fmax(pj->primitives.v[2], pi->primitives.limiter.v[2][1]);
+
+  pi->primitives.limiter.P[0] =
+      fmin(pj->primitives.P, pi->primitives.limiter.P[0]);
+  pi->primitives.limiter.P[1] =
+      fmax(pj->primitives.P, pi->primitives.limiter.P[1]);
+
+  pi->primitives.limiter.maxr = fmax(r, pi->primitives.limiter.maxr);
+}
+
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
+    struct part* p) {
+
+  float gradrho[3], gradv[3][3], gradP[3];
+  float gradtrue, gradmax, gradmin, alpha;
+
+  gradrho[0] = p->primitives.gradients.rho[0];
+  gradrho[1] = p->primitives.gradients.rho[1];
+  gradrho[2] = p->primitives.gradients.rho[2];
+
+  gradv[0][0] = p->primitives.gradients.v[0][0];
+  gradv[0][1] = p->primitives.gradients.v[0][1];
+  gradv[0][2] = p->primitives.gradients.v[0][2];
+
+  gradv[1][0] = p->primitives.gradients.v[1][0];
+  gradv[1][1] = p->primitives.gradients.v[1][1];
+  gradv[1][2] = p->primitives.gradients.v[1][2];
+
+  gradv[2][0] = p->primitives.gradients.v[2][0];
+  gradv[2][1] = p->primitives.gradients.v[2][1];
+  gradv[2][2] = p->primitives.gradients.v[2][2];
+
+  gradP[0] = p->primitives.gradients.P[0];
+  gradP[1] = p->primitives.gradients.P[1];
+  gradP[2] = p->primitives.gradients.P[2];
+
+  gradtrue = sqrtf(gradrho[0] * gradrho[0] + gradrho[1] * gradrho[1] +
+                   gradrho[2] * gradrho[2]);
+  if (gradtrue) {
+    gradtrue *= p->primitives.limiter.maxr;
+    gradmax = p->primitives.limiter.rho[1] - p->primitives.rho;
+    gradmin = p->primitives.rho - p->primitives.limiter.rho[0];
+    alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
+    p->primitives.gradients.rho[0] *= alpha;
+    p->primitives.gradients.rho[1] *= alpha;
+    p->primitives.gradients.rho[2] *= alpha;
+  }
+
+  gradtrue = sqrtf(gradv[0][0] * gradv[0][0] + gradv[0][1] * gradv[0][1] +
+                   gradv[0][2] * gradv[0][2]);
+  if (gradtrue) {
+    gradtrue *= p->primitives.limiter.maxr;
+    gradmax = p->primitives.limiter.v[0][1] - p->primitives.v[0];
+    gradmin = p->primitives.v[0] - p->primitives.limiter.v[0][0];
+    alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
+    p->primitives.gradients.v[0][0] *= alpha;
+    p->primitives.gradients.v[0][1] *= alpha;
+    p->primitives.gradients.v[0][2] *= alpha;
+  }
+
+  gradtrue = sqrtf(gradv[1][0] * gradv[1][0] + gradv[1][1] * gradv[1][1] +
+                   gradv[1][2] * gradv[1][2]);
+  if (gradtrue) {
+    gradtrue *= p->primitives.limiter.maxr;
+    gradmax = p->primitives.limiter.v[1][1] - p->primitives.v[1];
+    gradmin = p->primitives.v[1] - p->primitives.limiter.v[1][0];
+    alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
+    p->primitives.gradients.v[1][0] *= alpha;
+    p->primitives.gradients.v[1][1] *= alpha;
+    p->primitives.gradients.v[1][2] *= alpha;
+  }
+
+  gradtrue = sqrtf(gradv[2][0] * gradv[2][0] + gradv[2][1] * gradv[2][1] +
+                   gradv[2][2] * gradv[2][2]);
+  if (gradtrue) {
+    gradtrue *= p->primitives.limiter.maxr;
+    gradmax = p->primitives.limiter.v[2][1] - p->primitives.v[2];
+    gradmin = p->primitives.v[2] - p->primitives.limiter.v[2][0];
+    alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
+    p->primitives.gradients.v[2][0] *= alpha;
+    p->primitives.gradients.v[2][1] *= alpha;
+    p->primitives.gradients.v[2][2] *= alpha;
+  }
+
+  gradtrue =
+      sqrtf(gradP[0] * gradP[0] + gradP[1] * gradP[1] + gradP[2] * gradP[2]);
+  if (gradtrue) {
+    gradtrue *= p->primitives.limiter.maxr;
+    gradmax = p->primitives.limiter.P[1] - p->primitives.P;
+    gradmin = p->primitives.P - p->primitives.limiter.P[0];
+    alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
+    p->primitives.gradients.P[0] *= alpha;
+    p->primitives.gradients.P[1] *= alpha;
+    p->primitives.gradients.P[2] *= alpha;
+  }
+}

src/hydro/Gizmo/hydro_slope_limiters_face.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/>.
+ *
+ ******************************************************************************/
+
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_face(
+    float *Wi, float *Wj, float *dWi, float *dWj, float *xij_i, float *xij_j,
+    float r) {
+
+  float xij_i_norm;
+  float phi_i, phi_j;
+  float delta1, delta2;
+  float phiminus, phiplus;
+  float phimin, phimax;
+  float phibar;
+  /* free parameters, values from Hopkins */
+  float psi1 = 0.5, psi2 = 0.25;
+  float phi_mid0, phi_mid;
+  int k;
+
+  for (k = 0; k < 10; k++) {
+    if (k < 5) {
+      phi_i = Wi[k];
+      phi_j = Wj[k];
+      phi_mid0 = Wi[k] + dWi[k];
+      xij_i_norm = sqrtf(xij_i[0] * xij_i[0] + xij_i[1] * xij_i[1] +
+                         xij_i[2] * xij_i[2]);
+    } else {
+      phi_i = Wj[k - 5];
+      phi_j = Wi[k - 5];
+      phi_mid0 = Wj[k - 5] + dWj[k - 5];
+      xij_i_norm = sqrtf(xij_j[0] * xij_j[0] + xij_j[1] * xij_j[1] +
+                         xij_j[2] * xij_j[2]);
+    }
+
+    delta1 = psi1 * fabs(phi_i - phi_j);
+    delta2 = psi2 * fabs(phi_i - phi_j);
+
+    phimin = fmin(phi_i, phi_j);
+    phimax = fmax(phi_i, phi_j);
+
+    phibar = phi_i + xij_i_norm / r * (phi_j - phi_i);
+
+    /* if sign(phimax+delta1) == sign(phimax) */
+    if ((phimax + delta1) * phimax > 0.0f) {
+      phiplus = phimax + delta1;
+    } else {
+      phiplus = phimax / (1.0f + delta1 / fabs(phimax));
+    }
+
+    /* if sign(phimin-delta1) == sign(phimin) */
+    if ((phimin - delta1) * phimin > 0.0f) {
+      phiminus = phimin - delta1;
+    } else {
+      phiminus = phimin / (1.0f + delta1 / fabs(phimin));
+    }
+
+    if (phi_i == phi_j) {
+      phi_mid = phi_i;
+    } else {
+      if (phi_i < phi_j) {
+        phi_mid = fmax(phiminus, fmin(phibar + delta2, phi_mid0));
+      } else {
+        phi_mid = fmin(phiplus, fmax(phibar - delta2, phi_mid0));
+      }
+    }
+
+    if (k < 5) {
+      dWi[k] = phi_mid - phi_i;
+    } else {
+      dWj[k - 5] = phi_mid - phi_i;
+    }
+  }
+}