From f42f43f11b4f0d90007c9dd409bb13c54a1c84d7 Mon Sep 17 00:00:00 2001
From: Bert Vandenbroucke <bert.vandenbroucke@ugent.be>
Date: Wed, 10 Aug 2016 11:52:31 +0100
Subject: [PATCH] Extracted gradient reconstruction from hydro_iact as well.

---
 src/hydro/Gizmo/hydro_gradients.h | 191 +++++++++++++++++++++++++++-
 src/hydro/Gizmo/hydro_iact.h      | 199 +-----------------------------
 2 files changed, 192 insertions(+), 198 deletions(-)

diff --git a/src/hydro/Gizmo/hydro_gradients.h b/src/hydro/Gizmo/hydro_gradients.h
index a0054c7d5a..0f9ea0798f 100644
--- a/src/hydro/Gizmo/hydro_gradients.h
+++ b/src/hydro/Gizmo/hydro_gradients.h
@@ -20,7 +20,7 @@
 #ifndef SWIFT_HYDRO_GRADIENTS_H
 #define SWIFT_HYDRO_GRADIENTS_H
 
-#define SPH_GRADIENTS
+//#define SPH_GRADIENTS
 
 #if defined(SPH_GRADIENTS)
 #include "hydro_gradients_sph.h"
@@ -31,4 +31,193 @@
 #include "hydro_gradients_none.h"
 #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* pi, struct part* pj, float hi, float hj, float* dx, float r,
+    float* xij_i, float* Wi, float* Wj, float mindt) {
+
+  float dWi[5], dWj[5];
+  float xij_j[3];
+  int k;
+  float xfac;
+
+  /* perform gradient reconstruction in space and time */
+  /* space */
+  /* Compute interface position (relative to pj, since we don't need the actual
+   * position) */
+  /* eqn. (8) */
+  xfac = hj / (hi + hj);
+  for (k = 0; k < 3; k++) xij_j[k] = xfac * dx[k];
+
+  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];
+
+  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];
+
+  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;
+    }
+  }
+
+  /* time */
+  dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
+                           Wi[2] * pi->primitives.gradients.rho[1] +
+                           Wi[3] * pi->primitives.gradients.rho[2] +
+                           Wi[0] * (pi->primitives.gradients.v[0][0] +
+                                    pi->primitives.gradients.v[1][1] +
+                                    pi->primitives.gradients.v[2][2]));
+  dWi[1] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[0][0] +
+                           Wi[2] * pi->primitives.gradients.v[0][1] +
+                           Wi[3] * pi->primitives.gradients.v[0][2] +
+                           pi->primitives.gradients.P[0] / Wi[0]);
+  dWi[2] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[1][0] +
+                           Wi[2] * pi->primitives.gradients.v[1][1] +
+                           Wi[3] * pi->primitives.gradients.v[1][2] +
+                           pi->primitives.gradients.P[1] / Wi[0]);
+  dWi[3] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[2][0] +
+                           Wi[2] * pi->primitives.gradients.v[2][1] +
+                           Wi[3] * pi->primitives.gradients.v[2][2] +
+                           pi->primitives.gradients.P[2] / Wi[0]);
+  dWi[4] -=
+      0.5 * mindt * (Wi[1] * pi->primitives.gradients.P[0] +
+                     Wi[2] * pi->primitives.gradients.P[1] +
+                     Wi[3] * pi->primitives.gradients.P[2] +
+                     hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
+                                            pi->primitives.gradients.v[1][1] +
+                                            pi->primitives.gradients.v[2][2]));
+
+  dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
+                           Wj[2] * pj->primitives.gradients.rho[1] +
+                           Wj[3] * pj->primitives.gradients.rho[2] +
+                           Wj[0] * (pj->primitives.gradients.v[0][0] +
+                                    pj->primitives.gradients.v[1][1] +
+                                    pj->primitives.gradients.v[2][2]));
+  dWj[1] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[0][0] +
+                           Wj[2] * pj->primitives.gradients.v[0][1] +
+                           Wj[3] * pj->primitives.gradients.v[0][2] +
+                           pj->primitives.gradients.P[0] / Wj[0]);
+  dWj[2] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[1][0] +
+                           Wj[2] * pj->primitives.gradients.v[1][1] +
+                           Wj[3] * pj->primitives.gradients.v[1][2] +
+                           pj->primitives.gradients.P[1] / Wj[0]);
+  dWj[3] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[2][0] +
+                           Wj[2] * pj->primitives.gradients.v[2][1] +
+                           Wj[3] * pj->primitives.gradients.v[2][2] +
+                           pj->primitives.gradients.P[2] / Wj[0]);
+  dWj[4] -=
+      0.5 * mindt * (Wj[1] * pj->primitives.gradients.P[0] +
+                     Wj[2] * pj->primitives.gradients.P[1] +
+                     Wj[3] * pj->primitives.gradients.P[2] +
+                     hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
+                                            pj->primitives.gradients.v[1][1] +
+                                            pj->primitives.gradients.v[2][2]));
+
+  //    printf("WL: %g %g %g %g %g\n", Wi[0], Wi[1], Wi[2], Wi[3], Wi[4]);
+  //    printf("WR: %g %g %g %g %g\n", Wj[0], Wj[1], Wj[2], Wj[3], Wj[4]);
+
+  //    printf("dWL: %g %g %g %g %g\n", dWi[0], dWi[1], dWi[2], dWi[3], dWi[4]);
+  //    printf("dWR: %g %g %g %g %g\n", dWj[0], dWj[1], dWj[2], dWj[3], dWj[4]);
+
+  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];
+}
+
 #endif  // SWIFT_HYDRO_GRADIENTS_H
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index 4d0d0683aa..4e6965d5bb 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -125,13 +125,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   float xij_i[3], xfac, xijdotdx;
   float vmax, dvdotdx;
   float vi[3], vj[3], vij[3];
-  float Wi[5], Wj[5];  //, Whalf[5];
-#ifdef USE_GRADIENTS
-  float dWi[5], dWj[5];
-  float xij_j[3];
-#endif
-  //    float rhoe;
-  //    float flux[5][3];
+  float Wi[5], Wj[5];
   float dti, dtj, mindt;
   float n_unit[3];
 
@@ -159,10 +153,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   dti = pi->ti_end - pi->ti_begin;  // MATTHIEU
   dtj = pj->ti_end - pj->ti_begin;
 
-  //    if(dti > 1.e-7 || dtj > 1.e-7){
-  //        message("Timestep too large: %g %g!", dti, dtj);
-  //    }
-
   /* calculate the maximal signal velocity */
   vmax =
       sqrtf(hydro_gamma * Wi[4] / Wi[0]) + sqrtf(hydro_gamma * Wj[4] / Wj[0]);
@@ -252,192 +242,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   Wj[2] -= vij[1];
   Wj[3] -= vij[2];
 
-#ifdef USE_GRADIENTS
-  /* perform gradient reconstruction in space and time */
-  /* space */
-  /* Compute interface position (relative to pj, since we don't need the actual
-   * position) */
-  /* eqn. (8) */
-  xfac = hj / (hi + hj);
-  for (k = 0; k < 3; k++) xij_j[k] = xfac * dx[k];
-
-  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];
-
-  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];
-
-#ifdef PER_FACE_LIMITER
-
-  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;
-    }
-  }
-
-#endif
-
-  //    printf("dWL: %g %g %g %g %g\n", dWi[0], dWi[1], dWi[2], dWi[3], dWi[4]);
-  //    printf("dWR: %g %g %g %g %g\n", dWj[0], dWj[1], dWj[2], dWj[3], dWj[4]);
-
-  /* time */
-  dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
-                           Wi[2] * pi->primitives.gradients.rho[1] +
-                           Wi[3] * pi->primitives.gradients.rho[2] +
-                           Wi[0] * (pi->primitives.gradients.v[0][0] +
-                                    pi->primitives.gradients.v[1][1] +
-                                    pi->primitives.gradients.v[2][2]));
-  dWi[1] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[0][0] +
-                           Wi[2] * pi->primitives.gradients.v[0][1] +
-                           Wi[3] * pi->primitives.gradients.v[0][2] +
-                           pi->primitives.gradients.P[0] / Wi[0]);
-  dWi[2] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[1][0] +
-                           Wi[2] * pi->primitives.gradients.v[1][1] +
-                           Wi[3] * pi->primitives.gradients.v[1][2] +
-                           pi->primitives.gradients.P[1] / Wi[0]);
-  dWi[3] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[2][0] +
-                           Wi[2] * pi->primitives.gradients.v[2][1] +
-                           Wi[3] * pi->primitives.gradients.v[2][2] +
-                           pi->primitives.gradients.P[2] / Wi[0]);
-  dWi[4] -=
-      0.5 * mindt * (Wi[1] * pi->primitives.gradients.P[0] +
-                     Wi[2] * pi->primitives.gradients.P[1] +
-                     Wi[3] * pi->primitives.gradients.P[2] +
-                     hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
-                                            pi->primitives.gradients.v[1][1] +
-                                            pi->primitives.gradients.v[2][2]));
-
-  dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
-                           Wj[2] * pj->primitives.gradients.rho[1] +
-                           Wj[3] * pj->primitives.gradients.rho[2] +
-                           Wj[0] * (pj->primitives.gradients.v[0][0] +
-                                    pj->primitives.gradients.v[1][1] +
-                                    pj->primitives.gradients.v[2][2]));
-  dWj[1] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[0][0] +
-                           Wj[2] * pj->primitives.gradients.v[0][1] +
-                           Wj[3] * pj->primitives.gradients.v[0][2] +
-                           pj->primitives.gradients.P[0] / Wj[0]);
-  dWj[2] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[1][0] +
-                           Wj[2] * pj->primitives.gradients.v[1][1] +
-                           Wj[3] * pj->primitives.gradients.v[1][2] +
-                           pj->primitives.gradients.P[1] / Wj[0]);
-  dWj[3] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[2][0] +
-                           Wj[2] * pj->primitives.gradients.v[2][1] +
-                           Wj[3] * pj->primitives.gradients.v[2][2] +
-                           pj->primitives.gradients.P[2] / Wj[0]);
-  dWj[4] -=
-      0.5 * mindt * (Wj[1] * pj->primitives.gradients.P[0] +
-                     Wj[2] * pj->primitives.gradients.P[1] +
-                     Wj[3] * pj->primitives.gradients.P[2] +
-                     hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
-                                            pj->primitives.gradients.v[1][1] +
-                                            pj->primitives.gradients.v[2][2]));
-
-  //    printf("WL: %g %g %g %g %g\n", Wi[0], Wi[1], Wi[2], Wi[3], Wi[4]);
-  //    printf("WR: %g %g %g %g %g\n", Wj[0], Wj[1], Wj[2], Wj[3], Wj[4]);
-
-  //    printf("dWL: %g %g %g %g %g\n", dWi[0], dWi[1], dWi[2], dWi[3], dWi[4]);
-  //    printf("dWR: %g %g %g %g %g\n", dWj[0], dWj[1], dWj[2], dWj[3], dWj[4]);
-
-  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];
-#endif
-
-  /* apply slope limiter interface by interface */
-  /* ... to be done ... */
+  hydro_gradients_predict(pi, pj, hi, hj, dx, r, xij_i, Wi, Wj, mindt);
 
   /* we don't need to rotate, we can use the unit vector in the Riemann problem
    * itself (see GIZMO) */
-- 
GitLab