Pass an extra argument to the function returning the particle's softening....

Pass an extra argument to the function returning the particle's softening. Removed the field from the gpart themselves.

src/gravity.h

@@ -24,16 +24,17 @@
 
 /* Local headers. */
 #include "const.h"
-#include "engine.h"
 #include "inline.h"
 #include "part.h"
-#include "space.h"
 
 /* So far only one model here */
 /* Straight-forward import */
 #include "./gravity/Default/gravity.h"
 #include "./gravity/Default/gravity_iact.h"
 
+struct engine;
+struct space;
+
 void gravity_exact_force_ewald_init(double boxSize);
 void gravity_exact_force_ewald_free();
 void gravity_exact_force_compute(struct space *s, const struct engine *e);

src/gravity/Default/gravity.h

@@ -21,6 +21,7 @@
 #define SWIFT_DEFAULT_GRAVITY_H
 
 #include <float.h>
+
 #include "cosmology.h"
 #include "gravity_properties.h"
 #include "minmax.h"
@@ -42,9 +43,9 @@ __attribute__((always_inline)) INLINE static float gravity_get_mass(
  * @param gp The particle of interest
  */
 __attribute__((always_inline)) INLINE static float gravity_get_softening(
-    const struct gpart* restrict gp) {
+    const struct gpart* gp, const struct gravity_props* restrict grav_props) {
 
-  return gp->epsilon;
+  return grav_props->epsilon;
 }
 
 /**
@@ -102,8 +103,7 @@ gravity_compute_timestep_self(const struct gpart* const gp,
 
   const float ac_inv = (ac2 > 0.f) ? 1.f / sqrtf(ac2) : FLT_MAX;
 
-  // MATTHIEU cosmological evolution of the softening?
-  const float epsilon = gravity_get_softening(gp);
+  const float epsilon = gravity_get_softening(gp, grav_props);
 
   /* Note that 0.66666667 = 2. (from Gadget) / 3. (Plummer softening) */
   const float dt =
@@ -183,7 +183,6 @@ __attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
     struct gpart* gp, const struct gravity_props* grav_props) {
 
   gp->time_bin = 0;
-  gp->epsilon = grav_props->epsilon;
 
   gravity_init_gpart(gp);
 }

src/gravity/Default/gravity_debug.h

@@ -22,9 +22,9 @@
 __attribute__((always_inline)) INLINE static void gravity_debug_particle(
     const struct gpart* p) {
   printf(
-      "mass=%.3e epsilon=%.5e time_bin=%d\n"
+      "mass=%.3e time_bin=%d\n"
       "x=[%.5e,%.5e,%.5e], v_full=[%.5e,%.5e,%.5e], a=[%.5e,%.5e,%.5e]\n",
-      p->mass, p->epsilon, p->time_bin, p->x[0], p->x[1], p->x[2], p->v_full[0],
+      p->mass, p->time_bin, p->x[0], p->x[1], p->x[2], p->v_full[0],
       p->v_full[1], p->v_full[2], p->a_grav[0], p->a_grav[1], p->a_grav[2]);
 #ifdef SWIFT_DEBUG_CHECKS
   printf("num_interacted=%lld ti_drift=%lld ti_kick=%lld\n", p->num_interacted,

src/gravity/Default/gravity_part.h

@@ -22,35 +22,32 @@
 /* Gravity particle. */
 struct gpart {
 
-  /* Particle ID. If negative, it is the negative offset of the #part with
+  /*! Particle ID. If negative, it is the negative offset of the #part with
      which this gpart is linked. */
   long long id_or_neg_offset;
 
-  /* Particle position. */
+  /*! Particle position. */
   double x[3];
 
-  /* Offset between current position and position at last tree rebuild. */
+  /*! Offset between current position and position at last tree rebuild. */
   float x_diff[3];
 
-  /* Particle velocity. */
+  /*! Particle velocity. */
   float v_full[3];
 
-  /* Particle acceleration. */
+  /*! Particle acceleration. */
   float a_grav[3];
 
-  /* Particle mass. */
+  /*! Particle mass. */
   float mass;
 
-  /* Gravitational potential */
+  /*! Gravitational potential */
   float potential;
 
-  /* Softening length */
-  float epsilon;
-
-  /* Time-step length */
+  /*! Time-step length */
   timebin_t time_bin;
 
-  /* Type of the #gpart (DM, gas, star, ...) */
+  /*! Type of the #gpart (DM, gas, star, ...) */
   enum part_type type;
 
 #ifdef SWIFT_DEBUG_CHECKS

src/gravity_cache.h

@@ -156,7 +156,7 @@ __attribute__((always_inline)) INLINE static void gravity_cache_populate(
     timebin_t max_active_bin, struct gravity_cache *c,
     const struct gpart *restrict gparts, int gcount, int gcount_padded,
     const double shift[3], const float CoM[3], float r_max2, float theta_crit2,
-    const struct cell *cell) {
+    const struct cell *cell, const struct gravity_props *grav_props) {
 
   /* Make the compiler understand we are in happy vectorization land */
   swift_declare_aligned_ptr(float, x, c->x, SWIFT_CACHE_ALIGNMENT);
@@ -174,7 +174,7 @@ __attribute__((always_inline)) INLINE static void gravity_cache_populate(
     x[i] = (float)(gparts[i].x[0] - shift[0]);
     y[i] = (float)(gparts[i].x[1] - shift[1]);
     z[i] = (float)(gparts[i].x[2] - shift[2]);
-    epsilon[i] = gparts[i].epsilon;
+    epsilon[i] = gravity_get_softening(&gparts[i], grav_props);
     m[i] = gparts[i].mass;
     active[i] = (int)(gparts[i].time_bin <= max_active_bin);
 
@@ -226,7 +226,8 @@ gravity_cache_populate_no_mpole(timebin_t max_active_bin,
                                 struct gravity_cache *c,
                                 const struct gpart *restrict gparts, int gcount,
                                 int gcount_padded, const double shift[3],
-                                const struct cell *cell) {
+                                const struct cell *cell,
+                                const struct gravity_props *grav_props) {
 
   /* Make the compiler understand we are in happy vectorization land */
   swift_declare_aligned_ptr(float, x, c->x, SWIFT_CACHE_ALIGNMENT);
@@ -242,7 +243,7 @@ gravity_cache_populate_no_mpole(timebin_t max_active_bin,
     x[i] = (float)(gparts[i].x[0] - shift[0]);
     y[i] = (float)(gparts[i].x[1] - shift[1]);
     z[i] = (float)(gparts[i].x[2] - shift[2]);
-    epsilon[i] = gparts[i].epsilon;
+    epsilon[i] = gravity_get_softening(&gparts[i], grav_props);
     m[i] = gparts[i].mass;
     active[i] = (int)(gparts[i].time_bin <= max_active_bin);
   }

src/logger.c

@@ -218,12 +218,6 @@ void logger_log_gpart(struct gpart *p, unsigned int mask, size_t *offset,
     buff += 3 * sizeof(float);
   }
 
-  /* Particle smoothing length as a single float. */
-  if (mask & logger_mask_h) {
-    memcpy(buff, &p->epsilon, sizeof(float));
-    buff += sizeof(float);
-  }
-
   /* Particle constants, which is a bit more complicated. */
   if (mask & logger_mask_rho) {
     memcpy(buff, &p->mass, sizeof(float));
@@ -385,12 +379,6 @@ int logger_read_gpart(struct gpart *p, size_t *offset, const char *buff) {
     buff += 3 * sizeof(float);
   }
 
-  /* Particle smoothing length as a single float. */
-  if (mask & logger_mask_h) {
-    memcpy(&p->epsilon, buff, sizeof(float));
-    buff += sizeof(float);
-  }
-
   /* Particle constants, which is a bit more complicated. */
   if (mask & logger_mask_rho) {
     memcpy(&p->mass, buff, sizeof(float));

src/partition.c

@@ -48,6 +48,7 @@
 
 /* Local headers. */
 #include "debug.h"
+#include "engine.h"
 #include "error.h"
 #include "partition.h"
 #include "restart.h"

src/runner_doiact_grav.h

@@ -533,10 +533,10 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
   /* Fill the caches */
   gravity_cache_populate(e->max_active_bin, ci_cache, ci->gparts, gcount_i,
                          gcount_padded_i, shift_i, CoM_j, rmax2_j, theta_crit2,
-                         ci);
+                         ci, e->gravity_properties);
   gravity_cache_populate(e->max_active_bin, cj_cache, cj->gparts, gcount_j,
                          gcount_padded_j, shift_j, CoM_i, rmax2_i, theta_crit2,
-                         cj);
+                         cj, e->gravity_properties);
 
   /* Can we use the Newtonian version or do we need the truncated one ? */
   if (!periodic) {
@@ -675,7 +675,7 @@ static INLINE void runner_doself_grav_pp_full(struct runner *r,
   const int gcount_padded = gcount - (gcount % VEC_SIZE) + VEC_SIZE;
 
   gravity_cache_populate_no_mpole(e->max_active_bin, ci_cache, gparts, gcount,
-                                  gcount_padded, loc, c);
+                                  gcount_padded, loc, c, e->gravity_properties);
 
   /* Ok... Here we go ! */
 
@@ -805,7 +805,7 @@ static INLINE void runner_doself_grav_pp_truncated(struct runner *r,
   const int gcount_padded = gcount - (gcount % VEC_SIZE) + VEC_SIZE;
 
   gravity_cache_populate_no_mpole(e->max_active_bin, ci_cache, gparts, gcount,
-                                  gcount_padded, loc, c);
+                                  gcount_padded, loc, c, e->gravity_properties);
 
   /* Ok... Here we go ! */