Sink drift

examples/IsolatedGalaxy/IsolatedGalaxy_sink/getIC.sh

 #!/bin/bash
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/3e11-star-only-DM-halo-galaxy.hdf5
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/$1
+

examples/IsolatedGalaxy/IsolatedGalaxy_sink/isolated_galaxy.yml

@@ -8,6 +8,9 @@ InternalUnitSystem:
 
 # Parameters for the self-gravity scheme
 Gravity:
+  MAC:                           adaptive  # Choice of mulitpole acceptance criterion: 'adaptive' OR 'geometric'.
+  epsilon_fmm:                   0.001     # Tolerance parameter for the adaptive multipole acceptance criterion.
+  theta_cr:                      0.7       # Opening angle for the purely gemoetric criterion.
   eta:          0.025               # Constant dimensionless multiplier for time integration.
   theta:        0.7                 # Opening angle (Multipole acceptance criterion).
   max_physical_DM_softening: 0.7    # Physical softening length (in internal units).
@@ -39,4 +42,41 @@ InitialConditions:
   file_name:  isolated_galaxy.hdf5  # The file to read
   periodic:                    0    # Are we running with periodic ICs?
 
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:                    1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:                     0.1      # Courant-Friedrich-Levy condition for time integration.
+  use_mass_weighted_num_ngb:         0        # (Optional) Are we using the mass-weighted definition of the number of neighbours?
+  h_tolerance:                       1e-4     # (Optional) Relative accuracy of the Netwon-Raphson scheme for the smoothing lengths.
+  h_max:                             10.      # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
+  h_min_ratio:                       0.       # (Optional) Minimal allowed smoothing length in units of the softening. Defaults to 0 if unspecified.
+  max_volume_change:                 1.4      # (Optional) Maximal allowed change of kernel volume over one time-step.
+  max_ghost_iterations:              30       # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
+  particle_splitting:                1        # (Optional) Are we splitting particles that are too massive (default: 0)
+  particle_splitting_mass_threshold: 7e-4     # (Optional) Mass threshold for particle splitting (in internal units)
+  generate_random_ids:               0        # (Optional) When creating new particles via splitting, generate ids at random (1) or use new IDs beyond the current range (0) (default: 0)
+  initial_temperature:               0        # (Optional) Initial temperature (in internal units) to set the gas particles at start-up. Value is ignored if set to 0.
+  minimal_temperature:               0        # (Optional) Minimal temperature (in internal units) allowed for the gas particles. Value is ignored if set to 0.
+  H_mass_fraction:                   0.755    # (Optional) Hydrogen mass fraction used for initial conversion from temp to internal energy. Default value is derived from the physical constants.
+  H_ionization_temperature:          1e4      # (Optional) Temperature of the transition from neutral to ionized Hydrogen for primoridal gas.
+  viscosity_alpha:                   0.8      # (Optional) Override for the initial value of the artificial viscosity. In schemes that have a fixed AV, this remains as alpha throughout the run.
+  viscosity_alpha_max:               2.0      # (Optional) Maximal value for the artificial viscosity in schemes that allow alpha to vary.
+  viscosity_alpha_min:               0.1      # (Optional) Minimal value for the artificial viscosity in schemes that allow alpha to vary.
+  viscosity_length:                  0.1      # (Optional) Decay length for the artificial viscosity in schemes that allow alpha to vary.
+  diffusion_alpha:                   0.0      # (Optional) Override the initial value for the thermal diffusion coefficient in schemes with thermal diffusion.
+  diffusion_beta:                    0.01     # (Optional) Override the decay/rise rate tuning parameter for the thermal diffusion.
+  diffusion_alpha_max:               1.0      # (Optional) Override the maximal thermal diffusion coefficient that is allowed for a given particle.
+  diffusion_alpha_min:               0.0      # (Optional) Override the minimal thermal diffusion coefficient that is allowed for a given particle.
  
+# Hernquist potential parameters
+HernquistPotential:
+  useabspos:       0        # 0 -> positions based on centre, 1 -> absolute positions 
+  position:        [0.,0.,0.]    # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
+  idealizeddisk:   1        # Run with an idealized galaxy disk
+  M200:            137.0   # M200 of the galaxy disk
+  h:               0.704    # reduced Hubble constant (value does not specify the used units!)
+  concentration:   9.0      # concentration of the Halo
+  diskfraction:              0.040   # Disk mass fraction
+  bulgefraction:              0.014   # Bulge mass fraction
+  timestep_mult:   0.01     # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
+  epsilon:         0.2      # Softening size (internal units)

examples/IsolatedGalaxy/IsolatedGalaxy_sink/makeIC.py

@@ -20,14 +20,14 @@
 import h5py
 import numpy as np
 from shutil import copyfile
+import sys
 
 # Add sink particles to the isolated galaxy
 
-fileName = "3e11-star-only-DM-halo-galaxy.hdf5"
+fileName = sys.argv[-1]
 output = "isolated_galaxy.hdf5"
 
-L = 13756  # Size of the box (internal units)
-L_sink = 1000 # Size of the sink particle area (L_sink < L)
+L_sink = 50  # Size of the sink particle area
 N = 1000  # Number of sink particles
 max_vel = 100  # Maximal velocity of the sink particles (in km / s)
 mass = 0.000142  # Mass of a sink particle (internal units)
@@ -40,6 +40,11 @@ copyfile(fileName, output)
 # File
 file = h5py.File(output, 'a')
 
+# Get the box size (suppose that the galaxy is at the center of the box)
+L = 2 * file["PartType0/Coordinates"][:].mean()
+if (L < L_sink):
+    raise Exception("Error the size of the sink box is too large")
+
 pos = np.random.rand(N, 3) * L_sink + 0.5 * (L - L_sink)
 vel = 2 * (np .random.rand(N, 3) - 0.5) * max_vel
 m = mass * np.ones(N)

examples/IsolatedGalaxy/IsolatedGalaxy_sink/run.sh

 #!/bin/bash
 
-if [ ! -e 3e11-star-only-DM-halo-galaxy.hdf5 ]
+# filename=fid.hdf5
+# filename=f10.hdf5
+# filename=f90.hdf5
+filename=lowres8.hdf5
+# filename=lowres64.hdf5
+# filename=lowres512.hdf5
+# filename=highres6.hdf5
+
+if [ ! -e $filename ]
 then
     echo "Fetching initial conditons for the isolated galaxy with an external potential ..."
-    ./getIC.sh
+    ./getIC.sh $filename
 fi
 
 echo "Generate initial conditions"
-python3 makeIC.py
+python3 makeIC.py $filename
 
-../../swift --sinks --external-gravity --self-gravity --threads=16 isolated_galaxy.yml 2>&1 | tee output.log
+../../swift --hydro --sinks --external-gravity --self-gravity --threads=1 isolated_galaxy.yml 2>&1 | tee output.log

examples/main.c

@@ -1298,7 +1298,7 @@ int main(int argc, char *argv[]) {
     if (with_fof) engine_policies |= engine_policy_fof;
     if (with_logger) engine_policies |= engine_policy_logger;
     if (with_line_of_sight) engine_policies |= engine_policy_line_of_sight;
-    if (with_sink) engine_policies |= engine_policy_sink;
+    if (with_sink) engine_policies |= engine_policy_sinks;
     if (with_rt) engine_policies |= engine_policy_rt;
 
     /* Initialize the engine with the space and policies. */
@@ -1401,10 +1401,12 @@ int main(int argc, char *argv[]) {
 
   /* Legend */
   if (myrank == 0) {
-    printf("# %6s %14s %12s %12s %14s %9s %12s %12s %12s %12s %16s [%s] %6s\n",
-           "Step", "Time", "Scale-factor", "Redshift", "Time-step", "Time-bins",
-           "Updates", "g-Updates", "s-Updates", "b-Updates", "Wall-clock time",
-           clocks_getunit(), "Props");
+    printf(
+        "# %6s %14s %12s %12s %14s %9s %12s %12s %12s %12s %12s %16s [%s] "
+        "%6s\n",
+        "Step", "Time", "Scale-factor", "Redshift", "Time-step", "Time-bins",
+        "Updates", "g-Updates", "s-Updates", "sink-Updates", "b-Updates",
+        "Wall-clock time", clocks_getunit(), "Props");
     fflush(stdout);
   }
 
@@ -1548,19 +1550,21 @@ int main(int argc, char *argv[]) {
 
     /* Print some information to the screen */
     printf(
-        "  %6d %14e %12.7f %12.7f %14e %4d %4d %12lld %12lld %12lld %12lld"
+        "  %6d %14e %12.7f %12.7f %14e %4d %4d %12lld %12lld %12lld %12lld "
+        "%12lld"
         " %21.3f %6d\n",
         e.step, e.time, e.cosmology->a, e.cosmology->z, e.time_step,
         e.min_active_bin, e.max_active_bin, e.updates, e.g_updates, e.s_updates,
-        e.b_updates, e.wallclock_time, e.step_props);
+        e.sink_updates, e.b_updates, e.wallclock_time, e.step_props);
     fflush(stdout);
 
     fprintf(e.file_timesteps,
             "  %6d %14e %12.7f %12.7f %14e %4d %4d %12lld %12lld %12lld %12lld"
-            " %21.3f %6d\n",
+            " %12lld %21.3f %6d\n",
             e.step, e.time, e.cosmology->a, e.cosmology->z, e.time_step,
             e.min_active_bin, e.max_active_bin, e.updates, e.g_updates,
-            e.s_updates, e.b_updates, e.wallclock_time, e.step_props);
+            e.s_updates, e.sink_updates, e.b_updates, e.wallclock_time,
+            e.step_props);
     fflush(e.file_timesteps);
 
     /* Print information to the SFH logger */

src/cell.c

@@ -2065,6 +2065,7 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
   /* Store the counts and offsets. */
   for (int k = 0; k < 8; k++) {
     c->progeny[k]->sinks.count = bucket_count[k];
+    c->progeny[k]->sinks.count_total = c->progeny[k]->sinks.count;
     c->progeny[k]->sinks.parts = &c->sinks.parts[bucket_offset[k]];
   }
 
@@ -2294,6 +2295,44 @@ void cell_check_gpart_drift_point(struct cell *c, void *data) {
 #endif
 }
 
+/**
+ * @brief Checks that the #sink in a cell are at the
+ * current point in time
+ *
+ * Calls error() if the cell is not at the current time.
+ *
+ * @param c Cell to act upon
+ * @param data The current time on the integer time-line
+ */
+void cell_check_sink_drift_point(struct cell *c, void *data) {
+#ifdef SWIFT_DEBUG_CHECKS
+
+  const integertime_t ti_drift = *(integertime_t *)data;
+
+  /* Only check local cells */
+  if (c->nodeID != engine_rank) return;
+
+  /* Only check cells with content */
+  if (c->sinks.count == 0) return;
+
+  if (c->sinks.ti_old_part != ti_drift)
+    error(
+        "Cell in an incorrect time-zone! c->sinks.ti_old_part=%lld "
+        "ti_drift=%lld",
+        c->sinks.ti_old_part, ti_drift);
+
+  for (int i = 0; i < c->sinks.count; ++i)
+    if (c->sinks.parts[i].ti_drift != ti_drift &&
+        c->sinks.parts[i].time_bin != time_bin_inhibited)
+      error(
+          "sink-part in an incorrect time-zone! sink->ti_drift=%lld "
+          "ti_drift=%lld",
+          c->sinks.parts[i].ti_drift, ti_drift);
+#else
+  error("Calling debugging code without debugging flag activated.");
+#endif
+}
+
 /**
  * @brief Checks that the #spart in a cell are at the
  * current point in time
@@ -2633,7 +2672,8 @@ void cell_clear_drift_flags(struct cell *c, void *data) {
                          cell_flag_do_grav_drift | cell_flag_do_grav_sub_drift |
                          cell_flag_do_bh_drift | cell_flag_do_bh_sub_drift |
                          cell_flag_do_stars_drift |
-                         cell_flag_do_stars_sub_drift);
+                         cell_flag_do_stars_sub_drift |
+                         cell_flag_do_sink_drift | cell_flag_do_sink_sub_drift);
 }
 
 /**
@@ -3027,6 +3067,44 @@ void cell_activate_drift_bpart(struct cell *c, struct scheduler *s) {
   }
 }
 
+/**
+ * @brief Activate the #sink drifts on the given cell.
+ */
+void cell_activate_drift_sink(struct cell *c, struct scheduler *s) {
+
+  /* If this cell is already marked for drift, quit early. */
+  if (cell_get_flag(c, cell_flag_do_sink_drift)) return;
+
+  /* Mark this cell for drifting. */
+  cell_set_flag(c, cell_flag_do_sink_drift);
+
+  /* Set the do_sink_sub_drifts all the way up and activate the super
+     drift if this has not yet been done. */
+  if (c == c->hydro.super) {
+#ifdef SWIFT_DEBUG_CHECKS
+    if (c->sinks.drift == NULL)
+      error("Trying to activate un-existing c->sinks.drift");
+#endif
+    scheduler_activate(s, c->sinks.drift);
+  } else {
+    for (struct cell *parent = c->parent;
+         parent != NULL && !cell_get_flag(parent, cell_flag_do_sink_sub_drift);
+         parent = parent->parent) {
+      /* Mark this cell for drifting */
+      cell_set_flag(parent, cell_flag_do_sink_sub_drift);
+
+      if (parent == c->hydro.super) {
+#ifdef SWIFT_DEBUG_CHECKS
+        if (parent->sinks.drift == NULL)
+          error("Trying to activate un-existing parent->sinks.drift");
+#endif
+        scheduler_activate(s, parent->sinks.drift);
+        break;
+      }
+    }
+  }
+}
+
 /**
  * @brief Activate the drifts on the given cell.
  */
@@ -3514,6 +3592,102 @@ void cell_activate_subcell_black_holes_tasks(struct cell *ci, struct cell *cj,
   } /* Otherwise, pair interation */
 }
 
+/**
+ * @brief Traverse a sub-cell task and activate the sinks drift tasks that
+ * are required by a sinks task
+ *
+ * @param ci The first #cell we recurse in.
+ * @param cj The second #cell we recurse in.
+ * @param s The task #scheduler.
+ * @param with_timestep_sync Are we running with time-step synchronization on?
+ */
+void cell_activate_subcell_sinks_tasks(struct cell *ci, struct cell *cj,
+                                       struct scheduler *s,
+                                       const int with_timestep_sync) {
+  const struct engine *e = s->space->e;
+
+  /* Store the current dx_max and h_max values. */
+  ci->sinks.dx_max_part_old = ci->sinks.dx_max_part;
+  ci->hydro.dx_max_part_old = ci->hydro.dx_max_part;
+  ci->hydro.h_max_old = ci->hydro.h_max;
+
+  if (cj != NULL) {
+    cj->sinks.dx_max_part_old = cj->sinks.dx_max_part;
+    cj->hydro.dx_max_part_old = cj->hydro.dx_max_part;
+    cj->hydro.h_max_old = cj->hydro.h_max;
+  }
+
+  /* Self interaction? */
+  if (cj == NULL) {
+    /* Do anything? */
+    if (!cell_is_active_sinks(ci, e) || ci->hydro.count == 0 ||
+        ci->sinks.count == 0)
+      return;
+
+    /* Recurse? */
+    if (cell_can_recurse_in_self_sinks_task(ci)) {
+      /* Loop over all progenies and pairs of progenies */
+      for (int j = 0; j < 8; j++) {
+        if (ci->progeny[j] != NULL) {
+          cell_activate_subcell_sinks_tasks(ci->progeny[j], NULL, s,
+                                            with_timestep_sync);
+          for (int k = j + 1; k < 8; k++)
+            if (ci->progeny[k] != NULL)
+              cell_activate_subcell_sinks_tasks(ci->progeny[j], ci->progeny[k],
+                                                s, with_timestep_sync);
+        }
+      }
+    } else {
+      /* We have reached the bottom of the tree: activate drift */
+      cell_activate_drift_sink(ci, s);
+      cell_activate_drift_part(ci, s);
+    }
+  }
+
+  /* Otherwise, pair interation */
+  else {
+    error("Not implemented yet.");
+    /* /\* Should we even bother? *\/ */
+    /* if (!cell_is_active_sinks(ci, e) && */
+    /*     !cell_is_active_sinks(cj, e)) */
+    /*   return; */
+
+    /* /\* Get the orientation of the pair. *\/ */
+    /* double shift[3]; */
+    /* const int sid = space_getsid(s->space, &ci, &cj, shift); */
+
+    /* /\* recurse? *\/ */
+    /* if (cell_can_recurse_in_pair_sinks_task(ci, cj) && */
+    /*     cell_can_recurse_in_pair_sinks_task(cj, ci)) { */
+    /*   const struct cell_split_pair *csp = &cell_split_pairs[sid]; */
+    /*   for (int k = 0; k < csp->count; k++) { */
+    /*     const int pid = csp->pairs[k].pid; */
+    /*     const int pjd = csp->pairs[k].pjd; */
+    /*     if (ci->progeny[pid] != NULL && cj->progeny[pjd] != NULL) */
+    /*       cell_activate_subcell_sinks_tasks( */
+    /*           ci->progeny[pid], cj->progeny[pjd], s, with_timestep_sync); */
+    /*   } */
+    /* } */
+
+    /* /\* Otherwise, activate the drifts. *\/ */
+    /* else if (cell_is_active_sinks(ci, e) || */
+    /*          cell_is_active_sinks(cj, e)) { */
+
+    /*   /\* Activate the drifts if the cells are local. *\/ */
+    /*   if (ci->nodeID == engine_rank) cell_activate_drift_sinks(ci, s); */
+    /*   if (cj->nodeID == engine_rank) cell_activate_drift_part(cj, s); */
+    /*   if (cj->nodeID == engine_rank && with_timestep_sync) */
+    /*     cell_activate_sync_part(cj, s); */
+
+    /*   /\* Activate the drifts if the cells are local. *\/ */
+    /*   if (ci->nodeID == engine_rank) cell_activate_drift_part(ci, s); */
+    /*   if (cj->nodeID == engine_rank) cell_activate_drift_sinks(cj, s); */
+    /*   if (ci->nodeID == engine_rank && with_timestep_sync) */
+    /*     cell_activate_sync_part(ci, s); */
+    /* } */
+  } /* Otherwise, pair interation */
+}
+
 /**
  * @brief Traverse a sub-cell task and activate the gravity drift tasks that
  * are required by a self gravity task.
@@ -4609,6 +4783,38 @@ int cell_unskip_black_holes_tasks(struct cell *c, struct scheduler *s) {
   return rebuild;
 }
 
+/**
+ * @brief Un-skips all the sinks tasks associated with a given cell and
+ * checks if the space needs to be rebuilt.
+ *
+ * @param c the #cell.
+ * @param s the #scheduler.
+ *
+ * @return 1 If the space needs rebuilding. 0 otherwise.
+ */
+int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
+
+  struct engine *e = s->space->e;
+  const int nodeID = e->nodeID;
+  int rebuild = 0;
+
+  if (c->sinks.drift != NULL && cell_is_active_sinks(c, e)) {
+    cell_activate_drift_sink(c, s);
+  }
+
+  /* Unskip all the other task types. */
+  if (c->nodeID == nodeID && cell_is_active_sinks(c, e)) {
+    if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
+    if (c->kick2 != NULL) scheduler_activate(s, c->kick2);
+    if (c->timestep != NULL) scheduler_activate(s, c->timestep);
+#ifdef WITH_LOGGER
+    if (c->logger != NULL) scheduler_activate(s, c->logger);
+#endif
+  }
+
+  return rebuild;
+}
+
 /**
  * @brief Set the super-cell pointers for all cells in a hierarchy.
  *
@@ -5439,6 +5645,156 @@ void cell_drift_bpart(struct cell *c, const struct engine *e, int force) {
   cell_clear_flag(c, cell_flag_do_bh_drift | cell_flag_do_bh_sub_drift);
 }
 
+/**
+ * @brief Recursively drifts the #sinks in a cell hierarchy.
+ *
+ * @param c The #cell.
+ * @param e The #engine (to get ti_current).
+ * @param force Drift the particles irrespective of the #cell flags.
+ */
+void cell_drift_sink(struct cell *c, const struct engine *e, int force) {
+
+  const int periodic = e->s->periodic;
+  const double dim[3] = {e->s->dim[0], e->s->dim[1], e->s->dim[2]};
+  const int with_cosmology = (e->policy & engine_policy_cosmology);
+  const integertime_t ti_old_bpart = c->sinks.ti_old_part;
+  const integertime_t ti_current = e->ti_current;
+  struct sink *const sinks = c->sinks.parts;
+
+  float dx_max = 0.f, dx2_max = 0.f;
+
+  /* Drift irrespective of cell flags? */
+  force = (force || cell_get_flag(c, cell_flag_do_sink_drift));
+
+#ifdef SWIFT_DEBUG_CHECKS
+  /* Check that we only drift local cells. */
+  if (c->nodeID != engine_rank) error("Drifting a foreign cell is nope.");
+
+  /* Check that we are actually going to move forward. */
+  if (ti_current < ti_old_bpart) error("Attempt to drift to the past");
+#endif
+
+  /* Early abort? */
+  if (c->sinks.count == 0) {
+
+    /* Clear the drift flags. */
+    cell_clear_flag(c, cell_flag_do_sink_drift | cell_flag_do_sink_sub_drift);
+
+    /* Update the time of the last drift */
+    c->sinks.ti_old_part = ti_current;
+
+    return;
+  }
+
+  /* Ok, we have some particles somewhere in the hierarchy to drift */
+
+  /* Are we not in a leaf ? */
+  if (c->split && (force || cell_get_flag(c, cell_flag_do_sink_sub_drift))) {
+
+    /* Loop over the progeny and collect their data. */
+    for (int k = 0; k < 8; k++) {
+      if (c->progeny[k] != NULL) {
+        struct cell *cp = c->progeny[k];
+
+        /* Recurse */
+        cell_drift_sink(cp, e, force);
+
+        /* Update */
+        dx_max = max(dx_max, cp->sinks.dx_max_part);
+      }
+    }
+
+    /* Store the values */
+    c->sinks.dx_max_part = dx_max;
+
+    /* Update the time of the last drift */
+    c->sinks.ti_old_part = ti_current;
+
+  } else if (!c->split && force && ti_current > ti_old_bpart) {
+
+    /* Drift from the last time the cell was drifted to the current time */
+    double dt_drift;
+    if (with_cosmology) {
+      dt_drift =
+          cosmology_get_drift_factor(e->cosmology, ti_old_bpart, ti_current);
+    } else {
+      dt_drift = (ti_current - ti_old_bpart) * e->time_base;
+    }
+
+    /* Loop over all the star particles in the cell */
+    const size_t nr_sinks = c->sinks.count;
+    for (size_t k = 0; k < nr_sinks; k++) {
+
+      /* Get a handle on the sink. */
+      struct sink *const sink = &sinks[k];
+
+      /* Ignore inhibited particles */
+      if (sink_is_inhibited(sink, e)) continue;
+
+      /* Drift... */
+      drift_sink(sink, dt_drift, ti_old_bpart, ti_current);
+
+#ifdef SWIFT_DEBUG_CHECKS
+      /* Make sure the particle does not drift by more than a box length. */
+      if (fabs(sink->v[0] * dt_drift) > e->s->dim[0] ||
+          fabs(sink->v[1] * dt_drift) > e->s->dim[1] ||
+          fabs(sink->v[2] * dt_drift) > e->s->dim[2]) {
+        error("Particle drifts by more than a box length!");
+      }
+#endif
+
+      /* In non-periodic BC runs, remove particles that crossed the border */
+      if (!periodic) {
+
+        /* Did the particle leave the box?  */
+        if ((sink->x[0] > dim[0]) || (sink->x[0] < 0.) ||  // x
+            (sink->x[1] > dim[1]) || (sink->x[1] < 0.) ||  // y
+            (sink->x[2] > dim[2]) || (sink->x[2] < 0.)) {  // z
+
+          lock_lock(&e->s->lock);
+
+          /* Re-check that the particle has not been removed
+           * by another thread before we do the deed. */
+          if (!sink_is_inhibited(sink, e)) {
+
+            /* Remove the particle entirely */
+            // cell_remove_sink(e, c, bp);
+            error("TODO: loic implement cell_remove_sink");
+          }
+
+          if (lock_unlock(&e->s->lock) != 0)
+            error("Failed to unlock the space!");
+
+          continue;
+        }
+      }
+
+      /* Compute (square of) motion since last cell construction */
+      const float dx2 = sink->x_diff[0] * sink->x_diff[0] +
+                        sink->x_diff[1] * sink->x_diff[1] +
+                        sink->x_diff[2] * sink->x_diff[2];
+      dx2_max = max(dx2_max, dx2);
+
+      /* Get ready for a density calculation */
+      if (sink_is_active(sink, e)) {
+        sink_init_sink(sink);
+      }
+    }
+
+    /* Now, get the maximal particle motion from its square */
+    dx_max = sqrtf(dx2_max);
+
+    /* Store the values */
+    c->sinks.dx_max_part = dx_max;
+
+    /* Update the time of the last drift */
+    c->sinks.ti_old_part = ti_current;
+  }
+
+  /* Clear the drift flags. */
+  cell_clear_flag(c, cell_flag_do_sink_drift | cell_flag_do_sink_sub_drift);
+}
+
 /**
  * @brief Recursively drifts all multipoles in a cell hierarchy.
  *
@@ -5587,7 +5943,7 @@ void cell_check_timesteps(const struct cell *c, const integertime_t ti_current,
 
   if (c->hydro.ti_end_min == 0 && c->grav.ti_end_min == 0 &&
       c->stars.ti_end_min == 0 && c->black_holes.ti_end_min == 0 &&
-      c->nr_tasks > 0)
+      c->sinks.ti_end_min == 0 && c->nr_tasks > 0)
     error("Cell without assigned time-step");
 
   if (c->split) {

src/cell.h

@@ -304,10 +304,12 @@ enum cell_flags {
   cell_flag_do_stars_sub_drift = (1UL << 10),
   cell_flag_do_bh_drift = (1UL << 11),
   cell_flag_do_bh_sub_drift = (1UL << 12),
-  cell_flag_do_stars_resort = (1UL << 13),
-  cell_flag_has_tasks = (1UL << 14),
-  cell_flag_do_hydro_sync = (1UL << 15),
-  cell_flag_do_hydro_sub_sync = (1UL << 16)
+  cell_flag_do_sink_drift = (1UL << 13),
+  cell_flag_do_sink_sub_drift = (1UL << 14),
+  cell_flag_do_stars_resort = (1UL << 15),
+  cell_flag_has_tasks = (1UL << 16),
+  cell_flag_do_hydro_sync = (1UL << 17),
+  cell_flag_do_hydro_sub_sync = (1UL << 18)
 };
 
 /**
@@ -317,7 +319,7 @@ enum cell_flags {
  */
 struct cell {
 
-  /*! The cell location on the grid. */
+  /*! The cell location on the grid (corner nearest to the origin). */
   double loc[3];
 
   /*! The cell dimensions. */
@@ -764,12 +766,21 @@ struct cell {
     /*! Nr of #sink this cell can hold after addition of new one. */
     int count_total;