Merge branch 'gravity_multi_dt' into 'master'

Gravity multi dt



See merge request !333

README

@@ -27,6 +27,7 @@ Valid options are:
   -f    {int} Overwrite the CPU frequency (Hz) to be used for time measurements.
   -g          Run with an external gravitational potential.
   -G          Run with self-gravity.
+  -M          Reconstruct the multipoles every time-step.
   -n    {int} Execute a fixed number of time steps. When unset use the time_end parameter to stop. 
   -s          Run with hydrodynamics.
   -S          Run with stars.

configure.ac

@@ -215,6 +215,21 @@ elif test "$gravity_force_checks" != "no"; then
    AC_DEFINE_UNQUOTED([SWIFT_GRAVITY_FORCE_CHECKS], [$enableval] ,[Enable gravity brute-force checks])
 fi
 
+# Check if we want to zero the gravity forces for all particles below some ID.
+AC_ARG_ENABLE([no-gravity-below-id],
+   [AS_HELP_STRING([--enable-no-gravity-below-id],
+     [Zeros the gravitational acceleration of all particles with an ID smaller than @<:@N@:>@]
+   )],
+   [no_gravity_below_id="$enableval"],
+   [no_gravity_below_id="no"]
+)
+if test "$no_gravity_below_id" == "yes"; then
+   AC_MSG_ERROR(Need to specify the ID below which particles get zero forces when using --enable-no-gravity-below-id!)
+elif test "$no_gravity_below_id" != "no"; then
+   AC_DEFINE_UNQUOTED([SWIFT_NO_GRAVITY_BELOW_ID], [$enableval] ,[Particles with smaller ID than this will have zero gravity forces])
+fi
+
+
 # Define HAVE_POSIX_MEMALIGN if it works.
 AX_FUNC_POSIX_MEMALIGN
 
@@ -854,12 +869,14 @@ AC_MSG_RESULT([
    Adiabatic index    : $with_gamma
    Riemann solver     : $with_riemann
    Cooling function   : $with_cooling
-   External potential : $with_potential
-   Multipole order    : $with_multipole_order
-
-   Task debugging     : $enable_task_debugging
-   Debugging checks   : $enable_debugging_checks
-   Gravity checks     : $gravity_force_checks
+   
+   External potential  : $with_potential
+   Multipole order     : $with_multipole_order
+   No gravity below ID : $no_gravity_below_id
+
+   Task debugging    : $enable_task_debugging
+   Debugging checks  : $enable_debugging_checks
+   Gravity checks    : $gravity_force_checks
 ])
 
 # Make sure the latest git revision string gets included

examples/main.c

@@ -82,6 +82,8 @@ void print_help_message() {
          "Run with an external gravitational potential.");
   printf("  %2s %8s %s\n", "-F", "", "Run with feedback.");
   printf("  %2s %8s %s\n", "-G", "", "Run with self-gravity.");
+  printf("  %2s %8s %s\n", "-M", "",
+         "Reconstruct the multipoles every time-step.");
   printf("  %2s %8s %s\n", "-n", "{int}",
          "Execute a fixed number of time steps. When unset use the time_end "
          "parameter to stop.");
@@ -164,6 +166,7 @@ int main(int argc, char *argv[]) {
   int with_stars = 0;
   int with_fp_exceptions = 0;
   int with_drift_all = 0;
+  int with_mpole_reconstruction = 0;
   int verbose = 0;
   int nr_threads = 1;
   int with_verbose_timers = 0;
@@ -172,7 +175,8 @@ int main(int argc, char *argv[]) {
 
   /* Parse the parameters */
   int c;
-  while ((c = getopt(argc, argv, "acCdDef:FgGhn:sSt:Tv:y:")) != -1) switch (c) {
+  while ((c = getopt(argc, argv, "acCdDef:FgGhMn:sSt:Tv:y:")) != -1)
+    switch (c) {
       case 'a':
         with_aff = 1;
         break;
@@ -210,6 +214,9 @@ int main(int argc, char *argv[]) {
       case 'h':
         if (myrank == 0) print_help_message();
         return 0;
+      case 'M':
+        with_mpole_reconstruction = 1;
+        break;
       case 'n':
         if (sscanf(optarg, "%d", &nsteps) != 1) {
           if (myrank == 0) printf("Error parsing fixed number of steps.\n");
@@ -521,6 +528,8 @@ int main(int argc, char *argv[]) {
   /* Construct the engine policy */
   int engine_policies = ENGINE_POLICY | engine_policy_steal;
   if (with_drift_all) engine_policies |= engine_policy_drift_all;
+  if (with_mpole_reconstruction)
+    engine_policies |= engine_policy_reconstruct_mpoles;
   if (with_hydro) engine_policies |= engine_policy_hydro;
   if (with_self_gravity) engine_policies |= engine_policy_self_gravity;
   if (with_external_gravity) engine_policies |= engine_policy_external_gravity;

src/cell.c

@@ -1103,33 +1103,92 @@ void cell_reset_task_counters(struct cell *c) {
 }
 
 /**
- * @brief Checks whether the cells are direct neighbours ot not. Both cells have
- * to be of the same size
+ * @brief Recursively construct all the multipoles in a cell hierarchy.
  *
- * @param ci First #cell.
- * @param cj Second #cell.
- *
- * @todo Deal with periodicity.
+ * @param c The #cell.
  */
-int cell_are_neighbours(const struct cell *restrict ci,
-                        const struct cell *restrict cj) {
+void cell_make_multipoles(struct cell *c, integertime_t ti_current) {
 
-#ifdef SWIFT_DEBUG_CHECKS
-  if (ci->width[0] != cj->width[0]) error("Cells of different size !");
-#endif
+  /* Reset everything */
+  gravity_reset(c->multipole);
+
+  if (c->split) {
+
+    /* Compute CoM of all progenies */
+    double CoM[3] = {0., 0., 0.};
+    double mass = 0.;
+
+    for (int k = 0; k < 8; ++k) {
+      if (c->progeny[k] != NULL) {
+        const struct gravity_tensors *m = c->progeny[k]->multipole;
+        CoM[0] += m->CoM[0] * m->m_pole.M_000;
+        CoM[1] += m->CoM[1] * m->m_pole.M_000;
+        CoM[2] += m->CoM[2] * m->m_pole.M_000;
+        mass += m->m_pole.M_000;
+      }
+    }
+    c->multipole->CoM[0] = CoM[0] / mass;
+    c->multipole->CoM[1] = CoM[1] / mass;
+    c->multipole->CoM[2] = CoM[2] / mass;
+
+    /* Now shift progeny multipoles and add them up */
+    struct multipole temp;
+    double r_max = 0.;
+    for (int k = 0; k < 8; ++k) {
+      if (c->progeny[k] != NULL) {
+        const struct cell *cp = c->progeny[k];
+        const struct multipole *m = &cp->multipole->m_pole;
+
+        /* Contribution to multipole */
+        gravity_M2M(&temp, m, c->multipole->CoM, cp->multipole->CoM);
+        gravity_multipole_add(&c->multipole->m_pole, &temp);
+
+        /* Upper limit of max CoM<->gpart distance */
+        const double dx = c->multipole->CoM[0] - cp->multipole->CoM[0];
+        const double dy = c->multipole->CoM[1] - cp->multipole->CoM[1];
+        const double dz = c->multipole->CoM[2] - cp->multipole->CoM[2];
+        const double r2 = dx * dx + dy * dy + dz * dz;
+        r_max = max(r_max, cp->multipole->r_max + sqrt(r2));
+      }
+    }
+    /* Alternative upper limit of max CoM<->gpart distance */
+    const double dx = c->multipole->CoM[0] > c->loc[0] + c->width[0] / 2.
+                          ? c->multipole->CoM[0] - c->loc[0]
+                          : c->loc[0] + c->width[0] - c->multipole->CoM[0];
+    const double dy = c->multipole->CoM[1] > c->loc[1] + c->width[1] / 2.
+                          ? c->multipole->CoM[1] - c->loc[1]
+                          : c->loc[1] + c->width[1] - c->multipole->CoM[1];
+    const double dz = c->multipole->CoM[2] > c->loc[2] + c->width[2] / 2.
+                          ? c->multipole->CoM[2] - c->loc[2]
+                          : c->loc[2] + c->width[2] - c->multipole->CoM[2];
+
+    /* Take minimum of both limits */
+    c->multipole->r_max = min(r_max, sqrt(dx * dx + dy * dy + dz * dz));
 
-  /* Maximum allowed distance */
-  const double min_dist =
-      1.2 * ci->width[0]; /* 1.2 accounts for rounding errors */
+  } else {
 
-  /* (Manhattan) Distance between the cells */
-  for (int k = 0; k < 3; k++) {
-    const double center_i = ci->loc[k];
-    const double center_j = cj->loc[k];
-    if (fabs(center_i - center_j) > min_dist) return 0;
+    if (c->gcount > 0) {
+      gravity_P2M(c->multipole, c->gparts, c->gcount);
+      const double dx = c->multipole->CoM[0] > c->loc[0] + c->width[0] / 2.
+                            ? c->multipole->CoM[0] - c->loc[0]
+                            : c->loc[0] + c->width[0] - c->multipole->CoM[0];
+      const double dy = c->multipole->CoM[1] > c->loc[1] + c->width[1] / 2.
+                            ? c->multipole->CoM[1] - c->loc[1]
+                            : c->loc[1] + c->width[1] - c->multipole->CoM[1];
+      const double dz = c->multipole->CoM[2] > c->loc[2] + c->width[2] / 2.
+                            ? c->multipole->CoM[2] - c->loc[2]
+                            : c->loc[2] + c->width[2] - c->multipole->CoM[2];
+      c->multipole->r_max = sqrt(dx * dx + dy * dy + dz * dz);
+    } else {
+      gravity_multipole_init(&c->multipole->m_pole);
+      c->multipole->CoM[0] = c->loc[0] + c->width[0] / 2.;
+      c->multipole->CoM[1] = c->loc[1] + c->width[1] / 2.;
+      c->multipole->CoM[2] = c->loc[2] + c->width[2] / 2.;
+      c->multipole->r_max = 0.;
+    }
   }
 
-  return 1;
+  c->ti_old_multipole = ti_current;
 }
 
 /**
@@ -1145,6 +1204,8 @@ void cell_check_multipole(struct cell *c, void *data) {
   struct gravity_tensors ma;
   const double tolerance = 1e-3; /* Relative */
 
+  return;
+
   /* First recurse */
   if (c->split)
     for (int k = 0; k < 8; k++)

src/cell.h

@@ -351,8 +351,7 @@ int cell_link_gparts(struct cell *c, struct gpart *gparts);
 int cell_link_sparts(struct cell *c, struct spart *sparts);
 void cell_convert_hydro(struct cell *c, void *data);
 void cell_clean_links(struct cell *c, void *data);
-int cell_are_neighbours(const struct cell *restrict ci,
-                        const struct cell *restrict cj);
+void cell_make_multipoles(struct cell *c, integertime_t ti_current);
 void cell_check_multipole(struct cell *c, void *data);
 void cell_clean(struct cell *c);
 void cell_check_particle_drift_point(struct cell *c, void *data);

src/engine.c

@@ -1693,7 +1693,7 @@ void engine_make_self_gravity_tasks(struct engine *e) {
 
       /* Are the cells to close for a MM interaction ? */
       if (!gravity_multipole_accept(ci->multipole, cj->multipole,
-                                    theta_crit_inv))
+                                    theta_crit_inv, 1))
         scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0, ci,
                           cj, 1);
     }
@@ -3225,6 +3225,15 @@ void engine_step(struct engine *e) {
   /* Are we drifting everything (a la Gadget/GIZMO) ? */
   if (e->policy & engine_policy_drift_all) engine_drift_all(e);
 
+  /* Are we reconstructing the multipoles or drifting them ?*/
+  if (e->policy & engine_policy_self_gravity) {
+
+    if (e->policy & engine_policy_reconstruct_mpoles)
+      engine_reconstruct_multipoles(e);
+    else
+      engine_drift_top_multipoles(e);
+  }
+
   /* Print the number of active tasks ? */
   if (e->verbose) engine_print_task_counts(e);
 
@@ -3248,9 +3257,6 @@ void engine_step(struct engine *e) {
     gravity_exact_force_compute(e->s, e);
 #endif
 
-  /* Do we need to drift the top-level multipoles ? */
-  if (e->policy & engine_policy_self_gravity) engine_drift_top_multipoles(e);
-
   /* Start all the tasks. */
   TIMER_TIC;
   engine_launch(e, e->nr_threads);
@@ -3447,6 +3453,39 @@ void engine_drift_top_multipoles(struct engine *e) {
             clocks_getunit());
 }
 
+void engine_do_reconstruct_multipoles_mapper(void *map_data, int num_elements,
+                                             void *extra_data) {
+
+  struct engine *e = (struct engine *)extra_data;
+  struct cell *cells = (struct cell *)map_data;
+
+  for (int ind = 0; ind < num_elements; ind++) {
+    struct cell *c = &cells[ind];
+    if (c != NULL && c->nodeID == e->nodeID) {
+
+      /* Construct the multipoles in this cell hierarchy */
+      cell_make_multipoles(c, e->ti_current);
+    }
+  }
+}
+
+/**
+ * @brief Reconstruct all the multipoles at all the levels in the tree.
+ *
+ * @param e The #engine.
+ */
+void engine_reconstruct_multipoles(struct engine *e) {
+
+  const ticks tic = getticks();
+
+  threadpool_map(&e->threadpool, engine_do_reconstruct_multipoles_mapper,
+                 e->s->cells_top, e->s->nr_cells, sizeof(struct cell), 10, e);
+
+  if (e->verbose)
+    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+            clocks_getunit());
+}
+
 /**
  * @brief Create and fill the proxies.
  *

src/engine.h

@@ -66,9 +66,10 @@ enum engine_policy {
   engine_policy_external_gravity = (1 << 9),
   engine_policy_cosmology = (1 << 10),
   engine_policy_drift_all = (1 << 11),
-  engine_policy_cooling = (1 << 12),
-  engine_policy_sourceterms = (1 << 13),
-  engine_policy_stars = (1 << 14)
+  engine_policy_reconstruct_mpoles = (1 << 12),
+  engine_policy_cooling = (1 << 13),
+  engine_policy_sourceterms = (1 << 14),
+  engine_policy_stars = (1 << 15)
 };
 
 extern const char *engine_policy_names[];
@@ -256,6 +257,7 @@ void engine_compute_next_snapshot_time(struct engine *e);
 void engine_unskip(struct engine *e);
 void engine_drift_all(struct engine *e);
 void engine_drift_top_multipoles(struct engine *e);
+void engine_reconstruct_multipoles(struct engine *e);
 void engine_dump_snapshot(struct engine *e);
 void engine_init(struct engine *e, struct space *s,
                  const struct swift_params *params, int nr_nodes, int nodeID,

src/gravity.c

@@ -32,6 +32,13 @@
 #include "error.h"
 #include "version.h"
 
+struct exact_force_data {
+  const struct engine *e;
+  const struct space *s;
+  int counter_global;
+  double const_G;
+};
+
 /**
  * @brief Checks whether the file containing the exact accelerations for
  * the current choice of parameters already exists.
@@ -83,32 +90,23 @@ int gravity_exact_force_file_exits(const struct engine *e) {
 }
 
 /**
- * @brief Run a brute-force gravity calculation for a subset of particles.
- *
- * All gpart with ID modulo SWIFT_GRAVITY_FORCE_CHECKS will get their forces
- * computed.
- *
- * @param s The #space to use.
- * @param e The #engine (to access the current time).
+ * @brief Mapper function for the exact gravity calculation.
  */
-void gravity_exact_force_compute(struct space *s, const struct engine *e) {
-
+void gravity_exact_force_compute_mapper(void *map_data, int nr_gparts,
+                                        void *extra_data) {
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
 
-  const ticks tic = getticks();
-  const double const_G = e->physical_constants->const_newton_G;
+  /* Unpack the data */
+  struct gpart *restrict gparts = (struct gpart *)map_data;
+  struct exact_force_data *data = (struct exact_force_data *)extra_data;
+  const struct space *s = data->s;
+  const struct engine *e = data->e;
+  const double const_G = data->const_G;
   int counter = 0;
 
-  /* Let's start by checking whether we already computed these forces */
-  if (gravity_exact_force_file_exits(e)) {
-    message("Exact accelerations already computed. Skipping calculation.");
-    return;
-  }
-
-  /* No matching file present ? Do it then */
-  for (size_t i = 0; i < s->nr_gparts; ++i) {
+  for (int i = 0; i < nr_gparts; ++i) {
 
-    struct gpart *gpi = &s->gparts[i];
+    struct gpart *gpi = &gparts[i];
 
     /* Is the particle active and part of the subset to be tested ? */
     if (gpi->id_or_neg_offset % SWIFT_GRAVITY_FORCE_CHECKS == 0 &&
@@ -118,13 +116,13 @@ void gravity_exact_force_compute(struct space *s, const struct engine *e) {
       double a_grav[3] = {0., 0., 0.};
 
       /* Interact it with all other particles in the space.*/
-      for (size_t j = 0; j < s->nr_gparts; ++j) {
-
-        /* No self interaction */
-        if (i == j) continue;
+      for (int j = 0; j < (int)s->nr_gparts; ++j) {
 
         struct gpart *gpj = &s->gparts[j];
 
+        /* No self interaction */
+        if (gpi == gpj) continue;
+
         /* Compute the pairwise distance. */
         const double dx[3] = {gpi->x[0] - gpj->x[0],   // x
                               gpi->x[1] - gpj->x[1],   // y
@@ -173,9 +171,47 @@ void gravity_exact_force_compute(struct space *s, const struct engine *e) {
       counter++;
     }
   }
+  atomic_add(&data->counter_global, counter);
 
-  message("Computed exact gravity for %d gparts (took %.3f %s). ", counter,
-          clocks_from_ticks(getticks() - tic), clocks_getunit());
+#else
+  error("Gravity checking function called without the corresponding flag.");
+#endif
+}
+
+/**
+ * @brief Run a brute-force gravity calculation for a subset of particles.
+ *
+ * All gpart with ID modulo SWIFT_GRAVITY_FORCE_CHECKS will get their forces
+ * computed.
+ *
+ * @param s The #space to use.
+ * @param e The #engine (to access the current time).
+ */
+void gravity_exact_force_compute(struct space *s, const struct engine *e) {
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+
+  const ticks tic = getticks();
+
+  /* Let's start by checking whether we already computed these forces */
+  if (gravity_exact_force_file_exits(e)) {
+    message("Exact accelerations already computed. Skipping calculation.");
+    return;
+  }
+
+  /* No matching file present ? Do it then */
+  struct exact_force_data data;
+  data.e = e;
+  data.s = s;
+  data.counter_global = 0;
+  data.const_G = e->physical_constants->const_newton_G;
+
+  threadpool_map(&s->e->threadpool, gravity_exact_force_compute_mapper,
+                 s->gparts, s->nr_gparts, sizeof(struct gpart), 1000, &data);
+
+  message("Computed exact gravity for %d gparts (took %.3f %s). ",
+          data.counter_global, clocks_from_ticks(getticks() - tic),
+          clocks_getunit());
 
 #else
   error("Gravity checking function called without the corresponding flag.");

src/gravity_properties.c

@@ -59,7 +59,8 @@ void gravity_props_init(struct gravity_props *p,
 
 void gravity_props_print(const struct gravity_props *p) {
 
-  message("Self-gravity scheme: FMM-MM");
+  message("Self-gravity scheme: FMM-MM with m-poles of order %d",
+          SELF_GRAVITY_MULTIPOLE_ORDER);
 
   message("Self-gravity time integration: eta=%.4f", p->eta);
 
@@ -68,7 +69,7 @@ void gravity_props_print(const struct gravity_props *p) {
   message("Self-gravity softening:    epsilon=%.4f", p->epsilon);
 
   if (p->a_smooth != gravity_props_default_a_smooth)
-    message("Self-gravity smoothing-scale: a_smooth=%f", p->a_smooth);
+    message("Self-gravity MM smoothing-scale: a_smooth=%f", p->a_smooth);
 
   if (p->r_cut != gravity_props_default_r_cut)
     message("Self-gravity MM cut-off: r_cut=%f", p->r_cut);
@@ -81,6 +82,7 @@ void gravity_props_print_snapshot(hid_t h_grpgrav,
   io_write_attribute_f(h_grpgrav, "Time integration eta", p->eta);
   io_write_attribute_f(h_grpgrav, "Softening length", p->epsilon);
   io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit);
+  io_write_attribute_d(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
   io_write_attribute_f(h_grpgrav, "MM a_smooth", p->a_smooth);
   io_write_attribute_f(h_grpgrav, "MM r_cut", p->r_cut);
 }

src/multipole.h

@@ -176,9 +176,15 @@ struct gravity_tensors {
       /*! Centre of mass of the matter dsitribution */
       double CoM[3];
 
+      /*! Centre of mass of the matter dsitribution at the last rebuild */
+      double CoM_rebuild[3];
+
       /*! Upper limit of the CoM<->gpart distance */
       double r_max;
 
+      /*! Upper limit of the CoM<->gpart distance at the last rebuild */
+      double r_max_rebuild;
+
       /*! Multipole mass */
       struct multipole m_pole;
 
@@ -1232,12 +1238,10 @@ INLINE static void gravity_P2M(struct gravity_tensors *m,
  * @param m_b The #multipole to shift.
  * @param pos_a The position to which m_b will be shifted.
  * @param pos_b The current postion of the multipole to shift.
- * @param periodic Is the calculation periodic ?
  */
 INLINE static void gravity_M2M(struct multipole *m_a,
                                const struct multipole *m_b,
-                               const double pos_a[3], const double pos_b[3],
-                               int periodic) {
+                               const double pos_a[3], const double pos_b[3]) {
   /* Shift bulk velocity */
   m_a->vel[0] = m_b->vel[0];
   m_a->vel[1] = m_b->vel[1];
@@ -2521,17 +2525,24 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
  * @param ma The #multipole of the first #cell.
  * @param mb The #multipole of the second #cell.
  * @param theta_crit_inv The inverse of the critical opening angle.
+ * @param Are we using the current value of CoM or the ones from the last
+ * rebuild ?
  */
 __attribute__((always_inline)) INLINE static int gravity_multipole_accept(
     const struct gravity_tensors *ma, const struct gravity_tensors *mb,
-    double theta_crit_inv) {
-
-  const double r_crit_a = ma->r_max * theta_crit_inv;
-  const double r_crit_b = mb->r_max * theta_crit_inv;
-
-  const double dx = ma->CoM[0] - mb->CoM[0];
-  const double dy = ma->CoM[1] - mb->CoM[1];
-  const double dz = ma->CoM[2] - mb->CoM[2];
+    double theta_crit_inv, int rebuild) {
+
+  const double r_crit_a =
+      (rebuild ? ma->r_max_rebuild : ma->r_max) * theta_crit_inv;
+  const double r_crit_b =
+      (rebuild ? mb->r_max_rebuild : mb->r_max) * theta_crit_inv;
+
+  const double dx = rebuild ? ma->CoM_rebuild[0] - mb->CoM_rebuild[0]
+                            : ma->CoM[0] - mb->CoM[0];
+  const double dy = rebuild ? ma->CoM_rebuild[1] - mb->CoM_rebuild[1]
+                            : ma->CoM[1] - mb->CoM[1];
+  const double dz = rebuild ? ma->CoM_rebuild[2] - mb->CoM_rebuild[2]
+                            : ma->CoM[2] - mb->CoM[2];
 
   const double r2 = dx * dx + dy * dy + dz * dz;
 

src/runner.c

@@ -1356,9 +1356,8 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
 
       if (part_is_active(p, e)) {
 
-        /* First, finish the force loop */
+        /* Finish the force loop */
         hydro_end_force(p);
-        if (p->gpart != NULL) gravity_end_force(p->gpart, const_G);
       }
     }
 
@@ -1368,25 +1367,44 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
       /* Get a handle on the gpart. */
       struct gpart *restrict gp = &gparts[k];
 
-      if (gp->type == swift_type_dark_matter) {
+      if (gpart_is_active(gp, e)) {
 
-        if (gpart_is_active(gp, e)) gravity_end_force(gp, const_G);
-      }
+        /* Finish the force calculation */
+        gravity_end_force(gp, const_G);
+
+#ifdef SWIFT_NO_GRAVITY_BELOW_ID
+        /* Cancel gravity forces of these particles */
+        if ((gp->type == swift_type_dark_matter &&
+             gp->id_or_neg_offset < SWIFT_NO_GRAVITY_BELOW_ID) ||
+            (gp->type == swift_type_gas &&
+             parts[-gp->id_or_neg_offset].id < SWIFT_NO_GRAVITY_BELOW_ID) ||
+            (gp->type == swift_type_star &&
+             sparts[-gp->id_or_neg_offset].id < SWIFT_NO_GRAVITY_BELOW_ID)) {
+
+          /* Don't move ! */
+          gp->a_grav[0] = 0.f;
+          gp->a_grav[1] = 0.f;
+          gp->a_grav[2] = 0.f;
+        }
+#endif
 
 #ifdef SWIFT_DEBUG_CHECKS
-      if (e->policy & engine_policy_self_gravity && gpart_is_active(gp, e)) {
+        if (e->policy & engine_policy_self_gravity) {
 
-        /* Check that this gpart has interacted with all the other particles
-         * (via direct or multipoles) in the box */
-        gp->num_interacted++;
-        if (gp->num_interacted != (long long)e->s->nr_gparts)
-          error(
-              "g-particle (id=%lld, type=%d) did not interact gravitationally "
-              "with all other gparts gp->num_interacted=%lld, total_gparts=%zd",
-              gp->id_or_neg_offset, gp->type, gp->num_interacted,
-              e->s->nr_gparts);
-      }
+          /* Check that this gpart has interacted with all the other
+           * particles (via direct or multipoles) in the box */
+          gp->num_interacted++;
+          if (gp->num_interacted != (long long)e->s->nr_gparts)
+            error(
+                "g-particle (id=%lld, type=%d) did not interact "
+                "gravitationally "
+                "with all other gparts gp->num_interacted=%lld, "
+                "total_gparts=%zd",
+                gp->id_or_neg_offset, gp->type, gp->num_interacted,
+                e->s->nr_gparts);
+        }
 #endif
+      }
     }
 
     /* Loop over the star particles in this cell. */
@@ -1396,9 +1414,8 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
       struct spart *restrict sp = &sparts[k];