Cleaned up some equations in Gizmo's hydro_iact and hydro_gradients. Moved the...

Cleaned up some equations in Gizmo's hydro_iact and hydro_gradients. Moved the gradient prediction to before the velocity boost.

src/hydro/GizmoMFM/hydro_gradients.h

@@ -98,8 +98,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
   /* 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]};
+  const float xij_j[3] = {xij_i[0] + dx[0], xij_i[1] + dx[1], xij_i[2] + dx[2]};
 
   float dWi[5];
   dWi[0] = pi->primitives.gradients.rho[0] * xij_i[0] +

src/hydro/GizmoMFM/hydro_iact.h

@@ -375,11 +375,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
 
   /* 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};
+  const float vij[3] = {vi[0] + (vi[0] - vj[0]) * xfac,
+                        vi[1] + (vi[1] - vj[1]) * xfac,
+                        vi[2] + (vi[2] - vj[2]) * xfac};
 
   /* complete calculation of position of interface */
   /* NOTE: dx is not necessarily just pi->x - pj->x but can also contain
@@ -391,6 +389,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   //    for ( k = 0 ; k < 3 ; k++ )
   //      xij[k] += pi->x[k];
 
+  hydro_gradients_predict(pi, pj, hi, hj, dx, r, xij_i, Wi, Wj);
+
   /* 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];
@@ -400,8 +400,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   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) */
 

src/hydro/GizmoMFV/hydro_gradients.h

@@ -98,8 +98,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
   /* 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]};
+  const float xij_j[3] = {xij_i[0] + dx[0], xij_i[1] + dx[1], xij_i[2] + dx[2]};
 
   float dWi[5];
   dWi[0] = pi->primitives.gradients.rho[0] * xij_i[0] +

src/hydro/GizmoMFV/hydro_iact.h

@@ -375,21 +375,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
 
   /* 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];
+  const float vij[3] = {vi[0] + xfac * (vi[0] - vj[0]),
+                        vi[1] + xfac * (vi[1] - vj[1]),
+                        vi[2] + xfac * (vi[2] - vj[2])};
+
+  hydro_gradients_predict(pi, pj, hi, hj, dx, r, xij_i, Wi, Wj);
 
   /* Boost the primitive variables to the frame of reference of the interface */
   /* Note that velocities are indices 1-3 in W */
@@ -400,8 +390,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   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) */