diff --git a/src/hydro/Gizmo/hydro_gradients.h b/src/hydro/Gizmo/hydro_gradients.h
index 0f9ea0798f990a082f29269d4c7335817d13bc29..012e1de43b50091cdf45f280459dedf24c172d62 100644
--- a/src/hydro/Gizmo/hydro_gradients.h
+++ b/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] +
diff --git a/src/hydro/Gizmo/hydro_gradients_sph.h b/src/hydro/Gizmo/hydro_gradients_sph.h
index 71ceba532d1c321d65dbe77dbb1925faf5dbd8c0..8c9d7be052920c9fc86216c63b55179c4dfd2153 100644
--- a/src/hydro/Gizmo/hydro_gradients_sph.h
+++ b/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);
}
/**
diff --git a/src/hydro/Gizmo/hydro_slope_limiters.h b/src/hydro/Gizmo/hydro_slope_limiters.h
new file mode 100644
index 0000000000000000000000000000000000000000..bcb370fd0076e9e2092ea0c6840304f77bc74e06
--- /dev/null
+++ b/src/hydro/Gizmo/hydro_slope_limiters.h
@@ -0,0 +1,51 @@
+/*******************************************************************************
+ * 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
diff --git a/src/hydro/Gizmo/hydro_slope_limiters_cell.h b/src/hydro/Gizmo/hydro_slope_limiters_cell.h
new file mode 100644
index 0000000000000000000000000000000000000000..32df7c53625eea90f1c7f97fcf36c6c2c3f9be4e
--- /dev/null
+++ b/src/hydro/Gizmo/hydro_slope_limiters_cell.h
@@ -0,0 +1,153 @@
+/*******************************************************************************
+ * 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;
+ }
+}
diff --git a/src/hydro/Gizmo/hydro_slope_limiters_face.h b/src/hydro/Gizmo/hydro_slope_limiters_face.h
new file mode 100644
index 0000000000000000000000000000000000000000..7ef9319877fbb60309afdf7b92acf296f9e117b4
--- /dev/null
+++ b/src/hydro/Gizmo/hydro_slope_limiters_face.h
@@ -0,0 +1,88 @@
+/*******************************************************************************
+ * 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;
+ }
+ }
+}