diff --git a/src/hydro/GizmoMFM/hydro_gradients.h b/src/hydro/GizmoMFM/hydro_gradients.h
index 964a2adcfe09b95c2a221af540e5e3ff0830dd67..85617f05b6a5350e7b1b60de53f2c999b6ba7453 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 a1a82e514baa60d5895343ea84ef1c3aedc8a6b7..fb54202a86b3e83aad7b2759c18dd3ce548cda13 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.h b/src/hydro/GizmoMFV/hydro.h
index 6916fe33272692316354385b723ce9969606b6a2..c65f2bb5899b63dbaea4db2f108cf3d914ca78f0 100644
--- a/src/hydro/GizmoMFV/hydro.h
+++ b/src/hydro/GizmoMFV/hydro.h
@@ -591,12 +591,28 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->conserved.flux.momentum[2] * dt_drift / p->conserved.mass;
#if !defined(EOS_ISOTHERMAL_GAS)
- const float u = p->conserved.energy + p->conserved.flux.energy * dt_therm;
- p->primitives.P =
- hydro_gamma_minus_one * u * p->primitives.rho / p->conserved.mass;
+#ifdef GIZMO_TOTAL_ENERGY
+ const float Etot =
+ p->conserved.energy + p->conserved.flux.energy * dt_therm;
+ const float v2 = (p->primitives.v[0] * p->primitives.v[0] +
+ p->primitives.v[1] * p->primitives.v[1] +
+ p->primitives.v[2] * p->primitives.v[2]);
+ const float u = (Etot / p->conserved.mass - 0.5 * v2);
+#else
+ const float u =
+ (p->conserved.energy + p->conserved.flux.energy * dt_therm) /
+ p->conserved.mass;
+#endif
+ p->primitives.P = hydro_gamma_minus_one * u * p->primitives.rho;
#endif
}
+ /* we use a sneaky way to get the gravitational contribtuion to the
+ velocity update */
+ p->primitives.v[0] += p->v[0] - xp->v_full[0];
+ p->primitives.v[1] += p->v[1] - xp->v_full[1];
+ p->primitives.v[2] += p->v[2] - xp->v_full[2];
+
#ifdef SWIFT_DEBUG_CHECKS
if (p->h <= 0.) {
error("Zero or negative smoothing length (%g)!", p->h);
@@ -646,6 +662,34 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
float a_grav[3];
+ /* Add gravity. We only do this if we have gravity activated. */
+ if (p->gpart) {
+ /* Retrieve the current value of the gravitational acceleration from the
+ gpart. We are only allowed to do this because this is the kick. We still
+ need to check whether gpart exists though.*/
+ a_grav[0] = p->gpart->a_grav[0];
+ a_grav[1] = p->gpart->a_grav[1];
+ a_grav[2] = p->gpart->a_grav[2];
+
+#ifdef GIZMO_TOTAL_ENERGY
+ p->conserved.energy += dt * (p->conserved.momentum[0] * a_grav[0] +
+ p->conserved.momentum[1] * a_grav[1] +
+ p->conserved.momentum[2] * a_grav[2]);
+#endif
+
+ /* Kick the momentum for half a time step */
+ /* Note that this also affects the particle movement, as the velocity for
+ the particles is set after this. */
+ p->conserved.momentum[0] += p->conserved.mass * a_grav[0] * dt;
+ p->conserved.momentum[1] += p->conserved.mass * a_grav[1] * dt;
+ p->conserved.momentum[2] += p->conserved.mass * a_grav[2] * dt;
+
+ p->conserved.energy -=
+ 0.5f * dt *
+ (p->gravity.mflux[0] * a_grav[0] + p->gravity.mflux[1] * a_grav[1] +
+ p->gravity.mflux[2] * a_grav[2]);
+ }
+
/* Update conserved variables. */
p->conserved.mass += p->conserved.flux.mass * dt;
p->conserved.momentum[0] += p->conserved.flux.momentum[0] * dt;
@@ -688,28 +732,9 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
}
#endif
- /* Add gravity. We only do this if we have gravity activated. */
if (p->gpart) {
- /* Retrieve the current value of the gravitational acceleration from the
- gpart. We are only allowed to do this because this is the kick. We still
- need to check whether gpart exists though.*/
- a_grav[0] = p->gpart->a_grav[0];
- a_grav[1] = p->gpart->a_grav[1];
- a_grav[2] = p->gpart->a_grav[2];
-
/* Make sure the gpart knows the mass has changed. */
p->gpart->mass = p->conserved.mass;
-
- /* Kick the momentum for half a time step */
- /* Note that this also affects the particle movement, as the velocity for
- the particles is set after this. */
- p->conserved.momentum[0] += dt * p->conserved.mass * a_grav[0];
- p->conserved.momentum[1] += dt * p->conserved.mass * a_grav[1];
- p->conserved.momentum[2] += dt * p->conserved.mass * a_grav[2];
-
- p->conserved.energy += dt * (p->gravity.mflux[0] * a_grav[0] +
- p->gravity.mflux[1] * a_grav[1] +
- p->gravity.mflux[2] * a_grav[2]);
}
hydro_velocities_set(p, xp);
diff --git a/src/hydro/GizmoMFV/hydro_gradients.h b/src/hydro/GizmoMFV/hydro_gradients.h
index 387c263775ddfb37e4c7cb31a624ba5dc673beb2..4046e121bad9e329fecc30afb435bffca2815346 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 bb835094acd285b109383c1f5d04a6f5e2d936df..c766ce3cc9048f8da8b3438c3c27e6998dd5df7e 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) */
@@ -447,9 +435,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
if (mode == 1) {
/* Store mass flux */
const float mflux_j = totflux[0];
- pj->gravity.mflux[0] -= mflux_j * dx[0];
- pj->gravity.mflux[1] -= mflux_j * dx[1];
- pj->gravity.mflux[2] -= mflux_j * dx[2];
+ pj->gravity.mflux[0] += mflux_j * dx[0];
+ pj->gravity.mflux[1] += mflux_j * dx[1];
+ pj->gravity.mflux[2] += mflux_j * dx[2];
pj->conserved.flux.mass += totflux[0];
pj->conserved.flux.momentum[0] += totflux[1];
diff --git a/src/hydro/GizmoMFV/hydro_velocities.h b/src/hydro/GizmoMFV/hydro_velocities.h
index ea30d5a6c9c74d34ffd73fa6ab941640f37b02e4..a61a482683291cfc0662116ada99b94a4ccfe68d 100644
--- a/src/hydro/GizmoMFV/hydro_velocities.h
+++ b/src/hydro/GizmoMFV/hydro_velocities.h
@@ -103,7 +103,6 @@ __attribute__((always_inline)) INLINE static void hydro_velocities_set(
p->v[2] = p->conserved.momentum[2] * inverse_mass;
#ifdef GIZMO_STEER_MOTION
-
/* Add a correction to the velocity to keep particle positions close enough
to
the centroid of their mesh-free "cell". */