From 191c2cf297d7b00ec1e13f54255091d993f925fd Mon Sep 17 00:00:00 2001
From: Yolan Uyttenhove <yolan.uyttenhove@ugent.be>
Date: Wed, 31 May 2023 15:41:40 +0000
Subject: [PATCH] Use correct fluid velocity in time extrapolation (Gizmo).

---
 src/hydro/Gizmo/MFM/hydro_flux.h | 10 ++-------
 src/hydro/Gizmo/MFV/hydro_flux.h |  6 ------
 src/hydro/Gizmo/hydro.h          | 37 +++++++++++++++++++-------------
 src/hydro/Gizmo/hydro_getters.h  | 22 +++++++++++++++++--
 src/hydro/Gizmo/hydro_setters.h  | 16 ++++++++++++++
 5 files changed, 60 insertions(+), 31 deletions(-)

diff --git a/src/hydro/Gizmo/MFM/hydro_flux.h b/src/hydro/Gizmo/MFM/hydro_flux.h
index 68b590e9b6..5b409dcdf2 100644
--- a/src/hydro/Gizmo/MFM/hydro_flux.h
+++ b/src/hydro/Gizmo/MFM/hydro_flux.h
@@ -150,8 +150,7 @@ hydro_gizmo_mfv_density_drift_term(const float mass_flux, const float dt,
 /**
  * @brief Add the gravitational contribution to the fluid velocity drift.
  *
- * This just needs to reset the particle velocity, which was incorrectly
- * drifted, since this is Gizmo MFM.
+ * (MFV only)
  *
  * @param v (drifted) particle velocity.
  * @param fluid_v Fluid velocity.
@@ -163,12 +162,7 @@ __attribute__((always_inline)) INLINE static void
 hydro_gizmo_mfv_extra_velocity_drift(float* v, const float* restrict fluid_v,
                                      const float* restrict v_full,
                                      const float* restrict a_grav,
-                                     float dt_kick_grav) {
-  /* Just reset particle velocity */
-  v[0] = v_full[0];
-  v[1] = v_full[1];
-  v[2] = v_full[2];
-}
+                                     float dt_kick_grav) {}
 
 /**
  * @brief Get the term required to update the MFV energy due to the change in
diff --git a/src/hydro/Gizmo/MFV/hydro_flux.h b/src/hydro/Gizmo/MFV/hydro_flux.h
index bab922bd10..aa9fc37835 100644
--- a/src/hydro/Gizmo/MFV/hydro_flux.h
+++ b/src/hydro/Gizmo/MFV/hydro_flux.h
@@ -190,12 +190,6 @@ hydro_gizmo_mfv_extra_velocity_drift(float* restrict v, float* restrict fluid_v,
                                      const float* restrict v_full,
                                      const float* restrict a_grav,
                                      float dt_kick_grav) {
-
-  /* Reset particle velocity */
-  v[0] = v_full[0];
-  v[1] = v_full[1];
-  v[2] = v_full[2];
-
   /* Drift fluid velocity */
   fluid_v[0] += a_grav[0] * dt_kick_grav;
   fluid_v[1] += a_grav[1] * dt_kick_grav;
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 698141905c..85a16a4a64 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -69,16 +69,13 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
   float W[5];
   hydro_part_get_primitive_variables(p, W);
 
-  /* v_full is the actual velocity of the particle, v is its
-     hydrodynamical velocity. The time step depends on the relative difference
-     of the two. */
-  float vrel[3];
-  vrel[0] = W[1] - xp->v_full[0];
-  vrel[1] = W[2] - xp->v_full[1];
-  vrel[2] = W[3] - xp->v_full[2];
+  /* The time step depends on the relative difference of the fluid velocity and
+   * the particle velocity. */
+  float v_rel[3];
+  hydro_part_get_relative_fluid_velocity(p, v_rel);
   const float rhoinv = (W[0] > 0.0f) ? 1.0f / W[0] : 0.0f;
   float vmax =
-      sqrtf(vrel[0] * vrel[0] + vrel[1] * vrel[1] + vrel[2] * vrel[2]) +
+      sqrtf(v_rel[0] * v_rel[0] + v_rel[1] * v_rel[1] + v_rel[2] * v_rel[2]) +
       sqrtf(hydro_gamma * W[4] * rhoinv);
   vmax = max(vmax, p->timestepvars.vmax);
 
@@ -586,28 +583,38 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
   }
 #endif
 
+  /* Reset the particle velocity. (undo the drift) */
+  hydro_set_particle_velocity(p, xp->v_full);
+
   float W[5];
   hydro_part_get_primitive_variables(p, W);
+
+  /* Use the fluid velocity in the rest frame of the particle for the time
+   * extrapolation to preserve Galilean invariance. */
+  float v_rel[3];
+  hydro_part_get_relative_fluid_velocity(p, v_rel);
+
   float gradrho[3], gradvx[3], gradvy[3], gradvz[3], gradP[3];
   hydro_part_get_gradients(p, gradrho, gradvx, gradvy, gradvz, gradP);
 
   const float divv = gradvx[0] + gradvy[1] + gradvz[2];
 
   float Wprime[5];
-  Wprime[0] = W[0] - dt_therm * (W[0] * divv + W[1] * gradrho[0] +
-                                 W[2] * gradrho[1] + W[3] * gradrho[2]);
+  Wprime[0] = W[0] - dt_therm * (W[0] * divv + v_rel[0] * gradrho[0] +
+                                 v_rel[1] * gradrho[1] + v_rel[2] * gradrho[2]);
   if (W[0] != 0.0f) {
     const float rhoinv = 1.0f / W[0];
-    Wprime[1] = W[1] - dt_therm * (W[1] * divv + rhoinv * gradP[0]);
-    Wprime[2] = W[2] - dt_therm * (W[2] * divv + rhoinv * gradP[1]);
-    Wprime[3] = W[3] - dt_therm * (W[3] * divv + rhoinv * gradP[2]);
+    Wprime[1] = W[1] - dt_therm * (v_rel[0] * divv + rhoinv * gradP[0]);
+    Wprime[2] = W[2] - dt_therm * (v_rel[1] * divv + rhoinv * gradP[1]);
+    Wprime[3] = W[3] - dt_therm * (v_rel[2] * divv + rhoinv * gradP[2]);
   } else {
     Wprime[1] = 0.0f;
     Wprime[2] = 0.0f;
     Wprime[3] = 0.0f;
   }
-  Wprime[4] = W[4] - dt_therm * (hydro_gamma * W[4] * divv + W[1] * gradP[0] +
-                                 W[2] * gradP[1] + W[3] * gradP[2]);
+  Wprime[4] =
+      W[4] - dt_therm * (hydro_gamma * W[4] * divv + v_rel[0] * gradP[0] +
+                         v_rel[1] * gradP[1] + v_rel[2] * gradP[2]);
 
   W[0] = Wprime[0];
   W[1] = Wprime[1];
diff --git a/src/hydro/Gizmo/hydro_getters.h b/src/hydro/Gizmo/hydro_getters.h
index 7a83e09368..775e548e56 100644
--- a/src/hydro/Gizmo/hydro_getters.h
+++ b/src/hydro/Gizmo/hydro_getters.h
@@ -313,6 +313,20 @@ __attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
   }
 }
 
+/**
+ * @brief Compute the fluid velocity in the reference frame co-moving with the
+ * particle.
+ *
+ * @param p The #part
+ * @param v_rel (return) The relative fluid velocity.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_part_get_relative_fluid_velocity(const struct part* p, float* v_rel) {
+  v_rel[0] = p->fluid_v[0] - p->v[0];
+  v_rel[1] = p->fluid_v[1] - p->v[1];
+  v_rel[2] = p->fluid_v[2] - p->v[2];
+}
+
 /**
  * @brief Returns the time derivative of co-moving internal energy of a particle
  *
@@ -326,6 +340,9 @@ hydro_get_comoving_internal_energy_dt(const struct part* restrict p) {
   float W[5];
   hydro_part_get_primitive_variables(p, W);
 
+  float v_rel[3];
+  hydro_part_get_relative_fluid_velocity(p, v_rel);
+
   if (W[0] <= 0.0f) {
     return 0.0f;
   }
@@ -343,8 +360,9 @@ hydro_get_comoving_internal_energy_dt(const struct part* restrict p) {
                (gradP[i] - rho_inv * W[4] * gradrho[i]);
   }
 
-  const float du_dt = -(W[1] * gradu[0] + W[2] * gradu[1] + W[3] * gradu[2]) -
-                      rho_inv * W[4] * divv;
+  const float du_dt =
+      -(v_rel[0] * gradu[0] + v_rel[1] * gradu[1] + v_rel[2] * gradu[2]) -
+      rho_inv * W[4] * divv;
 
   return du_dt;
 }
diff --git a/src/hydro/Gizmo/hydro_setters.h b/src/hydro/Gizmo/hydro_setters.h
index 2715f7baff..d12532845c 100644
--- a/src/hydro/Gizmo/hydro_setters.h
+++ b/src/hydro/Gizmo/hydro_setters.h
@@ -415,4 +415,20 @@ hydro_set_init_internal_energy(struct part* p, float u_init) {
   p->conserved.energy = u_init;
 }
 
+/**
+ * @brief Set the *particle* velocity of a particle.
+ *
+ * This must be the velocity at which the particle is drifted during its
+ * timestep (i.e. xp.v_full).
+ *
+ * @param p The #part to write to.
+ * @param v The new particle velocity.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_particle_velocity(struct part* p, float* v) {
+  p->v[0] = v[0];
+  p->v[1] = v[1];
+  p->v[2] = v[2];
+}
+
 #endif /* SWIFT_GIZMO_HYDRO_SETTERS_H */
-- 
GitLab