diff --git a/examples/MultiTypes/multiTypes.yml b/examples/MultiTypes/multiTypes.yml
index 4d54f95fcdd09464b03d0f9987398cd2710b2e44..05d45595e9ec43b0849ed574f6b9922c19021a33 100644
--- a/examples/MultiTypes/multiTypes.yml
+++ b/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:
diff --git a/examples/MultiTypes/run.sh b/examples/MultiTypes/run.sh
index 508a5097f8961f446a51204e889875e33d4f634e..24c87b77a0bab83a8b5f153f2787382998f1e480 100755
--- a/examples/MultiTypes/run.sh
+++ b/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
diff --git a/examples/main.c b/examples/main.c
index cf0c37eda2ec5a4a14ebe3f5191030124d2fe6ed..9abbb5433666f3674e8c59389e72dc1d25dfe3c6 100644
--- a/examples/main.c
+++ b/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];
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 9681e880e9f4b44d4f77195eecf5b8642266ce52..31b948ea4c59de42bce7c3b80b8ae2526b10e0cf 100644
--- a/examples/parameter_example.yml
+++ b/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:
diff --git a/src/space.c b/src/space.c
index 802dc30d1bcd44d4cce46b2a803afade07f5d685..503727cfef81bffabd1d1cf568be6954bad175cf 100644
--- a/src/space.c
+++ b/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
*
diff --git a/src/space.h b/src/space.h
index a25149e8fe6971b24856a2a60cae23747fbc56ac..327f0225085663d1859bcdbd82a6ac2583fc2b43 100644
--- a/src/space.h
+++ b/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 */