Add fixed boundary particles to the code (boundary particles that will never have a velocity)

configure.ac

@@ -409,6 +409,22 @@ elif test "$boundary_particles" != "no"; then
    AC_DEFINE_UNQUOTED([SWIFT_BOUNDARY_PARTICLES], [$enableval] ,[Particles with smaller ID than this will be considered as boundaries.])
 fi
 
+# Check if we want to use fixed boundary particles.
+AC_ARG_ENABLE([fixed-boundary-particles],
+   [AS_HELP_STRING([--enable-fixed-boundary-particles],
+     [Set all particles with an ID smaller than @<:@N@:>@ as fixed boundary particles (i.e. receive zero gravity + hydro forces + zero velocity), this mode enables also --enable-boundary-particles and --enable-no-gravity-below-id.]
+   )],
+   [fixed_boundary_particles="$enableval"],
+   [fixed_boundary_particles="no"]
+)
+if test "$fixed_boundary_particles" == "yes"; then
+   AC_MSG_ERROR(Need to specify the ID below which particles get zero forces when using --enable-boundary-particles!)
+elif test "$fixed_boundary_particles" != "no"; then
+   AC_DEFINE_UNQUOTED([SWIFT_NO_GRAVITY_BELOW_ID], [$enableval] ,[Particles with smaller ID than this will have zero gravity forces])
+   AC_DEFINE_UNQUOTED([SWIFT_BOUNDARY_PARTICLES], [$enableval] ,[Particles with smaller ID than this will be considered as boundaries.])
+   AC_DEFINE_UNQUOTED([SWIFT_FIXED_BOUNDARY_PARTICLES], [$enableval] ,[Particles with smaller ID than this will be considered as boundaries.])
+fi
+
 # Check whether we have any of the ARM v8.1 tick timers
 AX_ASM_ARM_PMCCNTR
 AX_ASM_ARM_CNTVCT
@@ -2395,6 +2411,7 @@ AC_MSG_RESULT([
    Gravity checks              : $gravity_force_checks
    Custom icbrtf               : $enable_custom_icbrtf
    Boundary particles          : $boundary_particles
+   Fixed boundary particles    : $fixed_boundary_particles
 
    Particle Logger      : $with_logger
    Python enabled       : $have_python

doc/RTD/source/GettingStarted/special_modes.rst

@@ -29,6 +29,31 @@ in the sampling of the halo is reduced.
 Note also that this does not affect the hydrodynamic forces. This mode is
 purely designed for gravity-only accuracy tests.
 
+Besides the pure gravity mode, swift also has the boundary particle mode,
+this mode turns off both the gravity forces and hydro forces for all
+particles. Because gas particles only receive hydro this mode only impacts
+gas particles more strictly than other particles. This mode can be
+activated using ``--enable-boundary-particles=N``. This will zero the
+gravitational and hydrodynamic *accelerations* of all particles with ``id``
+(strictly) lower than ``N`` at every time-step. Still if particles have an
+initial velocity they will keep moving at that speed. This compilation
+option also activates ``--enable-no-gravity-below-id=N``. 
+
+Typical use cases are hydrodynamical experiments that have boundaries. 
+
+Both options discussed above only cancel accelerations and have no impact
+on what can happen in any code module that directly changes the velocity of
+the boundary or no gravity particles. An example of this is momentum
+injection of stellar winds for example. If we additionally want to keep the
+boundary particles fixed at the same position for the whole simulation we can
+use the ``--enable-fixed-boundary-particles=N`` compile option, this option
+explicitly sets the velocity of the boundary particles to zero every time
+step on top of also zeroing the accelerations.
+
+A typical use cases is an idealized setup in which we have a black hole in
+the centre of a galaxy that accretes gas but is not allowed to move from
+the momentum recoil of the gas it swallows.
+
 Gravity glasses
 ~~~~~~~~~~~~~~~
 

src/cell.c

@@ -5032,7 +5032,7 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
       if (gpart_is_inhibited(gp, e)) continue;
 
       /* Drift... */
-      drift_gpart(gp, dt_drift, ti_old_gpart, ti_current, grav_props);
+      drift_gpart(gp, dt_drift, ti_old_gpart, ti_current, grav_props, e);
 
 #ifdef SWIFT_DEBUG_CHECKS
       /* Make sure the particle does not drift by more than a box length. */

src/drift.h

@@ -27,6 +27,7 @@
 #include "const.h"
 #include "debug.h"
 #include "dimension.h"
+#include "engine.h"
 #include "entropy_floor.h"
 #include "hydro.h"
 #include "hydro_properties.h"
@@ -42,10 +43,12 @@
  * @param ti_old Integer start of time-step (for debugging checks).
  * @param ti_current Integer end of time-step (for debugging checks).
  * @param grav_props The properties of the gravity scheme.
+ * @param e the #engine
  */
 __attribute__((always_inline)) INLINE static void drift_gpart(
     struct gpart *restrict gp, double dt_drift, integertime_t ti_old,
-    integertime_t ti_current, const struct gravity_props *grav_props) {
+    integertime_t ti_current, const struct gravity_props *grav_props,
+    const struct engine *e) {
 
 #ifdef SWIFT_DEBUG_CHECKS
   if (gp->ti_drift != ti_old)
@@ -58,6 +61,29 @@ __attribute__((always_inline)) INLINE static void drift_gpart(
   gp->ti_drift = ti_current;
 #endif
 
+#ifdef SWIFT_FIXED_BOUNDARY_PARTICLES
+
+  /* Get the ID of the gpart */
+  long long id = 0;
+  if (gp->type == swift_type_gas)
+    id = e->s->parts[-gp->id_or_neg_offset].id;
+  else if (gp->type == swift_type_stars)
+    id = e->s->sparts[-gp->id_or_neg_offset].id;
+  else if (gp->type == swift_type_black_hole)
+    id = e->s->bparts[-gp->id_or_neg_offset].id;
+  else
+    id = gp->id_or_neg_offset;
+
+  /* Cancel the velocity of the particles */
+  if (id < SWIFT_FIXED_BOUNDARY_PARTICLES) {
+
+    /* Don't move! */
+    gp->v_full[0] = 0.f;
+    gp->v_full[1] = 0.f;
+    gp->v_full[2] = 0.f;
+  }
+#endif
+
   /* Drift... */
   gp->x[0] += gp->v_full[0] * dt_drift;
   gp->x[1] += gp->v_full[1] * dt_drift;
@@ -98,6 +124,21 @@ __attribute__((always_inline)) INLINE static void drift_part(
   p->ti_drift = ti_current;
 #endif
 
+#ifdef SWIFT_FIXED_BOUNDARY_PARTICLES
+
+  /* Get the ID of the gpart */
+  const long long id = p->id;
+
+  /* Cancel the velocity of the particles */
+  if (id < SWIFT_FIXED_BOUNDARY_PARTICLES) {
+
+    /* Don't move! */
+    xp->v_full[0] = 0.f;
+    xp->v_full[1] = 0.f;
+    xp->v_full[2] = 0.f;
+  }
+#endif
+
   /* Drift... */
   p->x[0] += xp->v_full[0] * dt_drift;
   p->x[1] += xp->v_full[1] * dt_drift;
@@ -121,6 +162,16 @@ __attribute__((always_inline)) INLINE static void drift_part(
     xp->x_diff[k] -= dx;
     xp->x_diff_sort[k] -= dx;
   }
+
+#ifdef SWIFT_FIXED_BOUNDARY_PARTICLES
+
+  /* Cancel the velocity of the particles */
+  if (id < SWIFT_FIXED_BOUNDARY_PARTICLES) {
+    p->v[0] = 0.f;
+    p->v[1] = 0.f;
+    p->v[2] = 0.f;
+  }
+#endif
 }
 
 /**
@@ -146,6 +197,21 @@ __attribute__((always_inline)) INLINE static void drift_spart(
   sp->ti_drift = ti_current;
 #endif
 
+#ifdef SWIFT_FIXED_BOUNDARY_PARTICLES
+
+  /* Get the ID of the gpart */
+  const long long id = sp->id;
+
+  /* Cancel the velocity of the particles */
+  if (id < SWIFT_FIXED_BOUNDARY_PARTICLES) {
+
+    /* Don't move! */
+    sp->v[0] = 0.f;
+    sp->v[1] = 0.f;
+    sp->v[2] = 0.f;
+  }
+#endif
+
   /* Drift... */
   sp->x[0] += sp->v[0] * dt_drift;
   sp->x[1] += sp->v[1] * dt_drift;
@@ -185,6 +251,21 @@ __attribute__((always_inline)) INLINE static void drift_bpart(
   bp->ti_drift = ti_current;
 #endif
 
+#ifdef SWIFT_FIXED_BOUNDARY_PARTICLES
+
+  /* Get the ID of the gpart */
+  const long long id = bp->id;
+
+  /* Cancel the velocity of the particles */
+  if (id < SWIFT_FIXED_BOUNDARY_PARTICLES) {
+
+    /* Don't move! */
+    bp->v[0] = 0.f;
+    bp->v[1] = 0.f;
+    bp->v[2] = 0.f;
+  }
+#endif
+
   /* Drift... */
   bp->x[0] += bp->v[0] * dt_drift;
   bp->x[1] += bp->v[1] * dt_drift;