Merge branch 'gizmo_mfv_clean' into 'master'

Gizmo mfv clean

See merge request !587

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.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);

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) */
 
@@ -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];

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". */