Merge branch 'nu_replicate' into 'master'

Allow space replication and id remapping for neutrino particles

Closes #837

See merge request !1673

doc/RTD/source/ParameterFiles/parameter_description.rst

@@ -706,6 +706,14 @@ creation of a compact range of IDs beyond which the new IDs generated by
 splitting can be safely drawn from. Note that, when ``remap_ids`` is
 switched on the ICs do not need to contain a ``ParticleIDs`` field.
 
+Both replication and remapping explicitly overwrite any particle IDs
+provided in the initial conditions. This may cause problems for runs
+with neutrino particles, as some models assume that that the particle
+ID was used as a random seed for the Fermi-Dirac momentum. In this case,
+the ``Neutrino:generate_ics`` option can be used to generate new initial
+conditions based on the replicated or remapped IDs. See :ref:`Neutrinos`
+for details.
+
 * Name of a HDF5 group to copy from the ICs file(s): ``metadata_group_name`` (default: ``ICs_parameters``)
 
 If the initial conditions generator writes a HDF5 group with the parameters

src/space.c

@@ -1129,8 +1129,6 @@ void space_init(struct space *s, struct swift_params *params,
   if (replicate > 1) {
     if (with_DM_background)
       error("Can't replicate the space if background DM particles are in use.");
-    if (neutrinos)
-      error("Can't replicate the space if neutrino DM particles are in use.");
 
     space_replicate(s, replicate, verbose);
     parts = s->parts;
@@ -1439,6 +1437,7 @@ void space_replicate(struct space *s, int replicate, int verbose) {
   const size_t nr_sparts = s->nr_sparts;
   const size_t nr_bparts = s->nr_bparts;
   const size_t nr_sinks = s->nr_sinks;
+  const size_t nr_nuparts = s->nr_nuparts;
   const size_t nr_dm = nr_gparts - nr_parts - nr_sparts - nr_bparts;
 
   s->size_parts = s->nr_parts = nr_parts * factor;
@@ -1446,6 +1445,7 @@ void space_replicate(struct space *s, int replicate, int verbose) {
   s->size_sparts = s->nr_sparts = nr_sparts * factor;
   s->size_bparts = s->nr_bparts = nr_bparts * factor;
   s->size_sinks = s->nr_sinks = nr_sinks * factor;
+  s->nr_nuparts = nr_nuparts * factor;
 
   /* Allocate space for new particles */
   struct part *parts = NULL;
@@ -1605,10 +1605,10 @@ void space_remap_ids(struct space *s, int nr_nodes, int verbose) {
   if (verbose) message("Remapping all the IDs");
 
   size_t local_nr_dm_background = 0;
-  size_t local_nr_neutrino = 0;
+  size_t local_nr_nuparts = 0;
   for (size_t i = 0; i < s->nr_gparts; ++i) {
     if (s->gparts[i].type == swift_type_neutrino)
-      local_nr_neutrino++;
+      local_nr_nuparts++;
     else if (s->gparts[i].type == swift_type_dark_matter_background)
       local_nr_dm_background++;
   }
@@ -1621,10 +1621,10 @@ 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 > 0 ? local_nr_gparts - local_nr_baryons -
-                                local_nr_neutrino - local_nr_dm_background
-                          : 0;
+  const size_t local_nr_dm = local_nr_gparts > 0
+                                 ? local_nr_gparts - local_nr_baryons -
+                                       local_nr_nuparts - local_nr_dm_background
+                                 : 0;
 
   /* Get the global offsets */
   long long offset_parts = 0;
@@ -1633,6 +1633,7 @@ void space_remap_ids(struct space *s, int nr_nodes, int verbose) {
   long long offset_bparts = 0;
   long long offset_dm = 0;
   long long offset_dm_background = 0;
+  long long offset_nuparts = 0;
 #ifdef WITH_MPI
   MPI_Exscan(&local_nr_parts, &offset_parts, 1, MPI_LONG_LONG_INT, MPI_SUM,
              MPI_COMM_WORLD);
@@ -1646,6 +1647,8 @@ void space_remap_ids(struct space *s, int nr_nodes, int verbose) {
              MPI_COMM_WORLD);
   MPI_Exscan(&local_nr_dm_background, &offset_dm_background, 1,
              MPI_LONG_LONG_INT, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Exscan(&local_nr_nuparts, &offset_nuparts, 1, MPI_LONG_LONG_INT, MPI_SUM,
+             MPI_COMM_WORLD);
 #endif
 
   /* Total number of particles of each kind */
@@ -1654,6 +1657,7 @@ void space_remap_ids(struct space *s, int nr_nodes, int verbose) {
   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;
+  long long total_nuparts = offset_nuparts + local_nr_nuparts;
   // long long total_dm_backgroud = offset_dm_background +
   // local_nr_dm_background;
 
@@ -1664,22 +1668,26 @@ void space_remap_ids(struct space *s, int nr_nodes, int verbose) {
   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);
+  MPI_Bcast(&total_nuparts, 1, MPI_LONG_LONG_INT, nr_nodes - 1, MPI_COMM_WORLD);
   // MPI_Bcast(&total_dm_background, 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 then DM background
-   * Note that we leave a large gap (10x the number of particles) in-between the
-   * regular particles and the background ones. This allow for particle
-   * splitting to keep a compact set of ids. */
+   * IDs will be DM then gas then sinks than stars then BHs then nus then
+   * DM background. Note that we leave a large gap (10x the number of particles)
+   * in-between the regular particles and the background ones. This allow for
+   * particle splitting to keep a compact set of ids. */
   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;
-  offset_dm_background += 1 + 10 * (total_dm * total_parts + total_sinks +
-                                    total_sparts + total_bparts);
+  offset_nuparts +=
+      1 + total_dm + total_parts + total_sinks + total_sparts + total_bparts;
+  offset_dm_background +=
+      1 + 10 * (total_dm * total_parts + total_sinks + total_sparts +
+                total_bparts + total_nuparts);
 
   /* We can now remap the IDs in the range [offset offset + local_nr] */
   for (size_t i = 0; i < local_nr_parts; ++i) {
@@ -1696,10 +1704,14 @@ void space_remap_ids(struct space *s, int nr_nodes, int verbose) {
   }
   size_t count_dm = 0;
   size_t count_dm_background = 0;
+  size_t count_nu = 0;
   for (size_t i = 0; i < s->nr_gparts; ++i) {
     if (s->gparts[i].type == swift_type_dark_matter) {
       s->gparts[i].id_or_neg_offset = offset_dm + count_dm;
       count_dm++;
+    } else if (s->gparts[i].type == swift_type_neutrino) {
+      s->gparts[i].id_or_neg_offset = offset_nuparts + count_nu;
+      count_nu++;
     } else if (s->gparts[i].type == swift_type_dark_matter_background) {
       s->gparts[i].id_or_neg_offset =
           offset_dm_background + count_dm_background;