Added the option to replicate the particles along the axis to create...

Added the option to replicate the particles along the axis to create artificially larger test cases. Driven by the parameter 'InitialConditions:replicate'.

examples/MultiTypes/multiTypes.yml

@@ -32,6 +32,7 @@ SPH:
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  ./multiTypes.hdf5     # The file to read
+  replicate:  2                     # Replicate all particles twice along each axis
 
 # External potential parameters
 PointMassPotential:

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 17 24 12
+    python makeIC.py 9 13 7
 fi
 
 ../swift -s -g -S -t 1 multiTypes.yml 2>&1 | tee output.log

examples/main.c

@@ -371,6 +371,8 @@ int main(int argc, char *argv[]) {
   /* Read particles and space information from (GADGET) ICs */
   char ICfileName[200] = "";
   parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
+  const int replicate =
+      parser_get_opt_param_int(params, "InitialConditions:replicate", 1);
   if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
   fflush(stdout);
 
@@ -441,7 +443,7 @@ int main(int argc, char *argv[]) {
   if (myrank == 0) clocks_gettime(&tic);
   struct space s;
   space_init(&s, params, dim, parts, gparts, sparts, Ngas, Ngpart, Nspart,
-             periodic, with_self_gravity, talking, dry_run);
+             periodic, replicate, with_self_gravity, talking, dry_run);
   if (myrank == 0) {
     clocks_gettime(&toc);
     message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
@@ -458,6 +460,7 @@ int main(int argc, char *argv[]) {
             s.cdim[1], s.cdim[2]);
     message("%zi parts in %i cells.", s.nr_parts, s.tot_cells);
     message("%zi gparts in %i cells.", s.nr_gparts, s.tot_cells);
+    message("%zi sparts in %i cells.", s.nr_sparts, s.tot_cells);
     message("maximum depth is %d.", s.maxdepth);
     fflush(stdout);
   }
@@ -521,6 +524,18 @@ int main(int argc, char *argv[]) {
   for (k = 0; k < runner_hist_N; k++) runner_hist_bins[k] = 0;
 #endif
 
+#if defined(WITH_MPI)
+  N_long[0] = s.nr_parts;
+  N_long[1] = s.nr_gparts;
+  N_long[2] = s.nr_sparts;
+  MPI_Reduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM, 0,
+             MPI_COMM_WORLD);
+#else
+  N_total[0] = s.nr_parts;
+  N_total[1] = s.nr_gparts;
+  N_total[2] = s.nr_sparts;
+#endif
+
   /* Get some info to the user. */
   if (myrank == 0) {
     long long N_DM = N_total[1] - N_total[2] - N_total[0];

examples/parameter_example.yml

@@ -55,6 +55,7 @@ InitialConditions:
   shift_x:    0.                    # (Optional) A shift to apply to all particles read from the ICs (in internal units).
   shift_y:    0.
   shift_z:    0.
+  replicate:  2                     # (Optional) Replicate all particles along each axis a given number of times. Default 1.
 
 # Parameters governing domain decomposition
 DomainDecomposition:

src/space.c

@@ -2358,6 +2358,7 @@ void space_init_sparts(struct space *s) {
  * @param Ngpart The number of Gravity particles in the space.
  * @param Nspart The number of star 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 gravity flag whether we are doing gravity or not.
  * @param verbose Print messages to stdout or not.
  * @param dry_run If 1, just initialise stuff, don't do anything with the parts.
@@ -2370,8 +2371,8 @@ void space_init_sparts(struct space *s) {
 void space_init(struct space *s, const struct swift_params *params,
                 double dim[3], struct part *parts, struct gpart *gparts,
                 struct spart *sparts, size_t Npart, size_t Ngpart,
-                size_t Nspart, int periodic, int gravity, int verbose,
-                int dry_run) {
+                size_t Nspart, int periodic, int replicate, int gravity,
+                int verbose, int dry_run) {
 
   /* Clean-up everything */
   bzero(s, sizeof(struct space));
@@ -2394,15 +2395,29 @@ void space_init(struct space *s, const struct swift_params *params,
   s->sparts = sparts;
   s->nr_queues = 1; /* Temporary value until engine construction */
 
+  /* Are we replicating the space ? */
+  if (replicate < 1)
+    error("Value of 'InitialConditions:replicate' (%d) is too small",
+          replicate);
+  if (replicate > 1) {
+    space_replicate(s, replicate, verbose);
+    parts = s->parts;
+    gparts = s->gparts;
+    sparts = s->sparts;
+    Npart = s->nr_parts;
+    Ngpart = s->nr_gparts;
+    Nspart = s->nr_sparts;
+  }
+
   /* Decide on the minimal top-level cell size */
-  const double dmax = max(max(dim[0], dim[1]), dim[2]);
+  const double dmax = max(max(s->dim[0], s->dim[1]), s->dim[2]);
   int maxtcells =
       parser_get_opt_param_int(params, "Scheduler:max_top_level_cells",
                                space_max_top_level_cells_default);
   s->cell_min = 0.99 * dmax / maxtcells;
 
   /* Check that it is big enough. */
-  const double dmin = min(min(dim[0], dim[1]), dim[2]);
+  const double dmin = min(min(s->dim[0], s->dim[1]), s->dim[2]);
   int needtcells = 3 * dmax / dmin;
   if (maxtcells < needtcells)
     error(
@@ -2464,13 +2479,13 @@ void space_init(struct space *s, const struct swift_params *params,
     if (periodic) {
       for (size_t k = 0; k < Npart; k++)
         for (int j = 0; j < 3; j++) {
-          while (parts[k].x[j] < 0) parts[k].x[j] += dim[j];
-          while (parts[k].x[j] >= dim[j]) parts[k].x[j] -= dim[j];
+          while (parts[k].x[j] < 0) parts[k].x[j] += s->dim[j];
+          while (parts[k].x[j] >= s->dim[j]) parts[k].x[j] -= s->dim[j];
         }
     } else {
       for (size_t k = 0; k < Npart; k++)
         for (int j = 0; j < 3; j++)
-          if (parts[k].x[j] < 0 || parts[k].x[j] >= dim[j])
+          if (parts[k].x[j] < 0 || parts[k].x[j] >= s->dim[j])
             error("Not all particles are within the specified domain.");
     }
 
@@ -2478,13 +2493,13 @@ void space_init(struct space *s, const struct swift_params *params,
     if (periodic) {
       for (size_t k = 0; k < Ngpart; k++)
         for (int j = 0; j < 3; j++) {
-          while (gparts[k].x[j] < 0) gparts[k].x[j] += dim[j];
-          while (gparts[k].x[j] >= dim[j]) gparts[k].x[j] -= dim[j];
+          while (gparts[k].x[j] < 0) gparts[k].x[j] += s->dim[j];
+          while (gparts[k].x[j] >= s->dim[j]) gparts[k].x[j] -= s->dim[j];
         }
     } else {
       for (size_t k = 0; k < Ngpart; k++)
         for (int j = 0; j < 3; j++)
-          if (gparts[k].x[j] < 0 || gparts[k].x[j] >= dim[j])
+          if (gparts[k].x[j] < 0 || gparts[k].x[j] >= s->dim[j])
             error("Not all g-particles are within the specified domain.");
     }
 
@@ -2492,13 +2507,13 @@ void space_init(struct space *s, const struct swift_params *params,
     if (periodic) {
       for (size_t k = 0; k < Nspart; k++)
         for (int j = 0; j < 3; j++) {
-          while (sparts[k].x[j] < 0) sparts[k].x[j] += dim[j];
-          while (sparts[k].x[j] >= dim[j]) sparts[k].x[j] -= dim[j];
+          while (sparts[k].x[j] < 0) sparts[k].x[j] += s->dim[j];
+          while (sparts[k].x[j] >= s->dim[j]) sparts[k].x[j] -= s->dim[j];
         }
     } else {
       for (size_t k = 0; k < Nspart; k++)
         for (int j = 0; j < 3; j++)
-          if (sparts[k].x[j] < 0 || sparts[k].x[j] >= dim[j])
+          if (sparts[k].x[j] < 0 || sparts[k].x[j] >= s->dim[j])
             error("Not all s-particles are within the specified domain.");
     }
   }
@@ -2524,6 +2539,127 @@ void space_init(struct space *s, const struct swift_params *params,
   if (!dry_run) space_regrid(s, verbose);
 }
 
+/**
+ * @brief Replicate the content of a space along each axis.
+ *
+ * Should only be called during initialisation.
+ *
+ * @param s The #space to replicate.
+ * @param replicate The number of copies along each axis.
+ * @param verbose Are we talkative ?
+ */
+void space_replicate(struct space *s, int replicate, int verbose) {
+
+  if (replicate < 1) error("Invalid replicate value: %d", replicate);
+
+  message("Replicating space %d times along each axis.", replicate);
+
+  const int factor = replicate * replicate * replicate;
+
+  /* Store the current values */
+  const size_t nr_parts = s->nr_parts;
+  const size_t nr_gparts = s->nr_gparts;
+  const size_t nr_sparts = s->nr_sparts;
+  const size_t nr_dm = nr_gparts - nr_parts - nr_sparts;
+
+  s->size_parts = s->nr_parts = nr_parts * factor;
+  s->size_gparts = s->nr_gparts = nr_gparts * factor;
+  s->size_sparts = s->nr_sparts = nr_sparts * factor;
+
+  /* Allocate space for new particles */
+  struct part *parts = NULL;
+  struct gpart *gparts = NULL;
+  struct spart *sparts = NULL;
+
+  if (posix_memalign((void *)&parts, part_align,
+                     s->nr_parts * sizeof(struct part)) != 0)
+    error("Failed to allocate new part array.");
+
+  if (posix_memalign((void *)&gparts, gpart_align,
+                     s->nr_gparts * sizeof(struct gpart)) != 0)
+    error("Failed to allocate new gpart array.");
+
+  if (posix_memalign((void *)&sparts, spart_align,
+                     s->nr_sparts * sizeof(struct spart)) != 0)
+    error("Failed to allocate new spart array.");
+
+  /* Replicate everything */
+  for (int i = 0; i < replicate; ++i) {
+    for (int j = 0; j < replicate; ++j) {
+      for (int k = 0; k < replicate; ++k) {
+        const size_t offset = i * replicate * replicate + j * replicate + k;
+
+        /* First copy the data */
+        memcpy(parts + offset * nr_parts, s->parts,
+               nr_parts * sizeof(struct part));
+        memcpy(sparts + offset * nr_sparts, s->sparts,
+               nr_sparts * sizeof(struct spart));
+        memcpy(gparts + offset * nr_gparts, s->gparts,
+               nr_gparts * sizeof(struct gpart));
+
+        /* Shift the positions */
+        const double shift[3] = {i * s->dim[0], j * s->dim[1], k * s->dim[2]};
+
+        for (size_t n = offset * nr_parts; n < (offset + 1) * nr_parts; ++n) {
+          parts[n].x[0] += shift[0];
+          parts[n].x[1] += shift[1];
+          parts[n].x[2] += shift[2];
+        }
+        for (size_t n = offset * nr_gparts; n < (offset + 1) * nr_gparts; ++n) {
+          gparts[n].x[0] += shift[0];
+          gparts[n].x[1] += shift[1];
+          gparts[n].x[2] += shift[2];
+        }
+        for (size_t n = offset * nr_sparts; n < (offset + 1) * nr_sparts; ++n) {
+          sparts[n].x[0] += shift[0];
+          sparts[n].x[1] += shift[1];
+          sparts[n].x[2] += shift[2];
+        }
+
+        /* Set the correct links (recall gpart are sorted by type at start-up):
+           first DM (unassociated gpart), then gas, then stars */
+        if (nr_parts > 0 && nr_gparts > 0) {
+          const size_t offset_part = offset * nr_parts;
+          const size_t offset_gpart = offset * nr_gparts + nr_dm;
+
+          for (size_t n = 0; n < nr_parts; ++n) {
+            parts[offset_part + n].gpart = &gparts[offset_gpart + n];
+            gparts[offset_gpart + n].id_or_neg_offset = -(offset_part + n);
+          }
+        }
+        if (nr_sparts > 0 && nr_gparts > 0) {
+          const size_t offset_spart = offset * nr_sparts;
+          const size_t offset_gpart = offset * nr_gparts + nr_dm + nr_parts;
+
+          for (size_t n = 0; n < nr_sparts; ++n) {
+            sparts[offset_spart + n].gpart = &gparts[offset_gpart + n];
+            gparts[offset_gpart + n].id_or_neg_offset = -(offset_spart + n);
+          }
+        }
+      }
+    }
+  }
+
+  /* Replace the content of the space */
+  free(s->parts);
+  free(s->gparts);
+  free(s->sparts);
+  s->parts = parts;
+  s->gparts = gparts;
+  s->sparts = sparts;
+
+  /* Finally, update the domain size */
+  s->dim[0] *= replicate;
+  s->dim[1] *= replicate;
+  s->dim[2] *= replicate;
+
+#ifdef SWIFT_DEBUG_CHECKS
+  /* Verify that everything is correct */
+  part_verify_links(s->parts, s->gparts, s->sparts, s->nr_parts, s->nr_gparts,
+                    s->nr_sparts, verbose);
+#endif
+}
+
 /**
  * @brief Cleans-up all the cell links in the space
  *

src/space.h

@@ -165,8 +165,8 @@ int space_getsid(struct space *s, struct cell **ci, struct cell **cj,
 void space_init(struct space *s, const struct swift_params *params,
                 double dim[3], struct part *parts, struct gpart *gparts,
                 struct spart *sparts, size_t Npart, size_t Ngpart,
-                size_t Nspart, int periodic, int gravity, int verbose,
-                int dry_run);
+                size_t Nspart, int periodic, int replicate, int gravity,
+                int verbose, int dry_run);
 void space_sanitize(struct space *s);
 void space_map_cells_pre(struct space *s, int full,
                          void (*fun)(struct cell *c, void *data), void *data);
@@ -206,6 +206,7 @@ void space_init_sparts(struct space *s);
 void space_link_cleanup(struct space *s);
 void space_check_drift_point(struct space *s, integertime_t ti_current);
 void space_check_timesteps(struct space *s);
+void space_replicate(struct space *s, int replicate, int verbose);
 void space_clean(struct space *s);
 
 #endif /* SWIFT_SPACE_H */