Common interface to access the mass and the density of particles.

src/common_io.c

@@ -42,6 +42,7 @@
 /* Local includes. */
 #include "const.h"
 #include "error.h"
+#include "hydro.h"
 #include "kernel_hydro.h"
 #include "part.h"
 #include "units.h"
@@ -614,7 +615,7 @@ void duplicate_hydro_gparts(struct part* const parts,
     gparts[i + Ndm].v_full[1] = parts[i].v[1];
     gparts[i + Ndm].v_full[2] = parts[i].v[2];
 
-    gparts[i + Ndm].mass = parts[i].mass;
+    gparts[i + Ndm].mass = hydro_get_mass(&parts[i]);
 
     /* Link the particles */
     gparts[i + Ndm].id_or_neg_offset = -i;

src/drift.h

@@ -65,8 +65,6 @@ __attribute__((always_inline)) INLINE static void drift_gpart(
 __attribute__((always_inline)) INLINE static void drift_part(
     struct part *restrict p, struct xpart *restrict xp, float dt,
     double timeBase, int ti_old, int ti_current) {
-  /* Useful quantity */
-  const float h_inv = 1.0f / p->h;
 
   /* Drift... */
   p->x[0] += xp->v_full[0] * dt;
@@ -78,22 +76,8 @@ __attribute__((always_inline)) INLINE static void drift_part(
   p->v[1] += p->a_hydro[1] * dt;
   p->v[2] += p->a_hydro[2] * dt;
 
-  /* Predict smoothing length */
-  const float w1 = p->force.h_dt * h_inv * dt;
-  if (fabsf(w1) < 0.2f)
-    p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
-  else
-    p->h *= expf(w1);
-
-  /* Predict density */
-  const float w2 = -hydro_dimension * w1;
-  if (fabsf(w2) < 0.2f)
-    p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
-  else
-    p->rho *= expf(w2);
-
   /* Predict the values of the extra fields */
-  hydro_predict_extra(p, xp, ti_old, ti_current, timeBase);
+  hydro_predict_extra(p, xp, dt, ti_old, ti_current, timeBase);
 
   /* Compute offset since last cell construction */
   xp->x_diff[0] -= xp->v_full[0] * dt;

src/hydro/Default/hydro.h

@@ -73,6 +73,28 @@ __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
   return p->force.soundspeed;
 }
 
+/**
+ * @brief Returns the density of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_density(
+    const struct part *restrict p) {
+
+  return p->rho;
+}
+
+/**
+ * @brief Returns the mass of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_mass(
+    const struct part *restrict p) {
+
+  return p->mass;
+}
+
 /**
  * @brief Modifies the thermal state of a particle to the imposed internal
  * energy
@@ -288,16 +310,31 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
  *
  * @param p The particle
  * @param xp The extended data of the particle
+ * @param dt The drift time-step.
  * @param t0 The time at the start of the drift
  * @param t1 The time at the end of the drift
  * @param timeBase The minimal time-step size
  */
 __attribute__((always_inline)) INLINE static void hydro_predict_extra(
-    struct part *restrict p, struct xpart *restrict xp, int t0, int t1,
-    double timeBase) {
+    struct part *restrict p, struct xpart *restrict xp, float dt, int t0,
+    int t1, double timeBase) {
   float u, w;
 
-  const float dt = (t1 - t0) * timeBase;
+  const float h_inv = 1.f / p->h;
+
+  /* Predict smoothing length */
+  const float w1 = p->force.h_dt * h_inv * dt;
+  if (fabsf(w1) < 0.2f)
+    p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
+  else
+    p->h *= expf(w1);
+
+  /* Predict density */
+  const float w2 = -hydro_dimension * w1;
+  if (fabsf(w2) < 0.2f)
+    p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
+  else
+    p->rho *= expf(w2);
 
   /* Predict internal energy */
   w = p->force.u_dt / p->u * dt;

src/hydro/Default/hydro_part.h

@@ -28,6 +28,8 @@ struct xpart {
   /* Velocity at the last full step. */
   float v_full[3];
 
+  float u_full;
+
   /* Old density. */
   float omega;
 

src/hydro/Gadget2/hydro.h

@@ -20,6 +20,7 @@
 #define SWIFT_GADGET2_HYDRO_H
 
 #include "adiabatic_index.h"
+#include "approx_math.h"
 #include "dimension.h"
 #include "equation_of_state.h"
 #include "hydro_properties.h"
@@ -77,6 +78,28 @@ __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
   return p->force.soundspeed;
 }
 
+/**
+ * @brief Returns the density of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_density(
+    const struct part *restrict p) {
+
+  return p->rho;
+}
+
+/**
+ * @brief Returns the mass of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_mass(
+    const struct part *restrict p) {
+
+  return p->mass;
+}
+
 /**
  * @brief Modifies the thermal state of a particle to the imposed internal
  * energy
@@ -285,13 +308,30 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
  *
  * @param p The particle
  * @param xp The extended data of the particle
+ * @param dt The drift time-step.
  * @param t0 The time at the start of the drift
  * @param t1 The time at the end of the drift
  * @param timeBase The minimal time-step size
  */
 __attribute__((always_inline)) INLINE static void hydro_predict_extra(
-    struct part *restrict p, const struct xpart *restrict xp, int t0, int t1,
-    double timeBase) {
+    struct part *restrict p, const struct xpart *restrict xp, float dt, int t0,
+    int t1, double timeBase) {
+
+  const float h_inv = 1.f / p->h;
+
+  /* Predict smoothing length */
+  const float w1 = p->force.h_dt * h_inv * dt;
+  if (fabsf(w1) < 0.2f)
+    p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
+  else
+    p->h *= expf(w1);
+
+  /* Predict density */
+  const float w2 = -hydro_dimension * w1;
+  if (fabsf(w2) < 0.2f)
+    p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
+  else
+    p->rho *= expf(w2);
 
   /* Drift the pressure */
   const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;

src/hydro/Gizmo/hydro.h

@@ -275,6 +275,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
  *
  * @param p Particle to act upon.
  * @param xp The extended particle data to act upon.
+ * @param dt The drift time-step.
  * @param t0 Integer start time of the drift interval.
  * @param t1 Integer end time of the drift interval.
  * @param timeBase Conversion factor between integer and physical time.
@@ -282,10 +283,16 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
 __attribute__((always_inline)) INLINE static void hydro_predict_extra(
     struct part* p, struct xpart* xp, int t0, int t1, double timeBase) {
 
-  // return;
-  float dt = (t1 - t0) * timeBase;
-  float h_inv = 1.0f / p->h;
-  float w = -hydro_dimension * p->force.h_dt * h_inv * dt;
+  const float dt = (t1 - t0) * timeBase;
+  const float h_inv = 1.0f / p->h;
+  const float w = -hydro_dimension * p->force.h_dt * h_inv * dt;
+
+  /* Predict smoothing length */
+  const float w1 = p->force.h_dt * h_inv * dt;
+  if (fabsf(w1) < 0.2f)
+    p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
+  else
+    p->h *= expf(w1);
 
   if (fabsf(w) < 0.2f) {
     p->primitives.rho *= approx_expf(w);
@@ -432,7 +439,6 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
  *
  * @param p The particle of interest.
  * @param dt Time since the last kick.
- * @return Internal energy of the particle.
  */
 __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
     const struct part* restrict p, float dt) {
@@ -445,7 +451,6 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
  *
  * @param p The particle of interest.
  * @param dt Time since the last kick.
- * @return Entropy of the particle.
  */
 __attribute__((always_inline)) INLINE static float hydro_get_entropy(
     const struct part* restrict p, float dt) {
@@ -458,7 +463,6 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
  *
  * @param p The particle of interest.
  * @param dt Time since the last kick.
- * @param Sound speed of the particle.
  */
 __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
     const struct part* restrict p, float dt) {
@@ -471,10 +475,31 @@ __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
  *
  * @param p The particle of interest
  * @param dt Time since the last kick
- * @param Pressure of the particle.
  */
 __attribute__((always_inline)) INLINE static float hydro_get_pressure(
     const struct part* restrict p, float dt) {
 
   return p->primitives.P;
 }
+
+/**
+ * @brief Returns the mass of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_mass(
+    const struct part* restrict p) {
+
+  return p->conserved.mass;
+}
+
+/**
+ * @brief Returns the density of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_density(
+    const struct part* restrict p) {
+
+  return p->primitives.rho;
+}

src/hydro/Gizmo/hydro_debug.h

@@ -52,8 +52,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       "div_v=%.3e, "
       "wcount_dh=%.3e, "
       "curl_v=[%.3e,%.3e,%.3e], "
-      "wcount=%.3e}, "
-      "mass=%.3e\n",
+      "wcount=%.3e}\n",
       p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->a_hydro[0],
       p->a_hydro[1], p->a_hydro[2], p->h, p->ti_begin, p->ti_end,
       p->primitives.v[0], p->primitives.v[1], p->primitives.v[2],
@@ -79,5 +78,5 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       p->geometry.matrix_E[2][1], p->geometry.matrix_E[2][2],
       p->timestepvars.vmax, p->density.div_v, p->density.wcount_dh,
       p->density.curl_v[0], p->density.curl_v[1], p->density.curl_v[2],
-      p->density.wcount, p->mass);
+      p->density.wcount);
 }

src/hydro/Gizmo/hydro_part.h

@@ -25,6 +25,7 @@ struct xpart {
 
   /* Velocity at the last full step. */
   float v_full[3];
+
 } __attribute__((aligned(xpart_align)));
 
 /* Data of a single particle. */

src/hydro/Minimal/hydro.h

@@ -35,6 +35,7 @@
  */
 
 #include "adiabatic_index.h"
+#include "approx_math.h"
 #include "dimension.h"
 #include "equation_of_state.h"
 #include "hydro_properties.h"
@@ -102,6 +103,28 @@ __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
   return gas_soundspeed_from_internal_energy(p->rho, u);
 }
 
+/**
+ * @brief Returns the density of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_density(
+    const struct part *restrict p) {
+
+  return p->rho;
+}
+
+/**
+ * @brief Returns the mass of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_mass(
+    const struct part *restrict p) {
+
+  return p->mass;
+}
+
 /**
  * @brief Modifies the thermal state of a particle to the imposed internal
  * energy
@@ -286,15 +309,32 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
  * Additional hydrodynamic quantites are drifted forward in time here. These
  * include thermal quantities (thermal energy or total energy or entropy, ...).
  *
- * @param p The particle
- * @param xp The extended data of the particle
- * @param t0 The time at the start of the drift (on the timeline)
- * @param t1 The time at the end of the drift (on the timeline)
- * @param timeBase The minimal time-step size
+ * @param p The particle.
+ * @param xp The extended data of the particle.
+ * @param dt The drift time-step.
+ * @param t0 The time at the start of the drift (on the timeline).
+ * @param t1 The time at the end of the drift (on the timeline).
+ * @param timeBase The minimal time-step size.
  */
 __attribute__((always_inline)) INLINE static void hydro_predict_extra(
-    struct part *restrict p, const struct xpart *restrict xp, int t0, int t1,
-    double timeBase) {
+    struct part *restrict p, const struct xpart *restrict xp, float dt, int t0,
+    int t1, double timeBase) {
+
+  const float h_inv = 1.f / p->h;
+
+  /* Predict smoothing length */
+  const float w1 = p->force.h_dt * h_inv * dt;
+  if (fabsf(w1) < 0.2f)
+    p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
+  else
+    p->h *= expf(w1);
+
+  /* Predict density */
+  const float w2 = -hydro_dimension * w1;
+  if (fabsf(w2) < 0.2f)
+    p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
+  else
+    p->rho *= expf(w2);
 
   /* Drift the pressure */
   const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;

src/runner.c

@@ -702,7 +702,8 @@ static void runner_do_drift(struct cell *c, struct engine *e) {
       const float v[3] = {xp->v_full[0] + p->a_hydro[0] * half_dt,
                           xp->v_full[1] + p->a_hydro[1] * half_dt,
                           xp->v_full[2] + p->a_hydro[2] * half_dt};
-      const float m = p->mass;
+
+      const float m = hydro_get_mass(p);
 
       /* Collect mass */
       mass += m;

src/tools.c

@@ -442,43 +442,6 @@ void pairs_single_grav(double *dim, long long int pid,
       aabs[2]);
 }
 
-/**
- * @brief Test the density function by dumping it for two random parts.
- *
- * @param N number of intervals in [0,1].
- */
-void density_dump(int N) {
-
-  int k;
-  float r2[4] = {0.0f, 0.0f, 0.0f, 0.0f}, hi[4], hj[4];
-  struct part /**pi[4],  *pj[4],*/ Pi[4], Pj[4];
-
-  /* Init the interaction parameters. */
-  for (k = 0; k < 4; k++) {
-    Pi[k].mass = 1.0f;
-    Pi[k].rho = 0.0f;
-    Pi[k].density.wcount = 0.0f;
-    Pi[k].id = k;
-    Pj[k].mass = 1.0f;
-    Pj[k].rho = 0.0f;
-    Pj[k].density.wcount = 0.0f;
-    Pj[k].id = k + 4;
-    hi[k] = 1.0;
-    hj[k] = 1.0;
-  }
-
-  for (k = 0; k <= N; k++) {
-    r2[3] = r2[2];
-    r2[2] = r2[1];
-    r2[1] = r2[0];
-    r2[0] = ((float)k) / N;
-    Pi[0].density.wcount = 0;
-    Pj[0].density.wcount = 0;
-    runner_iact_density(r2[0], NULL, hi[0], hj[0], &Pi[0], &Pj[0]);
-    printf(" %e %e %e", r2[0], Pi[0].density.wcount, Pj[0].density.wcount);
-  }
-}
-
 /**
  * @brief Compute the force on a single particle brute-force.
  */
@@ -520,7 +483,7 @@ void engine_single_density(double *dim, long long int pid,
   /* Dump the result. */
   hydro_end_density(&p, 0);
   message("part %lli (h=%e) has wcount=%e, rho=%e.", p.id, p.h,
-          p.density.wcount, p.rho);
+          p.density.wcount, hydro_get_density(&p));
   fflush(stdout);
 }