Merge branch 'correct_velocities_in_stats' into 'master'

Make the hydro schemes compute the dirfted velocities themselves when collecting statistics.

See merge request !474

src/hydro/Default/hydro.h

@@ -97,6 +97,23 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
   return p->mass;
 }
 
+/**
+ * @brief Returns the velocities drifted to the current time of a particle.
+ *
+ * @param p The particle of interest
+ * @param xp The extended data of the particle.
+ * @param dt The time since the last kick.
+ * @param v (return) The velocities at the current time.
+ */
+__attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
+    const struct part *restrict p, const struct xpart *xp, float dt,
+    float v[3]) {
+
+  v[0] = xp->v_full[0] + p->a_hydro[0] * dt;
+  v[1] = xp->v_full[1] + p->a_hydro[1] * dt;
+  v[2] = xp->v_full[2] + p->a_hydro[2] * dt;
+}
+
 /**
  * @brief Returns the time derivative of internal energy of a particle
  *

src/hydro/Gadget2/hydro.h

@@ -106,6 +106,23 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
   return p->mass;
 }
 
+/**
+ * @brief Returns the velocities drifted to the current time of a particle.
+ *
+ * @param p The particle of interest
+ * @param xp The extended data of the particle.
+ * @param dt The time since the last kick.
+ * @param v (return) The velocities at the current time.
+ */
+__attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
+    const struct part *restrict p, const struct xpart *xp, float dt,
+    float v[3]) {
+
+  v[0] = xp->v_full[0] + p->a_hydro[0] * dt;
+  v[1] = xp->v_full[1] + p->a_hydro[1] * dt;
+  v[2] = xp->v_full[2] + p->a_hydro[2] * dt;
+}
+
 /**
  * @brief Returns the time derivative of internal energy of a particle
  *

src/hydro/Gizmo/hydro.h

@@ -754,6 +754,32 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
   return p->conserved.mass;
 }
 
+/**
+ * @brief Returns the velocities drifted to the current time of a particle.
+ *
+ * @param p The particle of interest
+ * @param xp The extended data of the particle.
+ * @param dt The time since the last kick.
+ * @param v (return) The velocities at the current time.
+ */
+__attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
+    const struct part* restrict p, const struct xpart* xp, float dt,
+    float v[3]) {
+
+  if (p->conserved.mass > 0.) {
+    v[0] = p->primitives.v[0] +
+           p->conserved.flux.momentum[0] * dt / p->conserved.mass;
+    v[1] = p->primitives.v[1] +
+           p->conserved.flux.momentum[1] * dt / p->conserved.mass;
+    v[2] = p->primitives.v[2] +
+           p->conserved.flux.momentum[2] * dt / p->conserved.mass;
+  } else {
+    v[0] = p->primitives.v[0];
+    v[1] = p->primitives.v[1];
+    v[2] = p->primitives.v[2];
+  }
+}
+
 /**
  * @brief Returns the density of a particle
  *

src/hydro/Minimal/hydro.h

@@ -116,6 +116,23 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
   return p->mass;
 }
 
+/**
+ * @brief Returns the velocities drifted to the current time of a particle.
+ *
+ * @param p The particle of interest
+ * @param xp The extended data of the particle.
+ * @param dt The time since the last kick.
+ * @param v (return) The velocities at the current time.
+ */
+__attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
+    const struct part *restrict p, const struct xpart *xp, float dt,
+    float v[3]) {
+
+  v[0] = xp->v_full[0] + p->a_hydro[0] * dt;
+  v[1] = xp->v_full[1] + p->a_hydro[1] * dt;
+  v[2] = xp->v_full[2] + p->a_hydro[2] * dt;
+}
+
 /**
  * @brief Returns the time derivative of internal energy of a particle
  *

src/hydro/PressureEntropy/hydro.h

@@ -106,6 +106,23 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
   return p->mass;
 }
 
+/**
+ * @brief Returns the velocities drifted to the current time of a particle.
+ *
+ * @param p The particle of interest
+ * @param xp The extended data of the particle.
+ * @param dt The time since the last kick.
+ * @param v (return) The velocities at the current time.
+ */
+__attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
+    const struct part *restrict p, const struct xpart *xp, float dt,
+    float v[3]) {
+
+  v[0] = xp->v_full[0] + p->a_hydro[0] * dt;
+  v[1] = xp->v_full[1] + p->a_hydro[1] * dt;
+  v[2] = xp->v_full[2] + p->a_hydro[2] * dt;
+}
+
 /**
  * @brief Returns the time derivative of internal energy of a particle
  *

src/hydro/Shadowswift/hydro.h

@@ -544,6 +544,18 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
   return p->conserved.mass;
 }
 
+/**
+ * @brief Returns the velocities drifted to the current time of a particle.
+ *
+ * @param p The particle of interest
+ * @param xp The extended data of the particle.
+ * @param dt The time since the last kick.
+ * @param v (return) The velocities at the current time.
+ */
+__attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
+    const struct part* restrict p, const struct xpart* xp, float dt,
+    float v[3]) {}
+
 /**
  * @brief Returns the density of a particle
  *

src/statistics.c

@@ -135,22 +135,11 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
     const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
     const float dt = (ti_current - ((ti_begin + ti_end) / 2)) * timeBase;
 
-    /* Get the total acceleration */
-    float a_tot[3] = {p->a_hydro[0], p->a_hydro[1], p->a_hydro[2]};
-    if (gp != NULL) {
-      a_tot[0] += gp->a_grav[0];
-      a_tot[1] += gp->a_grav[1];
-      a_tot[2] += gp->a_grav[2];
-    }
-
-    /* Extrapolate velocities to current time */
-    const float v[3] = {xp->v_full[0] + a_tot[0] * dt,
-                        xp->v_full[1] + a_tot[1] * dt,
-                        xp->v_full[2] + a_tot[2] * dt};
-
-    const float m = hydro_get_mass(p);
+    float v[3];
+    hydro_get_drifted_velocities(p, xp, dt, v);
     const double x[3] = {p->x[0], p->x[1], p->x[2]};
-    const float rho = hydro_get_density(p);
+    const float m = hydro_get_mass(p);
+    const float entropy = hydro_get_entropy(p);
     const float u_int = hydro_get_internal_energy(p);
 
     /* Collect mass */
@@ -182,7 +171,7 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
     }
 
     /* Collect entropy */
-    stats.entropy += m * gas_entropy_from_internal_energy(rho, u_int);
+    stats.entropy += m * entropy;
   }
 
   /* Now write back to memory */