From da8b4cc2d3f43e02b959e4a03550093e92d6446f Mon Sep 17 00:00:00 2001
From: Stuart McAlpine <stuart.mcalpine@helsinki.fi>
Date: Mon, 27 Jul 2020 15:59:42 +0100
Subject: [PATCH] Don't load ParticleIDs when using remap_ids

---
 doc/RTD/source/InitialConditions/index.rst        |  2 ++
 .../ParameterFiles/parameter_description.rst      |  3 ++-
 examples/main.c                                   |  8 ++++----
 examples/main_fof.c                               |  6 +++---
 src/common_io.c                                   | 15 +++++++++++++++
 src/common_io.h                                   |  1 +
 src/parallel_io.c                                 | 12 ++++++++++--
 src/parallel_io.h                                 |  2 +-
 src/serial_io.c                                   | 13 +++++++++++--
 src/serial_io.h                                   |  2 +-
 src/single_io.c                                   | 12 ++++++++++--
 src/single_io.h                                   |  2 +-
 src/space.c                                       |  3 ++-
 13 files changed, 63 insertions(+), 18 deletions(-)

diff --git a/doc/RTD/source/InitialConditions/index.rst b/doc/RTD/source/InitialConditions/index.rst
index 66fb9f01f1..c22f38ad67 100644
--- a/doc/RTD/source/InitialConditions/index.rst
+++ b/doc/RTD/source/InitialConditions/index.rst
@@ -152,6 +152,8 @@ individual particle type (e.g. ``/PartType0/``) that have the following *dataset
   each particle. Note that these have to be unique to a particle, and cannot be
   the same even between particle types. The **IDs must be >= 0**. Negative
   IDs will be rejected by the code.
+  Note, however, that if the parameters to remap the IDs upon startup is switched
+  on (see :ref:`Parameters_ICs`), the IDs can be omitted entirely from the ICs.
 + ``Masses``, an array of length N that gives the masses of the particles.
 
 For ``PartType0`` (i.e. particles that interact through hydro-dynamics), you will
diff --git a/doc/RTD/source/ParameterFiles/parameter_description.rst b/doc/RTD/source/ParameterFiles/parameter_description.rst
index 87dcececee..568047e48b 100644
--- a/doc/RTD/source/ParameterFiles/parameter_description.rst
+++ b/doc/RTD/source/ParameterFiles/parameter_description.rst
@@ -596,7 +596,8 @@ 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.
+splitting can be safely drawn from. Note that, when ``remap_ids`` is
+switched on the ICs do not need to contain a ``ParticleIDs`` field.
 
 The full section to start a DM+hydro run from Gadget DM-only ICs would
 be:
diff --git a/examples/main.c b/examples/main.c
index cb9b0b273c..f1ce38458e 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -1032,15 +1032,15 @@ int main(int argc, char *argv[]) {
                      with_gravity, with_sink, with_stars, with_black_holes,
                      with_cosmology, cleanup_h, cleanup_sqrt_a, cosmo.h,
                      cosmo.a, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL,
-                     nr_threads, dry_run);
+                     nr_threads, dry_run, remap_ids);
 #else
     read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sinks, &sparts,
                    &bparts, &Ngas, &Ngpart, &Ngpart_background, &Nsink, &Nspart,
                    &Nbpart, &flag_entropy_ICs, with_hydro, with_gravity,
                    with_sink, with_stars, with_black_holes, with_cosmology,
                    cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank,
-                   nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
-                   dry_run);
+                   nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads, dry_run,
+                   remap_ids);
 #endif
 #else
     read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sinks, &sparts,
@@ -1048,7 +1048,7 @@ int main(int argc, char *argv[]) {
                    &Nbpart, &flag_entropy_ICs, with_hydro, with_gravity,
                    with_sink, with_stars, with_black_holes, with_cosmology,
                    cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads,
-                   dry_run);
+                   dry_run, remap_ids);
 #endif
 #endif
 
diff --git a/examples/main_fof.c b/examples/main_fof.c
index 2ff8f260e2..9d169ad778 100644
--- a/examples/main_fof.c
+++ b/examples/main_fof.c
@@ -466,7 +466,7 @@ int main(int argc, char *argv[]) {
                    /*with_grav=*/1, with_sinks, with_stars, with_black_holes,
                    with_cosmology, cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a,
                    myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
-                   /*dry_run=*/0);
+                   /*dry_run=*/0, /*remap_ids=*/0);
 #else
   read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sinks, &sparts,
                  &bparts, &Ngas, &Ngpart, &Ngpart_background, &Nsink, &Nspart,
@@ -474,7 +474,7 @@ int main(int argc, char *argv[]) {
                  /*with_grav=*/1, with_sinks, with_stars, with_black_holes,
                  with_cosmology, cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a,
                  myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
-                 /*dry_run=*/0);
+                 /*dry_run=*/0, /*remap_ids=*/0);
 #endif
 #else
   read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sinks, &sparts,
@@ -482,7 +482,7 @@ int main(int argc, char *argv[]) {
                  &Nbpart, &flag_entropy_ICs, with_hydro,
                  /*with_grav=*/1, with_sinks, with_stars, with_black_holes,
                  with_cosmology, cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a,
-                 nr_threads, /*dry_run=*/0);
+                 nr_threads, /*dry_run=*/0, /*remap_ids=*/0);
 #endif
 #endif
   if (myrank == 0) {
diff --git a/src/common_io.c b/src/common_io.c
index ef66a9ed4d..b44ce1bd09 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -2930,3 +2930,18 @@ int get_param_ptype(const char* name) {
    * an int after promising to do so... */
   return -1;
 }
+
+/**
+ * @brief Set all ParticleIDs for each gpart to 1.
+ *
+ * Function is called when remap_ids is 1.
+ *
+ * Note only the gparts IDs have to be set to 1, as other parttypes can survive
+ * as ParticleIDs=0 until the remapping routine.
+ *
+ * @param gparts The array of loaded gparts.
+ * @param Ngparts Number of loaded gparts.
+ */
+void set_ids_to_one(struct gpart* gparts, const size_t Ngparts) {
+  for (size_t i = 0; i < Ngparts; i++) gparts[i].id_or_neg_offset = 1;
+}
diff --git a/src/common_io.h b/src/common_io.h
index 920974db86..4dc227d104 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -198,5 +198,6 @@ int get_ptype_fields(const int ptype, struct io_props* list,
                      const int with_cosmology, const int with_fof,
                      const int with_stf);
 int get_param_ptype(const char* name);
+void set_ids_to_one(struct gpart* gparts, const size_t Ngparts);
 
 #endif /* SWIFT_COMMON_IO_H */
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 6a1d004e49..0a2bc934be 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -733,7 +733,7 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
                       int with_stars, int with_black_holes, int with_cosmology,
                       int cleanup_h, int cleanup_sqrt_a, double h, double a,
                       int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
-                      int n_threads, int dry_run) {
+                      int n_threads, int dry_run, int remap_ids) {
 
   hid_t h_file = 0, h_grp = 0;
   /* GADGET has only cubic boxes (in cosmological mode) */
@@ -1010,15 +1010,23 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
 
     /* Read everything */
     if (!dry_run)
-      for (int i = 0; i < num_fields; ++i)
+      for (int i = 0; i < num_fields; ++i) {
+        /* If we are remapping ParticleIDs later, don't need to read them. */
+        if (remap_ids && strcmp(list[i].name, "ParticleIDs") == 0) continue;
+
+        /* Read array. */
         read_array_parallel(h_grp, list[i], Nparticles, N_total[ptype],
                             mpi_rank, offset[ptype], internal_units, ic_units,
                             cleanup_h, cleanup_sqrt_a, h, a);
+      }
 
     /* Close particle group */
     H5Gclose(h_grp);
   }
 
+  /* If we are remapping ParticleIDs later, start by setting them to 1. */
+  if (remap_ids) set_ids_to_one(*gparts, *Ngparts);
+
   if (!dry_run && with_gravity) {
 
     /* Let's initialise a bit of thread parallelism here */
diff --git a/src/parallel_io.h b/src/parallel_io.h
index 9b5795fe28..07c013d99b 100644
--- a/src/parallel_io.h
+++ b/src/parallel_io.h
@@ -43,7 +43,7 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
                       int with_stars, int with_black_holes, int with_cosmology,
                       int cleanup_h, int cleanup_sqrt_a, double h, double a,
                       int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
-                      int nr_threads, int dry_run);
+                      int nr_threads, int dry_run, int remap_ids);
 
 void write_output_parallel(struct engine* e,
                            const struct unit_system* internal_units,
diff --git a/src/serial_io.c b/src/serial_io.c
index 16ed7e7801..b5559b91f7 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -510,7 +510,7 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
                     int with_black_holes, int with_cosmology, int cleanup_h,
                     int cleanup_sqrt_a, double h, double a, int mpi_rank,
                     int mpi_size, MPI_Comm comm, MPI_Info info, int n_threads,
-                    int dry_run) {
+                    int dry_run, int remap_ids) {
 
   hid_t h_file = 0, h_grp = 0;
   /* GADGET has only cubic boxes (in cosmological mode) */
@@ -811,10 +811,16 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
 
         /* Read everything */
         if (!dry_run)
-          for (int i = 0; i < num_fields; ++i)
+          for (int i = 0; i < num_fields; ++i) {
+            /* If we are remapping ParticleIDs later, don't need to read them.
+             */
+            if (remap_ids && strcmp(list[i].name, "ParticleIDs") == 0) continue;
+
+            /* Read array. */
             read_array_serial(h_grp, list[i], Nparticles, N_total[ptype],
                               offset[ptype], internal_units, ic_units,
                               cleanup_h, cleanup_sqrt_a, h, a);
+          }
 
         /* Close particle group */
         H5Gclose(h_grp);
@@ -828,6 +834,9 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
     MPI_Barrier(comm);
   }
 
+  /* If we are remapping ParticleIDs later, start by setting them to 1. */
+  if (remap_ids) set_ids_to_one(*gparts, *Ngparts);
+
   /* Duplicate the parts for gravity */
   if (!dry_run && with_gravity) {
 
diff --git a/src/serial_io.h b/src/serial_io.h
index 537340b0b1..c358281c39 100644
--- a/src/serial_io.h
+++ b/src/serial_io.h
@@ -45,7 +45,7 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
                     int with_black_holes, int with_cosmology, int cleanup_h,
                     int cleanup_sqrt_a, double h, double a, int mpi_rank,
                     int mpi_size, MPI_Comm comm, MPI_Info info, int n_threads,
-                    int dry_run);
+                    int dry_run, int remap_ids);
 
 void write_output_serial(struct engine* e,
                          const struct unit_system* internal_units,
diff --git a/src/single_io.c b/src/single_io.c
index 0657671f75..073d696248 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -418,7 +418,7 @@ void read_ic_single(const char* fileName,
                     int with_gravity, int with_sink, int with_stars,
                     int with_black_holes, int with_cosmology, int cleanup_h,
                     int cleanup_sqrt_a, double h, double a, int n_threads,
-                    int dry_run) {
+                    int dry_run, int remap_ids) {
 
   hid_t h_file = 0, h_grp = 0;
   /* GADGET has only cubic boxes (in cosmological mode) */
@@ -681,14 +681,22 @@ void read_ic_single(const char* fileName,
 
     /* Read everything */
     if (!dry_run)
-      for (int i = 0; i < num_fields; ++i)
+      for (int i = 0; i < num_fields; ++i) {
+        /* If we are remapping ParticleIDs later, don't need to read them. */
+        if (remap_ids && strcmp(list[i].name, "ParticleIDs") == 0) continue;
+
+        /* Read array. */
         read_array_single(h_grp, list[i], Nparticles, internal_units, ic_units,
                           cleanup_h, cleanup_sqrt_a, h, a);
+      }
 
     /* Close particle group */
     H5Gclose(h_grp);
   }
 
+  /* If we are remapping ParticleIDs later, start by setting them to 1. */
+  if (remap_ids) set_ids_to_one(*gparts, *Ngparts);
+
   /* Duplicate the parts for gravity */
   if (!dry_run && with_gravity) {
 
diff --git a/src/single_io.h b/src/single_io.h
index 5f03fb9669..c393aa5d36 100644
--- a/src/single_io.h
+++ b/src/single_io.h
@@ -40,7 +40,7 @@ void read_ic_single(const char* fileName,
                     int with_gravity, int with_sinks, int with_stars,
                     int with_black_holes, int with_cosmology, int cleanup_h,
                     int cleanup_sqrt_a, double h, double a, int nr_threads,
-                    int dry_run);
+                    int dry_run, int remap_ids);
 
 void write_output_single(struct engine* e,
                          const struct unit_system* internal_units,
diff --git a/src/space.c b/src/space.c
index 019b8a1e29..427e9ca9c2 100644
--- a/src/space.c
+++ b/src/space.c
@@ -5954,7 +5954,8 @@ void space_remap_ids(struct space *s, int nr_nodes, int verbose) {
   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;
+  const size_t local_nr_dm =
+      local_nr_gparts > 0 ? local_nr_gparts - local_nr_baryons : 0;
 
   /* Get the global offsets */
   long long offset_parts = 0;
-- 
GitLab