Added the option to zero the gravity accelerations of particles below a given ID.

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

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,45 @@ 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 +1415,8 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
       struct spart *restrict sp = &sparts[k];
       if (spart_is_active(sp, e)) {
 
-        /* First, finish the force loop */
+        /* Finish the force loop */
         star_end_force(sp);
-        gravity_end_force(sp->gpart, const_G);
       }
     }
   }