Merge branch 'remap_ids_in_ics' into 'master'

Cleverer ids for particle splitting

See merge request !1106

README

@@ -37,6 +37,7 @@ Parameters:
     -s, --hydro                       Run with hydrodynamics.
     -S, --stars                       Run with stars.
     -B, --black-holes                 Run with black holes.
+    -k, --sinks                       Run with sink particles. 
     -u, --fof                         Run Friends-of-Friends algorithm and 
                                       black holes seeding.
     -x, --velociraptor                Run with structure finding.

README.md

@@ -124,6 +124,7 @@ Parameters:
     -s, --hydro                       Run with hydrodynamics.
     -S, --stars                       Run with stars.
     -B, --black-holes                 Run with black holes.
+    -k, --sinks                       Run with sink particles.
     -u, --fof                         Run Friends-of-Friends algorithm to
                                       perform black hole seeding.
     -x, --velociraptor                Run with structure finding.

doc/RTD/source/CommandLineOptions/index.rst

@@ -34,6 +34,7 @@ can be found by typing ``./swift -h``:
     -s, --hydro                       Run with hydrodynamics.
     -S, --stars                       Run with stars.
     -B, --black-holes                 Run with black holes.
+    -k, --sinks                       Run with sink particles. 
     -u, --fof                         Run Friends-of-Friends algorithm to
                                       perform black hole seeding.
     -x, --velociraptor                Run with structure finding.

doc/RTD/source/ParameterFiles/parameter_description.rst

@@ -389,7 +389,13 @@ with a mass above a fixed threshold into two copies that are slightly shifted
 (by a randomly orientated vector of norm :math:`0.2h`). Their masses and other
 relevant particle-carried quantities are then halved. The mass threshold for
 splitting is set by the parameter ``particle_splitting_mass_threshold`` which is
-specified using the internal unit system.
+specified using the internal unit system. The IDs of the newly created particles
+can be either drawn randomly by setting the parameter ``generate_random_ids``
+(Default: 0) to :math:`1`. When this is activated, there is no check that the
+newly generated IDs do not clash with any other pre-existing particle. If this
+option is set to :math:`0` (the default setting) then the new IDs are created in
+increasing order from the maximal pre-existing value in the simulation, hence
+preventing any clash.
 
 The final set of parameters in this section determine the initial and minimum
 temperatures of the particles.
@@ -578,6 +584,8 @@ Finally, SWIFT also offers these options:
 * A factor to re-scale all the smoothing-lengths by a fixed amount: ``smoothing_length_scaling`` (default: ``1.``),
 * A shift to apply to all the particles: ``shift`` (default: ``[0.0,0.0,0.0]``),
 * Whether to replicate the box along each axis: ``replicate`` (default: ``1``).
+* Whether to re-map the IDs to the range ``[0, N]`` and hence discard
+  the original IDs from the IC file: ``remap_ids`` (default: ``0``).
   
 The shift is expressed in internal units. The option to replicate the
 box is especially useful for weak-scaling tests. When set to an
@@ -585,6 +593,11 @@ integer >1, the box size is multiplied by this integer along each axis
 and the particles are duplicated and shifted such as to create exact
 copies of the simulation volume.
 
+The remapping of IDs is especially useful in combination with the option to
+generate increasing IDs when splitting gas particles as it allows for the
+creation of a compact range of IDs beyond which the new IDs generated by
+splitting can be safely drawn from.
+
 The full section to start a DM+hydro run from Gadget DM-only ICs would
 be:
 

examples/EAGLE_ICs/EAGLE_100/eagle_100.yml

@@ -88,6 +88,7 @@ InitialConditions:
   cleanup_velocity_factors: 1        # Remove the sqrt(a) factor in the velocities inherited from Gadget
   generate_gas_in_ics: 1             # Generate gas particles from the DM-only ICs
   cleanup_smoothing_lengths: 1       # Since we generate gas, make use of the (expensive) cleaning-up procedure.
+  remap_ids: 1                       # Re-map the IDs to [1, N] to avoid collision problems when splitting
 
 # Impose primoridal metallicity
 EAGLEChemistry:

examples/EAGLE_ICs/EAGLE_12/eagle_12.yml

@@ -91,6 +91,7 @@ InitialConditions:
   cleanup_velocity_factors: 1        # Remove the sqrt(a) factor in the velocities inherited from Gadget
   generate_gas_in_ics: 1             # Generate gas particles from the DM-only ICs
   cleanup_smoothing_lengths: 1       # Since we generate gas, make use of the (expensive) cleaning-up procedure.
+  remap_ids: 1                       # Re-map the IDs to [1, N] to avoid collision problems when splitting
 
 # Parameters of the line-of-sight outputs
 LineOfSight:

examples/EAGLE_ICs/EAGLE_25/eagle_25.yml

@@ -91,6 +91,7 @@ InitialConditions:
   cleanup_velocity_factors: 1        # Remove the sqrt(a) factor in the velocities inherited from Gadget
   generate_gas_in_ics: 1             # Generate gas particles from the DM-only ICs
   cleanup_smoothing_lengths: 1       # Since we generate gas, make use of the (expensive) cleaning-up procedure.
+  remap_ids: 1                       # Re-map the IDs to [1, N] to avoid collision problems when splitting
   
 # Parameters of the line-of-sight outputs
 LineOfSight:

examples/EAGLE_ICs/EAGLE_50/eagle_50.yml

@@ -89,6 +89,7 @@ InitialConditions:
   cleanup_velocity_factors: 1        # Remove the sqrt(a) factor in the velocities inherited from Gadget
   generate_gas_in_ics: 1             # Generate gas particles from the DM-only ICs
   cleanup_smoothing_lengths: 1       # Since we generate gas, make use of the (expensive) cleaning-up procedure.
+  remap_ids: 1                       # Re-map the IDs to [1, N] to avoid collision problems when splitting
 
 # Parameters of the line-of-sight outputs
 LineOfSight:

examples/EAGLE_ICs/EAGLE_50_low_res/eagle_50.yml

@@ -88,6 +88,7 @@ InitialConditions:
   cleanup_velocity_factors: 1        # Remove the sqrt(a) factor in the velocities inherited from Gadget
   generate_gas_in_ics: 1             # Generate gas particles from the DM-only ICs
   cleanup_smoothing_lengths: 1       # Since we generate gas, make use of the (expensive) cleaning-up procedure.
+  remap_ids: 1                       # Re-map the IDs to [1, N] to avoid collision problems when splitting
 
 # Parameters of the line-of-sight outputs
 LineOfSight:

examples/EAGLE_low_z/EAGLE_6/eagle_6.yml

@@ -96,6 +96,7 @@ InitialConditions:
   periodic:   1
   cleanup_h_factors: 1               # Remove the h-factors inherited from Gadget
   cleanup_velocity_factors: 1        # Remove the sqrt(a) factor in the velocities inherited from Gadget
+  remap_ids: 1
   
 EAGLEChemistry: 	     # Solar abundances
   init_abundance_metal:      0.014

examples/main.c

@@ -903,6 +903,8 @@ int main(int argc, char *argv[]) {
         params, "InitialConditions:cleanup_velocity_factors", 0);
     const int generate_gas_in_ics = parser_get_opt_param_int(
         params, "InitialConditions:generate_gas_in_ics", 0);
+    const int remap_ids =
+        parser_get_opt_param_int(params, "InitialConditions:remap_ids", 0);
 
     /* Initialise the cosmology */
     if (with_cosmology)
@@ -1140,9 +1142,9 @@ int main(int argc, char *argv[]) {
     if (myrank == 0) clocks_gettime(&tic);
     space_init(&s, params, &cosmo, dim, &hydro_properties, parts, gparts, sinks,
                sparts, bparts, Ngas, Ngpart, Nsink, Nspart, Nbpart, periodic,
-               replicate, generate_gas_in_ics, with_hydro, with_self_gravity,
-               with_star_formation, with_DM_background_particles, talking,
-               dry_run, nr_nodes);
+               replicate, remap_ids, generate_gas_in_ics, with_hydro,
+               with_self_gravity, with_star_formation,
+               with_DM_background_particles, talking, dry_run, nr_nodes);
 
     /* Initialise the line of sight properties. */
     if (with_line_of_sight) los_init(s.dim, &los_properties, params);

examples/main_fof.c

@@ -546,7 +546,7 @@ int main(int argc, char *argv[]) {
   if (myrank == 0) clocks_gettime(&tic);
   space_init(&s, params, &cosmo, dim, /*hydro_props=*/NULL, parts, gparts,
              sinks, sparts, bparts, Ngas, Ngpart, Nsink, Nspart, Nbpart,
-             periodic, replicate,
+             periodic, replicate, /*remap_ids=*/0,
              /*generate_gas_in_ics=*/0, /*hydro=*/N_total[0] > 0, /*gravity=*/1,
              /*with_star_formation=*/0, with_DM_background_particles, talking,
              /*dry_run=*/0, nr_nodes);

examples/parameter_example.yml

@@ -38,6 +38,7 @@ SPH:
   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.
@@ -178,6 +179,7 @@ InitialConditions:
   smoothing_length_scaling:    1.   # (Optional) A scaling factor to apply to all smoothing lengths in the ICs.
   shift:      [0.0,0.0,0.0]         # (Optional) A shift to apply to all particles read from the ICs (in internal units).
   replicate:  2                     # (Optional) Replicate all particles along each axis a given integer number of times. Default 1.
+  remap_ids:  0                     # (Optional) Remap all the particle IDs to the range [1, NumPart].
 
 # Parameters controlling restarts
 Restarts:

src/atomic.h

@@ -168,6 +168,27 @@ __attribute__((always_inline)) INLINE static void atomic_max(
   } while (test_val != old_val);
 }
 
+/**
+ * @brief Atomic max operation on long long.
+ *
+ * This is a text-book implementation based on an atomic CAS.
+ *
+ * @param address The address to update.
+ * @param y The value to update the address with.
+ */
+__attribute__((always_inline)) INLINE static void atomic_max_ll(
+    volatile long long *const address, const long long y) {
+
+  int test_val, old_val, new_val;
+  old_val = *address;
+
+  do {
+    test_val = old_val;
+    new_val = max(old_val, y);
+    old_val = atomic_cas(address, test_val, new_val);
+  } while (test_val != old_val);
+}
+
 /**
  * @brief Atomic max operation on floats.
  *

src/engine.c

@@ -2170,6 +2170,21 @@ void engine_first_init_particles(struct engine *e) {
             clocks_getunit());
 }
 
+/**
+ * @brief Compute the maximal ID of any #part in the run.
+ *
+ * @param e The #engine.
+ */
+void engine_get_max_ids(struct engine *e) {
+
+  e->max_parts_id = space_get_max_parts_id(e->s);
+
+#ifdef WITH_MPI
+  MPI_Allreduce(MPI_IN_PLACE, &e->max_parts_id, 1, MPI_LONG_LONG_INT, MPI_MAX,
+                MPI_COMM_WORLD);
+#endif
+}
+
 /**
  * @brief Initialises the particles and set them in a state ready to move
  *forward in time.
@@ -2431,6 +2446,9 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
                     e->verbose);
 #endif
 
+  /* Gather the max IDs at this stage */
+  engine_get_max_ids(e);
+
   /* Ready to go */
   e->step = 0;
   e->forcerebuild = 1;

src/engine.h

@@ -275,6 +275,10 @@ struct engine {
   long long count_inhibited_bparts;
 #endif
 
+  /* Maximal ID of the parts (used for the generation of new IDs when splitting)
+   */
+  long long max_parts_id;
+
   /* Total mass in the simulation */
   double total_mass;
 

src/engine_split_particles.c

@@ -43,6 +43,7 @@ struct data_count {
   const struct engine *const e;
   const float mass_threshold;
   size_t counter;
+  long long max_id;
 };
 
 /**
@@ -51,8 +52,11 @@ struct data_count {
 struct data_split {
   const struct engine *const e;
   const float mass_threshold;
+  const int generate_random_ids;
   size_t *const k_parts;
   size_t *const k_gparts;
+  long long offset_id;
+  long long *count_id;
   swift_lock_type lock;
 };
 
@@ -70,6 +74,7 @@ void engine_split_gas_particle_count_mapper(void *restrict map_data, int count,
   const float mass_threshold = data->mass_threshold;
 
   size_t counter = 0;
+  long long max_id = 0;
 
   for (int i = 0; i < count; ++i) {
 
@@ -81,10 +86,14 @@ void engine_split_gas_particle_count_mapper(void *restrict map_data, int count,
     /* Is the mass of this particle larger than the threshold? */
     const float gas_mass = hydro_get_mass(p);
     if (gas_mass > mass_threshold) ++counter;
+
+    /* Get the maximal id */
+    max_id = max(max_id, p->id);
   }
 
   /* Increment the global counter */
   atomic_add(&data->counter, counter);
+  atomic_max_ll(&data->max_id, max_id);
 }
 
 /**
@@ -97,6 +106,9 @@ void engine_split_gas_particle_split_mapper(void *restrict map_data, int count,
   struct part *parts = (struct part *)map_data;
   struct data_split *data = (struct data_split *)extra_data;
   const float mass_threshold = data->mass_threshold;
+  const int generate_random_ids = data->generate_random_ids;
+  const long long offset_id = data->offset_id;
+  long long *count_id = (long long *)data->count_id;
   const struct engine *e = data->e;
   const int with_gravity = (e->policy & engine_policy_self_gravity) ||
                            (e->policy & engine_policy_external_gravity);
@@ -158,10 +170,14 @@ void engine_split_gas_particle_split_mapper(void *restrict map_data, int count,
         memcpy(&global_gparts[k_gparts], gp, sizeof(struct gpart));
       }
 
-      /* Update the IDs.
-       * The gas IDs are always odd, so we multiply by two here to
+      /* Update the IDs. */
+      if (generate_random_ids) {
+        /* The gas IDs are always odd, so we multiply by two here to
          * repsect the parity. */
         global_parts[k_parts].id += 2 * (long long)rand_r(&seedp);
+      } else {
+        global_parts[k_parts].id = offset_id + 2 * atomic_inc(count_id);
+      }
 
       /* Re-link everything */
       if (with_gravity) {
@@ -264,6 +280,7 @@ void engine_split_gas_particles(struct engine *e) {
                            (e->policy & engine_policy_external_gravity);
   const float mass_threshold =
       e->hydro_properties->particle_splitting_mass_threshold;
+  const int generate_random_ids = e->hydro_properties->generate_random_ids;
   const size_t nr_parts_old = s->nr_parts;
 
   /* Quick check to avoid problems */
@@ -274,12 +291,43 @@ void engine_split_gas_particles(struct engine *e) {
 
   /* Start by counting how many particles are above the threshold
    * for splitting (this is done in parallel over the threads) */
-  struct data_count data_count = {e, mass_threshold, 0};
+  struct data_count data_count = {e, mass_threshold, 0, 0};
   threadpool_map(&e->threadpool, engine_split_gas_particle_count_mapper,
                  s->parts, nr_parts_old, sizeof(struct part), 0, &data_count);
   const size_t counter = data_count.counter;
 
-  /* Early abort? */
+  /* Verify that nothing wrong happened with the IDs */
+  if (data_count.max_id > e->max_parts_id) {
+    error("Found a gas particle with an ID larger than the current max!");
+  }
+
+  /* Be verbose about this. This is an important event */
+  message("Splitting %zd particles above the mass threshold", counter);
+
+  /* Number of particles to create */
+  const long long count_new_gas = counter * particle_split_factor;
+
+  /* Get the global offset to generate new IDs (the *2 is to respect the ID
+   * parity) */
+  long long expected_count_id = 2 * counter * (particle_split_factor - 1);
+  long long offset_id = 0;
+#ifdef WITH_MPI
+  MPI_Exscan(&expected_count_id, &offset_id, 1, MPI_LONG_LONG_INT, MPI_SUM,
+             MPI_COMM_WORLD);
+#endif
+  offset_id += e->max_parts_id + 1;
+
+  /* Store the new maximal id */
+  e->max_parts_id = offset_id + expected_count_id;
+#ifdef WITH_MPI
+  MPI_Bcast(&e->max_parts_id, 1, MPI_LONG_LONG_INT, e->nr_nodes - 1,
+            MPI_COMM_WORLD);
+#endif
+
+  /* Each node now has a unique range of IDs [offset_id, offset_id + count_id]
+   */
+
+  /* Early abort (i.e. no particle to split on this MPI rank) ? */
   if (counter == 0) {
 
     if (e->verbose)
@@ -288,12 +336,6 @@ void engine_split_gas_particles(struct engine *e) {
     return;
   }
 
-  /* Be verbose about this. This is an important event */
-  message("Splitting %zd particles above the mass threshold", counter);
-
-  /* Number of particles to create */
-  const size_t count_new_gas = counter * particle_split_factor;
-
   /* Do we need to reallocate the gas array for the new particles? */
   if (s->nr_parts + count_new_gas > s->size_parts) {
 
@@ -376,12 +418,22 @@ void engine_split_gas_particles(struct engine *e) {
   size_t k_gparts = s->nr_gparts;
 
   /* Loop over the particles again to split them */
-  struct data_split data_split = {e, mass_threshold, &k_parts, &k_gparts, 0};
+  long long local_count_id = 0;
+  struct data_split data_split = {
+      e,         mass_threshold, generate_random_ids, &k_parts,
+      &k_gparts, offset_id,      &local_count_id,     0};
   lock_init(&data_split.lock);
   threadpool_map(&e->threadpool, engine_split_gas_particle_split_mapper,
                  s->parts, nr_parts_old, sizeof(struct part), 0, &data_split);
   if (lock_destroy(&data_split.lock) != 0) error("Error destroying lock");
 
+  /* Check that ID assignment went well */
+  if (!generate_random_ids && 2 * local_count_id != expected_count_id)
+    error(
+        "Something went wrong when assigning new IDs expected count=%lld "
+        "actual count=%lld",
+        expected_count_id, local_count_id);
+
   /* Update the local counters */
   s->nr_parts = k_parts;
   s->nr_gparts = k_gparts;

src/hydro_properties.c

@@ -203,6 +203,9 @@ void hydro_props_init(struct hydro_props *p,
   if (p->particle_splitting) {
     p->particle_splitting_mass_threshold =
         parser_get_param_float(params, "SPH:particle_splitting_mass_threshold");
+
+    p->generate_random_ids = parser_get_opt_param_int(
+        params, "SPH:particle_splitting_generate_random_ids", 0);
   }
 }
 

src/hydro_properties.h

@@ -120,6 +120,9 @@ struct hydro_props {
   /*! Mass above which particles get split (internal units) */
   float particle_splitting_mass_threshold;
 
+  /*! Are we generating random IDs when splitting particles? */
+  int generate_random_ids;
+
   /* ------ Viscosity and diffusion ---------------- */
 
   /*! Artificial viscosity parameters */

src/space.c

@@ -5393,6 +5393,7 @@ void space_convert_quantities(struct space *s, int verbose) {
  * @param Nbpart The number of black hole particles in the space.
  * @param periodic flag whether the domain is periodic or not.
  * @param replicate How many replications along each direction do we want?
+ * @param remap_ids Are we remapping the IDs from 1 to N?
  * @param generate_gas_in_ics Are we generating gas particles from the gparts?
  * @param hydro flag whether we are doing hydro or not?
  * @param self_gravity flag whether we are doing gravity or not?
@@ -5413,9 +5414,9 @@ void space_init(struct space *s, struct swift_params *params,
                 struct gpart *gparts, struct sink *sinks, struct spart *sparts,
                 struct bpart *bparts, size_t Npart, size_t Ngpart, size_t Nsink,
                 size_t Nspart, size_t Nbpart, int periodic, int replicate,
-                int generate_gas_in_ics, int hydro, int self_gravity,
-                int star_formation, int DM_background, int verbose, int dry_run,
-                int nr_nodes) {
+                int remap_ids, int generate_gas_in_ics, int hydro,
+                int self_gravity, int star_formation, int DM_background,
+                int verbose, int dry_run, int nr_nodes) {
 
   /* Clean-up everything */
   bzero(s, sizeof(struct space));
@@ -5467,6 +5468,11 @@ void space_init(struct space *s, struct swift_params *params,
   /* Initiate some basic randomness */
   srand(42);
 
+  /* Are we remapping the IDs to the range [1, NumPart]? */
+  if (remap_ids) {
+    space_remap_ids(s, nr_nodes, verbose);
+  }
+
   /* Are we generating gas from the DM-only ICs? */
   if (generate_gas_in_ics) {
     space_generate_gas(s, cosmo, hydro_properties, periodic, DM_background, dim,
@@ -5926,6 +5932,94 @@ void space_replicate(struct space *s, int replicate, int verbose) {
 #endif
 }
 
+/**
+ * @brief Remaps the IDs of the particles to the range [1, N]
+ *
+ * The IDs are unique accross all MPI ranks and are generated
+ * in ther order DM, gas, sinks, stars, BHs.
+ *
+ * @param s The current #space object.
+ * @param nr_nodes The number of MPI ranks used in the run.
+ * @param verbose Are we talkative?
+ */
+void space_remap_ids(struct space *s, int nr_nodes, int verbose) {
+
+  if (verbose) message("Remapping all the IDs");
+
+  /* Get the current local number of particles */
+  const size_t local_nr_parts = s->nr_parts;
+  const size_t local_nr_sinks = s->nr_sinks;
+  const size_t local_nr_gparts = s->nr_gparts;
+  const size_t local_nr_sparts = s->nr_sparts;
+  const size_t local_nr_bparts = s->nr_bparts;
+  const size_t local_nr_baryons =
+      local_nr_parts + local_nr_sinks + local_nr_sparts + local_nr_bparts;
+  const size_t local_nr_dm = local_nr_gparts - local_nr_baryons;
+
+  /* Get the global offsets */
+  long long offset_parts = 0;
+  long long offset_sinks = 0;
+  long long offset_sparts = 0;
+  long long offset_bparts = 0;
+  long long offset_dm = 0;
+#ifdef WITH_MPI
+  MPI_Exscan(&local_nr_parts, &offset_parts, 1, MPI_LONG_LONG_INT, MPI_SUM,
+             MPI_COMM_WORLD);
+  MPI_Exscan(&local_nr_sinks, &offset_sinks, 1, MPI_LONG_LONG_INT, MPI_SUM,
+             MPI_COMM_WORLD);
+  MPI_Exscan(&local_nr_sparts, &offset_sparts, 1, MPI_LONG_LONG_INT, MPI_SUM,
+             MPI_COMM_WORLD);
+  MPI_Exscan(&local_nr_bparts, &offset_bparts, 1, MPI_LONG_LONG_INT, MPI_SUM,
+             MPI_COMM_WORLD);
+  MPI_Exscan(&local_nr_dm, &offset_dm, 1, MPI_LONG_LONG_INT, MPI_SUM,
+             MPI_COMM_WORLD);
+#endif
+
+  /* Total number of particles of each kind */
+  long long total_dm = offset_dm + local_nr_dm;
+  long long total_parts = offset_parts + local_nr_parts;
+  long long total_sinks = offset_sinks + local_nr_sinks;
+  long long total_sparts = offset_sparts + local_nr_sparts;
+  // long long total_bparts = offset_bparts + local_nr_bparts;
+
+#ifdef WITH_MPI
+  /* The last rank now has the correct total, let's broadcast this back */
+  MPI_Bcast(&total_dm, 1, MPI_LONG_LONG_INT, nr_nodes - 1, MPI_COMM_WORLD);
+  MPI_Bcast(&total_parts, 1, MPI_LONG_LONG_INT, nr_nodes - 1, MPI_COMM_WORLD);
+  MPI_Bcast(&total_sinks, 1, MPI_LONG_LONG_INT, nr_nodes - 1, MPI_COMM_WORLD);
+  MPI_Bcast(&total_sparts, 1, MPI_LONG_LONG_INT, nr_nodes - 1, MPI_COMM_WORLD);
+  // MPI_Bcast(&total_bparts, 1, MPI_LONG_LONG_INT, nr_nodes - 1,
+  // MPI_COMM_WORLD);
+#endif
+
+  /* Let's order the particles
+   * IDs will be DM then gas then sinks than stars then BHs */
+  offset_dm += 1;
+  offset_parts += 1 + total_dm;
+  offset_sinks += 1 + total_dm + total_parts;
+  offset_sparts += 1 + total_dm + total_parts + total_sinks;
+  offset_bparts += 1 + total_dm + total_parts + total_sinks + total_sparts;
+
+  /* We can now remap the IDs in the range [offset offset + local_nr] */
+  for (size_t i = 0; i < local_nr_parts; ++i) {
+    s->parts[i].id = offset_parts + i;
+  }
+  for (size_t i = 0; i < local_nr_sinks; ++i) {
+    s->sinks[i].id = offset_sinks + i;
+  }
+  for (size_t i = 0; i < local_nr_sparts; ++i) {
+    s->sparts[i].id = offset_sparts + i;
+  }
+  for (size_t i = 0; i < local_nr_bparts; ++i) {
+    s->bparts[i].id = offset_bparts + i;
+  }
+  for (size_t i = 0; i < local_nr_dm; ++i) {
+    if (s->gparts[i].type == swift_type_dark_matter ||
+        s->gparts[i].type == swift_type_dark_matter_background)
+      s->gparts[i].id_or_neg_offset = offset_dm + i;
+  }
+}
+
 /**
  * @brief Duplicate all the dark matter particles to create the same number
  * of gas particles with mass ratios given by the cosmology.
@@ -6168,6 +6262,29 @@ void space_check_cosmology(struct space *s, const struct cosmology *cosmo,
   }
 }
 
+/**
+ * @brief Compute the max id of any #part in this space.
+ *
+ * This function is inefficient. Don't call often.
+ *
+ * @param s The #space.
+ */
+long long space_get_max_parts_id(struct space *s) {
+
+  long long max_id = -1;
+  for (size_t i = 0; i < s->nr_parts; ++i) max_id = max(max_id, s->parts[i].id);
+  for (size_t i = 0; i < s->nr_sinks; ++i) max_id = max(max_id, s->sinks[i].id);
+  for (size_t i = 0; i < s->nr_sparts; ++i)
+    max_id = max(max_id, s->sparts[i].id);
+  for (size_t i = 0; i < s->nr_bparts; ++i)
+    max_id = max(max_id, s->bparts[i].id);
+  for (size_t i = 0; i < s->nr_gparts; ++i)
+    if (s->gparts[i].type == swift_type_dark_matter ||
+        s->gparts[i].type == swift_type_dark_matter_background)
+      max_id = max(max_id, s->gparts[i].id_or_neg_offset);
+  return max_id;
+}
+
 /**
  * @brief Cleans-up all the cell links in the space
  *