diff --git a/src/parallel_io.c b/src/parallel_io.c
index 929e0bb43ac15526bffe66a9563e38e5ea7af9a4..57d6235fea3040dccd0ff8f69227704e3b7070fa 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -704,8 +704,8 @@ void writeArray(struct engine* e, hid_t grp, char* fileName,
void read_ic_parallel(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,
+ size_t* Ngas, size_t* Ngparts, size_t *Ngparts_background, size_t* Nstars,
+ size_t* Nblackholes, int* flag_entropy, float* gpart_mass, 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 mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
@@ -721,9 +721,10 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
long long offset[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); */
@@ -875,11 +876,14 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
/* Allocate memory to store gravity particles */
if (with_gravity) {
- Ndm = N[1];
+ 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");
@@ -928,6 +932,13 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
}
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;
@@ -968,17 +979,27 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
/* Prepare the DM particles */
io_prepare_dm_gparts(&tp, *gparts, Ndm);
+ /* Record the mass of the DM particles */
+ const float local_gpart_mass = (*gparts)[0].mass;
+ float global_gpart_mass;
+ MPI_Allreduce(&local_gpart_mass, &global_gpart_mass, 1, MPI_FLOAT, MPI_MAX,
+ comm);
+ *gpart_mass = global_gpart_mass;
+
+ /* 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 stars 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 stars 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);
}
@@ -1236,6 +1257,18 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
}
break;
+ case swift_type_dark_matter_background:
+ darkmatter_write_particles(gparts, list, &num_fields);
+ if (with_fof) {
+ num_fields += fof_write_gparts(gparts, list + num_fields);
+ }
+ if (with_stf) {
+ num_fields += velociraptor_write_gparts(e->s->gpart_group_data,
+ list + num_fields);
+ }
+ break;
+
+
case swift_type_stars:
stars_write_particles(sparts, list, &num_fields, with_cosmology);
num_fields += chemistry_write_sparticles(sparts, list + num_fields);
@@ -1329,6 +1362,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
const int with_cooling = e->policy & engine_policy_cooling;
const int with_temperature = e->policy & engine_policy_temperature;
const int with_fof = e->policy & engine_policy_fof;
+ const int with_DM_background = e->s->with_DM_background;
#ifdef HAVE_VELOCIRAPTOR
const int with_stf = (e->policy & engine_policy_structure_finding) &&
(e->s->gpart_group_data != NULL);
@@ -1344,6 +1378,11 @@ void write_output_parallel(struct engine* e, const char* baseName,
// const size_t Nbaryons = Ngas + Nstars;
// const size_t Ndm = Ntot > 0 ? Ntot - Nbaryons : 0;
+ size_t Ndm_background = 0;
+ if (with_DM_background) {
+ Ndm_background = io_count_dm_background_gparts(gparts, Ntot);
+ }
+
/* Number of particles that we will write */
const size_t Ntot_written =
e->s->nr_gparts - e->s->nr_inhibited_gparts - e->s->nr_extra_gparts;
@@ -1356,10 +1395,10 @@ void write_output_parallel(struct engine* e, const char* baseName,
const size_t Nbaryons_written =
Ngas_written + Nstars_written + Nblackholes_written;
const size_t Ndm_written =
- Ntot_written > 0 ? Ntot_written - Nbaryons_written : 0;
+ Ntot_written > 0 ? Ntot_written - Nbaryons_written - Ndm_background : 0;
/* Compute offset in the file and total number of particles */
- size_t N[swift_type_count] = {Ngas_written, Ndm_written, 0, 0,
+ size_t N[swift_type_count] = {Ngas_written, Ndm_written, Ndm_background, 0,
Nstars_written, Nblackholes_written};
long long N_total[swift_type_count] = {0};
long long offset[swift_type_count] = {0};
diff --git a/src/parallel_io.h b/src/parallel_io.h
index 6e6b52b74ebf166068f8dff110ea0ca0fd8464a2..11c1e895aefe5c7972dc904b4601ade3b9241bd5 100644
--- a/src/parallel_io.h
+++ b/src/parallel_io.h
@@ -36,8 +36,9 @@
void read_ic_parallel(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* Nsparts,
- size_t* Nbparts, int* flag_entropy, int with_hydro,
+ size_t* Ngas, size_t* Ngparts, size_t* Ngparts_background, size_t* Nsparts,
+ size_t* Nbparts, int* flag_entropy, float* gpart_mass,
+ 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 mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
diff --git a/src/serial_io.c b/src/serial_io.c
index 05a1839572cfcea8a81559cbdb197dbed26e605d..8917690ca6f4c3ac1ed2e066bdca02fdcb52e2fb 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -677,7 +677,7 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
/* Allocate memory to store all gravity particles */
if (with_gravity) {
- Ndm = N[1];
+ Ndm = N[swift_dark_matter];
Ndm_background = N[swift_type_dark_matter_background];
*Ngparts = (with_hydro ? N[swift_type_gas] : 0) +
N[swift_type_dark_matter] +
@@ -806,7 +806,6 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
float global_gpart_mass;
MPI_Allreduce(&local_gpart_mass, &global_gpart_mass, 1, MPI_FLOAT, MPI_MAX,
comm);
- /* Record the mass of the DM particles */
*gpart_mass = global_gpart_mass;
/* Prepare the DM background particles */