diff --git a/examples/main.c b/examples/main.c
index b3577b0e2adec7650bb6acd73b128552071f47d8..2c78d6b1aab0d6fbad630266a4762957a7810bda 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -861,7 +861,7 @@ int main(int argc, char *argv[]) {
fflush(stdout);
/* Get ready to read particles of all kinds */
- size_t Ngas = 0, Ngpart = 0, Nspart = 0, Nbpart = 0;
+ size_t Ngas = 0, Ngpart = 0, Ngpart_background = 0, Nspart = 0, Nbpart = 0;
double dim[3] = {0., 0., 0.};
if (myrank == 0) clocks_gettime(&tic);
#if defined(HAVE_HDF5)
@@ -883,10 +883,11 @@ int main(int argc, char *argv[]) {
#endif
#else
read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
- &Ngas, &Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs,
- with_hydro, (with_external_gravity || with_self_gravity),
- with_stars, with_black_holes, cleanup_h, cleanup_sqrt_a,
- cosmo.h, cosmo.a, nr_threads, dry_run);
+ &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
+ &flag_entropy_ICs, with_hydro,
+ (with_external_gravity || with_self_gravity), with_stars,
+ with_black_holes, cleanup_h, cleanup_sqrt_a, cosmo.h,
+ cosmo.a, nr_threads, dry_run);
#endif
#endif
if (myrank == 0) {
@@ -918,23 +919,35 @@ int main(int argc, char *argv[]) {
#endif
/* Get the total number of particles across all nodes. */
- long long N_total[4] = {0, 0, 0, 0};
+ long long N_total[swift_type_count + 1] = {0};
+ const long long Nbaryons = Ngas + Nspart + Nbpart;
#if defined(WITH_MPI)
- long long N_long[4] = {Ngas, Ngpart, Nspart, Nbpart};
- MPI_Allreduce(&N_long, &N_total, 4, MPI_LONG_LONG_INT, MPI_SUM,
- MPI_COMM_WORLD);
+ long long N_long[swift_type_count + 1] = {0};
+ N_long[swift_type_gas] = Ngas;
+ N_long[swift_type_dark_matter] = Ngpart - Ngpart_background - Nbaryons;
+ N_long[swift_type_dark_matter_background] = Ngpart_background;
+ N_long[swift_type_stars] = Nspart;
+ N_long[swift_type_black_hole] = Nbpart;
+ N_long[swift_type_count] = Ngpart;
+ MPI_Allreduce(&N_long, &N_total, swift_type_count + 1, MPI_LONG_LONG_INT,
+ MPI_SUM, MPI_COMM_WORLD);
#else
- N_total[0] = Ngas;
- N_total[1] = Ngpart;
- N_total[2] = Nspart;
- N_total[3] = Nbpart;
+ N_total[swift_type_gas] = Ngas;
+ N_total[swift_type_dark_matter] = Ngpart - Ngpart_background - Nbaryons;
+ N_total[swift_type_dark_matter_background] = Ngpart_background;
+ N_total[swift_type_stars] = Nspart;
+ N_total[swift_type_black_hole] = Nbpart;
+ N_total[swift_type_count] = Ngpart;
#endif
if (myrank == 0)
message(
"Read %lld gas particles, %lld stars particles, %lld black hole "
- "particles and %lld gparts from the ICs.",
- N_total[0], N_total[2], N_total[3], N_total[1]);
+ "particles, %lld DM particles and %lld DM background particles from "
+ "the ICs.",
+ N_total[swift_type_gas], N_total[swift_type_stars],
+ N_total[swift_type_black_hole], N_total[swift_type_dark_matter],
+ N_total[swift_type_dark_matter_background]);
/* Verify that the fields to dump actually exist */
if (myrank == 0) io_check_output_fields(params, N_total);
diff --git a/src/common_io.c b/src/common_io.c
index 945badd14cb814e5c06c5e0edf2b983669259850..88d67d4785f248e9a71cc8b9eef3139663ab8386 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -24,13 +24,11 @@
/* This object's header. */
#include "common_io.h"
-/* First common header */
-#include "engine.h"
-
/* Local includes. */
#include "black_holes_io.h"
#include "chemistry_io.h"
#include "cooling_io.h"
+#include "engine.h"
#include "error.h"
#include "fof_io.h"
#include "gravity_io.h"
@@ -1551,6 +1549,42 @@ void io_prepare_dm_gparts(struct threadpool* tp, struct gpart* const gparts,
sizeof(struct gpart), 0, NULL);
}
+void io_prepare_dm_background_gparts_mapper(void* restrict data, int Ndm,
+ void* dummy) {
+
+ struct gpart* restrict gparts = (struct gpart*)data;
+
+ /* Let's give all these gparts a negative id */
+ for (int i = 0; i < Ndm; ++i) {
+
+ /* Negative ids are not allowed */
+ if (gparts[i].id_or_neg_offset < 0)
+ error("Negative ID for DM particle %i: ID=%lld", i,
+ gparts[i].id_or_neg_offset);
+
+ /* Set gpart type */
+ gparts[i].type = swift_type_dark_matter_background;
+ }
+}
+
+/**
+ * @brief Prepare the DM backgorund particles (in gparts) read in
+ * for the addition of the other particle types
+ *
+ * This function assumes that the DM particles are all at the start of the
+ * gparts array and that the background particles directly follow them.
+ *
+ * @param tp The current #threadpool.
+ * @param gparts The array of #gpart freshly read in.
+ * @param Ndm The number of DM particles read in.
+ */
+void io_prepare_dm_background_gparts(struct threadpool* tp,
+ struct gpart* const gparts, size_t Ndm) {
+
+ threadpool_map(tp, io_prepare_dm_background_gparts_mapper, gparts, Ndm,
+ sizeof(struct gpart), 0, NULL);
+}
+
struct duplication_data {
struct part* parts;
@@ -1895,7 +1929,7 @@ void io_collect_gparts_to_write(
* @param N_total The total number of each particle type.
*/
void io_check_output_fields(const struct swift_params* params,
- const long long N_total[3]) {
+ const long long N_total[swift_type_count]) {
/* Copy N_total to array with length == 6 */
const long long nr_total[swift_type_count] = {N_total[0], N_total[1], 0,
@@ -1908,7 +1942,7 @@ void io_check_output_fields(const struct swift_params* params,
struct io_props list[100];
/* Don't do anything if no particle of this kind */
- if (nr_total[ptype] == 0) continue;
+ if (N_total[ptype] == 0) continue;
/* Gather particle fields from the particle structures */
switch (ptype) {
@@ -1932,6 +1966,12 @@ void io_check_output_fields(const struct swift_params* params,
num_fields += velociraptor_write_gparts(NULL, list + num_fields);
break;
+ case swift_type_dark_matter_background:
+ darkmatter_write_particles(NULL, list, &num_fields);
+ num_fields += fof_write_gparts(NULL, list + num_fields);
+ num_fields += velociraptor_write_gparts(NULL, list + num_fields);
+ break;
+
case swift_type_stars:
stars_write_particles(NULL, list, &num_fields, /*with_cosmology=*/1);
num_fields += chemistry_write_sparticles(NULL, list + num_fields);
@@ -2045,6 +2085,12 @@ void io_write_output_field_parameter(const char* filename) {
num_fields += velociraptor_write_gparts(NULL, list + num_fields);
break;
+ case swift_type_dark_matter_background:
+ darkmatter_write_particles(NULL, list, &num_fields);
+ num_fields += fof_write_gparts(NULL, list + num_fields);
+ num_fields += velociraptor_write_gparts(NULL, list + num_fields);
+ break;
+
case swift_type_stars:
stars_write_particles(NULL, list, &num_fields, /*with_cosmology=*/1);
num_fields += chemistry_write_sparticles(NULL, list + num_fields);
diff --git a/src/common_io.h b/src/common_io.h
index 03437d1a38a09e59e6abd4b9ca762f181ffc3731..ca28d60901535422dd842f4554480d083c05e622 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -127,6 +127,8 @@ void io_collect_gparts_to_write(const struct gpart* restrict gparts,
const size_t Ngparts_written, int with_stf);
void io_prepare_dm_gparts(struct threadpool* tp, struct gpart* const gparts,
size_t Ndm);
+void io_prepare_dm_background_gparts(struct threadpool* tp,
+ struct gpart* const gparts, size_t Ndm);
void io_duplicate_hydro_gparts(struct threadpool* tp, struct part* const parts,
struct gpart* const gparts, size_t Ngas,
size_t Ndm);
diff --git a/src/part.c b/src/part.c
index 7d967caad18caf0c45d2a422fb4766d34f559d00..a6475983dcf5d8c3dfca239b85a56d32ef80a771 100644
--- a/src/part.c
+++ b/src/part.c
@@ -188,6 +188,15 @@ void part_verify_links(struct part *parts, struct gpart *gparts,
error("DM gpart particle linked to something !");
}
+ /* We have a background DM particle */
+ if (gparts[k].type == swift_type_dark_matter_background &&
+ gparts[k].time_bin != time_bin_not_created) {
+
+ /* Check that it's not linked */
+ if (gparts[k].id_or_neg_offset <= 0)
+ error("DM gpart particle linked to something !");
+ }
+
/* We have a gas particle */
else if (gparts[k].type == swift_type_gas) {
diff --git a/src/part_type.c b/src/part_type.c
index 75a3f270d06576831833cf75c9bed9b7fde2c68d..5a2eee2ee2e4296f4abdd1dd3ae4927748484b15 100644
--- a/src/part_type.c
+++ b/src/part_type.c
@@ -21,4 +21,4 @@
#include "part_type.h"
const char* part_type_names[swift_type_count] = {
- "Gas", "DM", "DM_boundary", "Dummy", "Stars", "BH"};
+ "Gas", "DM", "DM_background", "Dummy", "Stars", "BH"};
diff --git a/src/part_type.h b/src/part_type.h
index edab07a17080f90bee99507cd5ee3d6ed08054aa..8311f5d17a7b838b9577a029387a5f01fb1a0801 100644
--- a/src/part_type.h
+++ b/src/part_type.h
@@ -27,7 +27,7 @@
enum part_type {
swift_type_gas = 0,
swift_type_dark_matter = 1,
- swift_type_dark_matter_boundary = 2,
+ swift_type_dark_matter_background = 2,
swift_type_stars = 4,
swift_type_black_hole = 5,
swift_type_count
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index c6885746a29fd7b6bd828496316f8dad01c1b7da..d308d4712354b8266db9fb002c0cb7f3e44737db 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -1382,6 +1382,8 @@ static INLINE void runner_dopair_grav_mm_progenies(struct runner *r,
struct cell *restrict ci,
struct cell *restrict cj) {
+ return;
+
/* Loop over all pairs of progenies */
for (int i = 0; i < 8; i++) {
if (ci->progeny[i] != NULL) {
@@ -1503,6 +1505,8 @@ static INLINE void runner_dopair_recursive_grav(struct runner *r,
struct cell *ci,
struct cell *cj, int gettimer) {
+ return;
+
/* Some constants */
const struct engine *e = r->e;
const int nodeID = e->nodeID;
@@ -1657,6 +1661,8 @@ static INLINE void runner_dopair_recursive_grav(struct runner *r,
static INLINE void runner_doself_recursive_grav(struct runner *r,
struct cell *c, int gettimer) {
+ return;
+
/* Some constants */
const struct engine *e = r->e;
diff --git a/src/single_io.c b/src/single_io.c
index 76a5630120df5e946093e9f046503f73bce706b6..6be604ddf215a5be527930f45c7dc483c0cd2b9f 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -404,11 +404,11 @@ void read_ic_single(const char* fileName,
const struct unit_system* internal_units, double dim[3],
struct part** parts, struct gpart** gparts,
struct spart** sparts, struct bpart** bparts, size_t* Ngas,
- size_t* Ngparts, size_t* Nstars, size_t* Nblackholes,
- int* flag_entropy, int with_hydro, int with_gravity,
- int with_stars, int with_black_holes, int cleanup_h,
- int cleanup_sqrt_a, double h, double a, int n_threads,
- int dry_run) {
+ size_t* Ngparts, size_t* Ngparts_background, size_t* Nstars,
+ size_t* Nblackholes, int* flag_entropy, int with_hydro,
+ int with_gravity, int with_stars, int with_black_holes,
+ int cleanup_h, int cleanup_sqrt_a, double h, double a,
+ int n_threads, int dry_run) {
hid_t h_file = 0, h_grp = 0;
/* GADGET has only cubic boxes (in cosmological mode) */
@@ -419,9 +419,11 @@ void read_ic_single(const char* fileName,
size_t N[swift_type_count] = {0};
int dimension = 3; /* Assume 3D if nothing is specified */
size_t Ndm = 0;
+ size_t Ndm_background = 0;
/* Initialise counters */
- *Ngas = 0, *Ngparts = 0, *Nstars = 0, *Nblackholes = 0;
+ *Ngas = 0, *Ngparts = 0, *Ngparts_background = 0, *Nstars = 0,
+ *Nblackholes = 0;
/* Open file */
/* message("Opening file '%s' as IC.", fileName); */
@@ -563,10 +565,13 @@ void read_ic_single(const char* fileName,
/* Allocate memory to store all gravity particles */
if (with_gravity) {
Ndm = N[swift_type_dark_matter];
+ Ndm_background = N[swift_type_dark_matter_background];
*Ngparts = (with_hydro ? N[swift_type_gas] : 0) +
N[swift_type_dark_matter] +
+ N[swift_type_dark_matter_background] +
(with_stars ? N[swift_type_stars] : 0) +
(with_black_holes ? N[swift_type_black_hole] : 0);
+ *Ngparts_background = Ndm_background;
if (swift_memalign("gparts", (void**)gparts, gpart_align,
*Ngparts * sizeof(struct gpart)) != 0)
error("Error while allocating memory for gravity particles");
@@ -615,6 +620,13 @@ void read_ic_single(const char* fileName,
}
break;
+ case swift_type_dark_matter_background:
+ if (with_gravity) {
+ Nparticles = Ndm_background;
+ darkmatter_read_particles(*gparts + Ndm, list, &num_fields);
+ }
+ break;
+
case swift_type_stars:
if (with_stars) {
Nparticles = *Nstars;
@@ -653,17 +665,23 @@ void read_ic_single(const char* fileName,
/* Prepare the DM particles */
io_prepare_dm_gparts(&tp, *gparts, Ndm);
+ /* Prepare the DM background particles */
+ io_prepare_dm_background_gparts(&tp, *gparts + Ndm, Ndm_background);
+
/* Duplicate the hydro particles into gparts */
- if (with_hydro) io_duplicate_hydro_gparts(&tp, *parts, *gparts, *Ngas, Ndm);
+ if (with_hydro)
+ io_duplicate_hydro_gparts(&tp, *parts, *gparts, *Ngas,
+ Ndm + Ndm_background);
/* Duplicate the star particles into gparts */
if (with_stars)
- io_duplicate_stars_gparts(&tp, *sparts, *gparts, *Nstars, Ndm + *Ngas);
+ io_duplicate_stars_gparts(&tp, *sparts, *gparts, *Nstars,
+ Ndm + Ndm_background + *Ngas);
/* Duplicate the black hole particles into gparts */
if (with_black_holes)
io_duplicate_black_holes_gparts(&tp, *bparts, *gparts, *Nblackholes,
- Ndm + *Ngas + *Nstars);
+ Ndm + Ndm_background + *Ngas + *Nstars);
threadpool_clean(&tp);
}
diff --git a/src/single_io.h b/src/single_io.h
index 9ff04c893378d7f8c01621e7dbf387dc939d0157..47b575272c2fbe5b170eb4370743cd6830b99ae0 100644
--- a/src/single_io.h
+++ b/src/single_io.h
@@ -34,11 +34,11 @@ void read_ic_single(const char* fileName,
const struct unit_system* internal_units, double dim[3],
struct part** parts, struct gpart** gparts,
struct spart** sparts, struct bpart** bparts, size_t* Ngas,
- size_t* Ndm, size_t* Nstars, size_t* Nblackholes,
- int* flag_entropy, int with_hydro, int with_gravity,
- int with_stars, int with_black_holes, int cleanup_h,
- int cleanup_sqrt_a, double h, double a, int nr_threads,
- int dry_run);
+ size_t* Ndm, size_t* Ndm_background, size_t* Nstars,
+ size_t* Nblackholes, int* flag_entropy, int with_hydro,
+ int with_gravity, int with_stars, int with_black_holes,
+ int cleanup_h, int cleanup_sqrt_a, double h, double a,
+ int nr_threads, int dry_run);
void write_output_single(struct engine* e, const char* baseName,
const struct unit_system* internal_units,
diff --git a/src/space.c b/src/space.c
index 49c714743728a3b93f631434d87de69a61dd716f..b25750b0d648a5123454a8d2551fb8999c6fa4cb 100644
--- a/src/space.c
+++ b/src/space.c
@@ -3945,6 +3945,9 @@ void space_synchronize_particle_positions_mapper(void *map_data, int nr_gparts,
if (gp->type == swift_type_dark_matter)
continue;
+ else if (gp->type == swift_type_dark_matter_background)
+ continue;
+
else if (gp->type == swift_type_gas) {
/* Get its gassy friend */
@@ -3988,6 +3991,10 @@ void space_synchronize_particle_positions_mapper(void *map_data, int nr_gparts,
gp->mass = bp->mass;
}
+
+ else {
+ error("Invalid type!");
+ }
}
}
diff --git a/src/velociraptor_interface.c b/src/velociraptor_interface.c
index 47c107d047e3abfcac32571f780cf52649a9c38d..2599c7b6ec4df05c0ae701e33474ffb126f3d6e3 100644
--- a/src/velociraptor_interface.c
+++ b/src/velociraptor_interface.c
@@ -283,6 +283,13 @@ void velociraptor_convert_particles_mapper(void *map_data, int nr_gparts,
swift_parts[i].T = 0.f;
break;
+ case swift_type_dark_matter_background:
+
+ swift_parts[i].id = gparts[i].id_or_neg_offset;
+ swift_parts[i].u = 0.f;
+ swift_parts[i].T = 0.f;
+ break;
+
default:
error("Particle type not handled by VELOCIraptor.");
}