From bc8abdbe28831212df2aae1e513fa95d2a850f8c Mon Sep 17 00:00:00 2001
From: Bert Vandenbroucke <bert.vandenbroucke@gmail.com>
Date: Fri, 25 May 2018 10:34:21 +0100
Subject: [PATCH] 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 | 3 +--
src/hydro/GizmoMFM/hydro_iact.h | 12 +++++-------
src/hydro/GizmoMFV/hydro_gradients.h | 3 +--
src/hydro/GizmoMFV/hydro_iact.h | 22 +++++-----------------
4 files changed, 12 insertions(+), 28 deletions(-)
diff --git a/src/hydro/GizmoMFM/hydro_gradients.h b/src/hydro/GizmoMFM/hydro_gradients.h
index 964a2adcfe..85617f05b6 100644
--- a/src/hydro/GizmoMFM/hydro_gradients.h
+++ b/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] +
diff --git a/src/hydro/GizmoMFM/hydro_iact.h b/src/hydro/GizmoMFM/hydro_iact.h
index a1a82e514b..fb54202a86 100644
--- a/src/hydro/GizmoMFM/hydro_iact.h
+++ b/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) */
diff --git a/src/hydro/GizmoMFV/hydro_gradients.h b/src/hydro/GizmoMFV/hydro_gradients.h
index 387c263775..4046e121ba 100644
--- a/src/hydro/GizmoMFV/hydro_gradients.h
+++ b/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] +
diff --git a/src/hydro/GizmoMFV/hydro_iact.h b/src/hydro/GizmoMFV/hydro_iact.h
index bb835094ac..61a8df34f4 100644
--- a/src/hydro/GizmoMFV/hydro_iact.h
+++ b/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) */
--
GitLab