Merge branch 'star_particles' into 'master'

Star particles and gparts links over MPI

A big bundle of changes:

 - Introduce star particles. At the moment they are empty shells but they are ready to be used.
 - An `spart` is always linked to a `gpart` and will receive its gravitational forces from there.
 - A `gpart` now contains an enum containing information about whether it is linked to a `part`, to an `spart` or to nothing.
 - Added a more comprehensive test to check that the part<->gpart and spart<->gpart links are correct after every rebuild or MPI transaction.
 - `spart` get kicked and drifted and get their time-step computed. 
 - `spart` are read-in and written to snapshots (only if running with `-S` for now). 
 - Only the particle types that will be used are now read from ICs. Speeds-up the starting up process.
 - Corrects many incorrect communication of `gpart` over MPI. 

With this in we can run the EAGLE_25 over MPI with gpart and spart activated (e.g. `swift_mpi -s -g -S -t 16`). Note that this does not fix the remaining bug in #256 (Now fixed independently).

See merge request !310

README

@@ -28,6 +28,7 @@ Valid options are:
   -G          Run with self-gravity
   -n    {int} Execute a fixed number of time steps. When unset use the time_end parameter to stop. 
   -s          Run with SPH
+  -S          Run with stars
   -t    {int} The number of threads to use on each MPI rank. Defaults to 1 if not specified.
   -v     [12] Increase the level of verbosity
   	      1: MPI-rank 0 writes

doc/Doxyfile.in

@@ -762,6 +762,7 @@ WARN_LOGFILE           =
 INPUT                  =  @top_srcdir@ @top_srcdir@/src @top_srcdir@/tests @top_srcdir@/examples
 INPUT		       += @top_srcdir@/src/hydro/Minimal
 INPUT		       += @top_srcdir@/src/gravity/Default
+INPUT		       += @top_srcdir@/src/stars/Default
 INPUT		       += @top_srcdir@/src/riemann
 INPUT		       += @top_srcdir@/src/potential/point_mass
 INPUT		       += @top_srcdir@/src/cooling/const_du

examples/MultiTypes/makeIC.py

@@ -36,6 +36,9 @@ eta = 1.2349             # 48 ngbs with cubic spline kernel
 rhoDM = 1.
 Ldm = int(sys.argv[2])  # Number of particles along one axis
 
+massStars = 0.1
+Lstars = int(sys.argv[3])  # Number of particles along one axis
+
 fileName = "multiTypes.hdf5"
 
 #---------------------------------------------------
@@ -46,6 +49,10 @@ internalEnergy = P / ((gamma - 1.)*rhoGas)
 numDM = Ldm**3
 massDM = boxSize**3 * rhoDM / numDM
 
+numStars = Lstars**3
+massStars = massDM * massStars
+
+
 #--------------------------------------------------
 
 #File
@@ -54,9 +61,9 @@ file = h5py.File(fileName, 'w')
 # Header
 grp = file.create_group("/Header")
 grp.attrs["BoxSize"] = boxSize
-grp.attrs["NumPart_Total"] =  [numGas, numDM, 0, 0, 0, 0]
+grp.attrs["NumPart_Total"] =  [numGas, numDM, 0, 0, numStars, 0]
 grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_ThisFile"] = [numGas, numDM, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numGas, numDM, 0, 0, numStars, 0]
 grp.attrs["Time"] = 0.0
 grp.attrs["NumFilesPerSnapshot"] = 1
 grp.attrs["MassTable"] = [0.0, massDM, 0.0, 0.0, 0.0, 0.0]
@@ -142,4 +149,33 @@ coords[:,2] = x[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
 ds = grp.create_dataset('Coordinates', (numDM, 3), 'd')
 ds[()] = coords
 
+
+
+# Star Particle group
+grp = file.create_group("/PartType4")
+
+v  = zeros((numStars, 3))
+ds = grp.create_dataset('Velocities', (numStars, 3), 'f')
+ds[()] = v
+v = zeros(1)
+
+m = full((numStars, 1), massStars)
+ds = grp.create_dataset('Masses', (numStars,1), 'f')
+ds[()] = m
+m = zeros(1)
+
+ids = linspace(0, numStars, numStars, endpoint=False).reshape((numStars,1))
+ds = grp.create_dataset('ParticleIDs', (numStars, 1), 'L')
+ds[()] = ids + Lgas**3 + 1
+x      = ids % Ldm;
+y      = ((ids - x) / Ldm) % Ldm;
+z      = (ids - x - Ldm * y) / Ldm**2;
+coords = zeros((numStars, 3))
+coords[:,0] = z[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
+coords[:,1] = y[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
+coords[:,2] = x[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
+ds = grp.create_dataset('Coordinates', (numStars, 3), 'd')
+ds[()] = coords
+
+
 file.close()

examples/MultiTypes/run.sh

@@ -4,7 +4,7 @@
 if [ ! -e multiTypes.hdf5 ]
 then
     echo "Generating initial conditions for the multitype box example..."
-    python makeIC.py 50 60
+    python makeIC.py 17 24 12
 fi
 
-../swift -s -g -t 16 multiTypes.yml 2>&1 | tee output.log
+../swift -s -g -S -t 1 multiTypes.yml 2>&1 | tee output.log

examples/Stellar_Disk/makeIC.py

+import numpy as np
+import h5py as h5
+

examples/Stellar_Disk/stellar_disk.yml

+IsothermalPotential
\ No newline at end of file

examples/main.c

@@ -83,6 +83,7 @@ void print_help_message() {
          "Execute a fixed number of time steps. When unset use the time_end "
          "parameter to stop.");
   printf("  %2s %8s %s\n", "-s", "", "Run with SPH");
+  printf("  %2s %8s %s\n", "-S", "", "Run with stars");
   printf("  %2s %8s %s\n", "-t", "{int}",
          "The number of threads to use on each MPI rank. Defaults to 1 if not "
          "specified.");
@@ -156,6 +157,7 @@ int main(int argc, char *argv[]) {
   int with_cooling = 0;
   int with_self_gravity = 0;
   int with_hydro = 0;
+  int with_stars = 0;
   int with_fp_exceptions = 0;
   int with_drift_all = 0;
   int verbose = 0;
@@ -165,7 +167,7 @@ int main(int argc, char *argv[]) {
 
   /* Parse the parameters */
   int c;
-  while ((c = getopt(argc, argv, "acCdDef:FgGhn:st:v:y:")) != -1) switch (c) {
+  while ((c = getopt(argc, argv, "acCdDef:FgGhn:sSt:v:y:")) != -1) switch (c) {
       case 'a':
         with_aff = 1;
         break;
@@ -213,6 +215,9 @@ int main(int argc, char *argv[]) {
       case 's':
         with_hydro = 1;
         break;
+      case 'S':
+        with_stars = 1;
+        break;
       case 't':
         if (sscanf(optarg, "%d", &nr_threads) != 1) {
           if (myrank == 0)
@@ -269,6 +274,9 @@ int main(int argc, char *argv[]) {
   /* Genesis 1.1: And then, there was time ! */
   clocks_set_cpufreq(cpufreq);
 
+  /* How vocal are we ? */
+  const int talking = (verbose == 1 && myrank == 0) || (verbose == 2);
+
   if (myrank == 0 && dry_run)
     message(
         "Executing a dry run. No i/o or time integration will be performed.");
@@ -281,7 +289,7 @@ int main(int argc, char *argv[]) {
 
 /* Report host name(s). */
 #ifdef WITH_MPI
-  if (myrank == 0 || verbose > 1) {
+  if (talking) {
     message("Rank %d running on: %s", myrank, hostname());
   }
 #else
@@ -305,14 +313,12 @@ int main(int argc, char *argv[]) {
   if (myrank == 0) {
     message("sizeof(struct part)  is %4zi bytes.", sizeof(struct part));
     message("sizeof(struct xpart) is %4zi bytes.", sizeof(struct xpart));
+    message("sizeof(struct spart) is %4zi bytes.", sizeof(struct spart));
     message("sizeof(struct gpart) is %4zi bytes.", sizeof(struct gpart));
     message("sizeof(struct task)  is %4zi bytes.", sizeof(struct task));
     message("sizeof(struct cell)  is %4zi bytes.", sizeof(struct cell));
   }
 
-  /* How vocal are we ? */
-  const int talking = (verbose == 1 && myrank == 0) || (verbose == 2);
-
   /* Read the parameter file */
   struct swift_params *params = malloc(sizeof(struct swift_params));
   if (params == NULL) error("Error allocating memory for the parameter file.");
@@ -368,26 +374,32 @@ int main(int argc, char *argv[]) {
   if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
   fflush(stdout);
 
+  /* Get ready to read particles of all kinds */
   struct part *parts = NULL;
   struct gpart *gparts = NULL;
-  size_t Ngas = 0, Ngpart = 0;
+  struct spart *sparts = NULL;
+  size_t Ngas = 0, Ngpart = 0, Nspart = 0;
   double dim[3] = {0., 0., 0.};
   int periodic = 0;
   int flag_entropy_ICs = 0;
   if (myrank == 0) clocks_gettime(&tic);
 #if defined(WITH_MPI)
 #if defined(HAVE_PARALLEL_HDF5)
-  read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &Ngas, &Ngpart,
-                   &periodic, &flag_entropy_ICs, myrank, nr_nodes,
-                   MPI_COMM_WORLD, MPI_INFO_NULL, dry_run);
+  read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
+                   &Ngpart, &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
+                   (with_external_gravity || with_self_gravity), with_stars,
+                   myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, dry_run);
 #else
-  read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &Ngas, &Ngpart,
-                 &periodic, &flag_entropy_ICs, myrank, nr_nodes, MPI_COMM_WORLD,
-                 MPI_INFO_NULL, dry_run);
+  read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas, &Ngpart,
+                 &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
+                 (with_external_gravity || with_self_gravity), with_stars,
+                 myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, dry_run);
 #endif
 #else
-  read_ic_single(ICfileName, &us, dim, &parts, &gparts, &Ngas, &Ngpart,
-                 &periodic, &flag_entropy_ICs, dry_run);
+  read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas, &Ngpart,
+                 &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
+                 (with_external_gravity || with_self_gravity), with_stars,
+                 dry_run);
 #endif
   if (myrank == 0) {
     clocks_gettime(&toc);
@@ -396,40 +408,40 @@ int main(int argc, char *argv[]) {
     fflush(stdout);
   }
 
-  /* Discard gparts if we don't have gravity
-   * (Better implementation of i/o will come)*/
-  if (!with_external_gravity && !with_self_gravity) {
-    free(gparts);
-    gparts = NULL;
-    for (size_t k = 0; k < Ngas; ++k) parts[k].gpart = NULL;
-    Ngpart = 0;
+#ifdef SWIFT_DEBUG_CHECKS
+  /* Check once and for all that we don't have unwanted links */
+  if (!with_stars) {
+    for (size_t k = 0; k < Ngpart; ++k)
+      if (gparts[k].type == swift_type_star) error("Linking problem");
   }
   if (!with_hydro) {
-    free(parts);
-    parts = NULL;
     for (size_t k = 0; k < Ngpart; ++k)
-      if (gparts[k].id_or_neg_offset < 0) error("Linking problem");
-    Ngas = 0;
+      if (gparts[k].type == swift_type_gas) error("Linking problem");
   }
+#endif
 
   /* Get the total number of particles across all nodes. */
-  long long N_total[2] = {0, 0};
+  long long N_total[3] = {0, 0, 0};
 #if defined(WITH_MPI)
-  long long N_long[2] = {Ngas, Ngpart};
-  MPI_Reduce(&N_long, &N_total, 2, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD);
+  long long N_long[3] = {Ngas, Ngpart, Nspart};
+  MPI_Reduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM, 0,
+             MPI_COMM_WORLD);
 #else
   N_total[0] = Ngas;
   N_total[1] = Ngpart;
+  N_total[2] = Nspart;
 #endif
   if (myrank == 0)
-    message("Read %lld gas particles and %lld gparts from the ICs.", N_total[0],
-            N_total[1]);
+    message(
+        "Read %lld gas particles, %lld star particles and %lld gparts from the "
+        "ICs.",
+        N_total[0], N_total[2], N_total[1]);
 
   /* Initialize the space with these data. */
   if (myrank == 0) clocks_gettime(&tic);
   struct space s;
-  space_init(&s, params, dim, parts, gparts, Ngas, Ngpart, periodic,
-             with_self_gravity, talking, dry_run);
+  space_init(&s, params, dim, parts, gparts, sparts, Ngas, Ngpart, Nspart,
+             periodic, with_self_gravity, talking, dry_run);
   if (myrank == 0) {
     clocks_gettime(&toc);
     message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
@@ -489,6 +501,7 @@ int main(int argc, char *argv[]) {
   if (with_cosmology) engine_policies |= engine_policy_cosmology;
   if (with_cooling) engine_policies |= engine_policy_cooling;
   if (with_sourceterms) engine_policies |= engine_policy_sourceterms;
+  if (with_stars) engine_policies |= engine_policy_stars;
 
   /* Initialize the engine with the space and policies. */
   if (myrank == 0) clocks_gettime(&tic);
@@ -510,11 +523,16 @@ int main(int argc, char *argv[]) {
 
   /* Get some info to the user. */
   if (myrank == 0) {
+    long long N_DM = N_total[1] - N_total[2] - N_total[0];
+    message(
+        "Running on %lld gas particles, %lld star particles and %lld DM "
+        "particles (%lld gravity particles)",
+        N_total[0], N_total[2], N_total[1] > 0 ? N_DM : 0, N_total[1]);
     message(
-        "Running on %lld gas particles and %lld DM particles from t=%.3e until "
-        "t=%.3e with %d threads and %d queues (dt_min=%.3e, dt_max=%.3e)...",
-        N_total[0], N_total[1], e.timeBegin, e.timeEnd, e.nr_threads,
-        e.sched.nr_queues, e.dt_min, e.dt_max);
+        "from t=%.3e until t=%.3e with %d threads and %d queues (dt_min=%.3e, "
+        "dt_max=%.3e)...",
+        e.timeBegin, e.timeEnd, e.nr_threads, e.sched.nr_queues, e.dt_min,
+        e.dt_max);
     fflush(stdout);
   }
 
@@ -545,8 +563,9 @@ int main(int argc, char *argv[]) {
 
   /* Legend */
   if (myrank == 0)
-    printf("# %6s %14s %14s %10s %10s %16s [%s]\n", "Step", "Time", "Time-step",
-           "Updates", "g-Updates", "Wall-clock time", clocks_getunit());
+    printf("# %6s %14s %14s %10s %10s %10s %16s [%s]\n", "Step", "Time",
+           "Time-step", "Updates", "g-Updates", "s-Updates", "Wall-clock time",
+           clocks_getunit());
 
   /* Main simulation loop */
   for (int j = 0; !engine_is_done(&e) && e.step != nsteps; j++) {

src/Makefile.am

@@ -60,7 +60,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
 nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
 		 kernel_long_gravity.h vector.h cache.h runner_doiact.h runner_doiact_vec.h runner_doiact_grav.h runner_doiact_fft.h \
                  runner_doiact_nosort.h units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
-		 dimension.h equation_of_state.h \
+		 dimension.h equation_of_state.h part_type.h \
 		 gravity.h gravity_io.h \
 		 gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
 		 gravity/Default/gravity_debug.h gravity/Default/gravity_part.h  \
@@ -78,6 +78,9 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h
                  hydro/Gizmo/hydro_debug.h hydro/Gizmo/hydro_part.h \
 	         riemann/riemann_hllc.h riemann/riemann_trrs.h \
 		 riemann/riemann_exact.h riemann/riemann_vacuum.h \
+	 	 stars.h stars_io.h \
+		 stars/Default/star.h stars/Default/star_iact.h stars/Default/star_io.h \
+		 stars/Default/star_debug.h stars/Default/star_part.h  \
 	         potential/none/potential.h potential/point_mass/potential.h \
                  potential/isothermal/potential.h potential/disc_patch/potential.h \
 		 cooling/none/cooling.h cooling/none/cooling_struct.h \

src/active.h

@@ -141,4 +141,28 @@ __attribute__((always_inline)) INLINE static int gpart_is_active(
   return (ti_end == ti_current);
 }
 
+/**
+ * @brief Is this s-particle active ?
+ *
+ * @param sp The #spart.
+ * @param e The #engine containing information about the current time.
+ * @return 1 if the #spart is active, 0 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int spart_is_active(
+    const struct spart *sp, const struct engine *e) {
+
+  const integertime_t ti_current = e->ti_current;
+  const integertime_t ti_end = get_integer_time_end(ti_current, sp->time_bin);
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (ti_end < ti_current)
+    error(
+        "s-particle in an impossible time-zone! gp->ti_end=%lld "
+        "e->ti_current=%lld",
+        ti_end, ti_current);
+#endif
+
+  return (ti_end == ti_current);
+}
+
 #endif /* SWIFT_ACTIVE_H */

src/cache.h

@@ -162,4 +162,22 @@ __attribute__((always_inline)) INLINE void cache_read_particles(
 #endif
 }
 
+/**
+ * @brief Clean the memory allocated by a #cache object.
+ *
+ * @param c The #cache to clean.
+ */
+static INLINE void cache_clean(struct cache *c) {
+  if (c->count > 0) {
+    free(c->x);
+    free(c->y);
+    free(c->z);
+    free(c->m);
+    free(c->vx);
+    free(c->vy);
+    free(c->vz);
+    free(c->h);
+  }
+}
+
 #endif /* SWIFT_CACHE_H */

src/cell.c

@@ -102,6 +102,7 @@ int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {
   c->ti_old = pc->ti_old;
   c->count = pc->count;
   c->gcount = pc->gcount;
+  c->scount = pc->scount;
   c->tag = pc->tag;
 
   /* Number of new cells created. */
@@ -114,6 +115,7 @@ int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {
       space_getcells(s, 1, &temp);
       temp->count = 0;
       temp->gcount = 0;
+      temp->scount = 0;
       temp->loc[0] = c->loc[0];
       temp->loc[1] = c->loc[1];
       temp->loc[2] = c->loc[2];
@@ -194,6 +196,31 @@ int cell_link_gparts(struct cell *c, struct gpart *gparts) {
   return c->gcount;
 }
 
+/**
+ * @brief Link the cells recursively to the given #spart array.
+ *
+ * @param c The #cell.
+ * @param sparts The #spart array.
+ *
+ * @return The number of particles linked.
+ */
+int cell_link_sparts(struct cell *c, struct spart *sparts) {
+
+  c->sparts = sparts;
+
+  /* Fill the progeny recursively, depth-first. */
+  if (c->split) {
+    int offset = 0;
+    for (int k = 0; k < 8; k++) {
+      if (c->progeny[k] != NULL)
+        offset += cell_link_sparts(c->progeny[k], &sparts[offset]);
+    }
+  }
+
+  /* Return the total number of linked particles. */
+  return c->scount;
+}
+
 /**
  * @brief Pack the data of the given cell and all it's sub-cells.
  *
@@ -214,6 +241,7 @@ int cell_pack(struct cell *c, struct pcell *pc) {
   pc->ti_old = c->ti_old;
   pc->count = c->count;
   pc->gcount = c->gcount;
+  pc->scount = c->scount;
   c->tag = pc->tag = atomic_inc(&cell_next_tag) % cell_max_tag;
 
   /* Fill in the progeny, depth-first recursion. */
@@ -425,6 +453,70 @@ int cell_glocktree(struct cell *c) {
   }
 }
 
+/**
+ * @brief Lock a cell for access to its array of #spart and hold its parents.
+ *
+ * @param c The #cell.
+ * @return 0 on success, 1 on failure
+ */
+int cell_slocktree(struct cell *c) {
+
+  TIMER_TIC
+
+  /* First of all, try to lock this cell. */
+  if (c->shold || lock_trylock(&c->slock) != 0) {
+    TIMER_TOC(timer_locktree);
+    return 1;
+  }
+
+  /* Did somebody hold this cell in the meantime? */
+  if (c->shold) {
+
+    /* Unlock this cell. */
+    if (lock_unlock(&c->slock) != 0) error("Failed to unlock cell.");
+
+    /* Admit defeat. */
+    TIMER_TOC(timer_locktree);
+    return 1;
+  }
+
+  /* Climb up the tree and lock/hold/unlock. */
+  struct cell *finger;
+  for (finger = c->parent; finger != NULL; finger = finger->parent) {
+
+    /* Lock this cell. */
+    if (lock_trylock(&finger->slock) != 0) break;
+
+    /* Increment the hold. */
+    atomic_inc(&finger->shold);
+
+    /* Unlock the cell. */
+    if (lock_unlock(&finger->slock) != 0) error("Failed to unlock cell.");
+  }
+
+  /* If we reached the top of the tree, we're done. */
+  if (finger == NULL) {
+    TIMER_TOC(timer_locktree);
+    return 0;
+  }
+
+  /* Otherwise, we hit a snag. */
+  else {
+
+    /* Undo the holds up to finger. */
+    for (struct cell *finger2 = c->parent; finger2 != finger;
+         finger2 = finger2->parent)
+      atomic_dec(&finger2->shold);
+
+    /* Unlock this cell. */
+    if (lock_unlock(&c->slock) != 0) error("Failed to unlock cell.");
+
+    /* Admit defeat. */
+    TIMER_TOC(timer_locktree);
+    return 1;
+  }
+}
+
 /**
  * @brief Unlock a cell's parents for access to #part array.
  *
@@ -463,24 +555,49 @@ void cell_gunlocktree(struct cell *c) {
   TIMER_TOC(timer_locktree);
 }
 
+/**
+ * @brief Unlock a cell's parents for access to #spart array.
+ *
+ * @param c The #cell.
+ */
+void cell_sunlocktree(struct cell *c) {
+
+  TIMER_TIC
+
+  /* First of all, try to unlock this cell. */
+  if (lock_unlock(&c->slock) != 0) error("Failed to unlock cell.");
+
+  /* Climb up the tree and unhold the parents. */
+  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
+    atomic_dec(&finger->shold);
+
+  TIMER_TOC(timer_locktree);
+}
+
 /**
  * @brief Sort the parts into eight bins along the given pivots.
  *
  * @param c The #cell array to be sorted.
  * @param parts_offset Offset of the cell parts array relative to the
  *        space's parts array, i.e. c->parts - s->parts.
+ * @param sparts_offset Offset of the cell sparts array relative to the
+ *        space's sparts array, i.e. c->sparts - s->sparts.
  * @param buff A buffer with at least max(c->count, c->gcount) entries,
  *        used for sorting indices.
+ * @param sbuff A buffer with at least max(c->scount, c->gcount) entries,
+ *        used for sorting indices for the sparts.
  * @param gbuff A buffer with at least max(c->count, c->gcount) entries,
  *        used for sorting indices for the gparts.
  */
-void cell_split(struct cell *c, ptrdiff_t parts_offset, struct cell_buff *buff,
+void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
+                struct cell_buff *buff, struct cell_buff *sbuff,
                 struct cell_buff *gbuff) {
 
-  const int count = c->count, gcount = c->gcount;
+  const int count = c->count, gcount = c->gcount, scount = c->scount;
   struct part *parts = c->parts;
   struct xpart *xparts = c->xparts;
   struct gpart *gparts = c->gparts;
+  struct spart *sparts = c->sparts;
   const double pivot[3] = {c->loc[0] + c->width[0] / 2,
                            c->loc[1] + c->width[1] / 2,
                            c->loc[2] + c->width[2] / 2};
@@ -494,6 +611,16 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, struct cell_buff *buff,
         buff[k].x[2] != parts[k].x[2])
       error("Inconsistent buff contents.");
   }
+  for (int k = 0; k < gcount; k++) {
+    if (gbuff[k].x[0] != gparts[k].x[0] || gbuff[k].x[1] != gparts[k].x[1] ||
+        gbuff[k].x[2] != gparts[k].x[2])
+      error("Inconsistent gbuff contents.");
+  }
+  for (int k = 0; k < scount; k++) {
+    if (sbuff[k].x[0] != sparts[k].x[0] || sbuff[k].x[1] != sparts[k].x[1] ||
+        sbuff[k].x[2] != sparts[k].x[2])
+      error("Inconsistent sbuff contents.");
+  }
 #endif /* SWIFT_DEBUG_CHECKS */
 
   /* Fill the buffer with the indices. */
@@ -547,7 +674,8 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, struct cell_buff *buff,
   }
 
   /* Re-link the gparts. */
-  if (count > 0 && gcount > 0) part_relink_gparts(parts, count, parts_offset);
+  if (count > 0 && gcount > 0)
+    part_relink_gparts_to_parts(parts, count, parts_offset);
 
 #ifdef SWIFT_DEBUG_CHECKS
   /* Check that the buffs are OK. */
@@ -611,7 +739,60 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, struct cell_buff *buff,
       error("Sorting failed (progeny=7).");
 #endif
 
-  /* Now do the same song and dance for the gparts. */
+  /* Now do the same song and dance for the sparts. */
+  for (int k = 0; k < 8; k++) bucket_count[k] = 0;
+
+  /* Fill the buffer with the indices. */
+  for (int k = 0; k < scount; k++) {
+    const int bid = (sbuff[k].x[0] > pivot[0]) * 4 +
+                    (sbuff[k].x[1] > pivot[1]) * 2 + (sbuff[k].x[2] > pivot[2]);
+    bucket_count[bid]++;
+    sbuff[k].ind = bid;
+  }
+
+  /* Set the buffer offsets. */
+  bucket_offset[0] = 0;
+  for (int k = 1; k <= 8; k++) {
+    bucket_offset[k] = bucket_offset[k - 1] + bucket_count[k - 1];
+    bucket_count[k - 1] = 0;
+  }
+
+  /* Run through the buckets, and swap particles to their correct spot. */
+  for (int bucket = 0; bucket < 8; bucket++) {
+    for (int k = bucket_offset[bucket] + bucket_count[bucket];
+         k < bucket_offset[bucket + 1]; k++) {
+      int bid = sbuff[k].ind;
+      if (bid != bucket) {
+        struct spart spart = sparts[k];
+