diff --git a/examples/main.c b/examples/main.c index 2c78d6b1aab0d6fbad630266a4762957a7810bda..5127ab9ef169be30153ea2c033c17a23c8b5461f 100644 --- a/examples/main.c +++ b/examples/main.c @@ -949,6 +949,10 @@ int main(int argc, char *argv[]) { N_total[swift_type_black_hole], N_total[swift_type_dark_matter], N_total[swift_type_dark_matter_background]); + /* Do we have background DM particles? */ + const int with_DM_background_particles = + N_total[swift_type_dark_matter_background] > 0; + /* Verify that the fields to dump actually exist */ if (myrank == 0) io_check_output_fields(params, N_total); @@ -956,8 +960,8 @@ int main(int argc, char *argv[]) { if (myrank == 0) clocks_gettime(&tic); space_init(&s, params, &cosmo, dim, parts, gparts, sparts, bparts, Ngas, Ngpart, Nspart, Nbpart, periodic, replicate, generate_gas_in_ics, - with_hydro, with_self_gravity, with_star_formation, talking, - dry_run); + with_hydro, with_self_gravity, with_star_formation, + with_DM_background_particles, talking, dry_run); if (myrank == 0) { clocks_gettime(&toc); diff --git a/src/cell.c b/src/cell.c index f4903f024b10c7cc7867291b597a3bdf42997c27..80c5798c2b5a13bcb3c283eb16fa58d537d63bf6 100644 --- a/src/cell.c +++ b/src/cell.c @@ -5363,6 +5363,9 @@ void cell_remove_gpart(const struct engine *e, struct cell *c, if (c->nodeID != e->nodeID) error("Can't remove a particle in a foreign cell."); + if (gp->type == swift_type_dark_matter_background) + error("Can't remove a DM background particle!"); + /* Mark the particle as inhibited */ gp->time_bin = time_bin_inhibited; diff --git a/src/common_io.c b/src/common_io.c index 88d67d4785f248e9a71cc8b9eef3139663ab8386..ee548c99435aca359a131e3dbb4ec55205a4b11e 100644 --- a/src/common_io.c +++ b/src/common_io.c @@ -1585,6 +1585,20 @@ void io_prepare_dm_background_gparts(struct threadpool* tp, sizeof(struct gpart), 0, NULL); } +size_t io_count_dm_background_gparts(const struct gpart* const gparts, + const size_t Ndm) { + + swift_declare_aligned_ptr(const struct gpart, gparts_array, gparts, + gpart_align); + + size_t count = 0; + for (size_t i = 0; i < Ndm; ++i) { + if (gparts[i].type == swift_type_dark_matter_background) ++count; + } + + return count; +} + struct duplication_data { struct part* parts; @@ -1880,7 +1894,8 @@ void io_collect_bparts_to_write(const struct bpart* restrict bparts, } /** - * @brief Copy every non-inhibited DM #gpart into the gparts_written array. + * @brief Copy every non-inhibited regulat DM #gpart into the gparts_written + * array. * * @param gparts The array of #gpart containing all particles. * @param vr_data The array of gpart-related VELOCIraptor output. @@ -1922,6 +1937,50 @@ void io_collect_gparts_to_write( count, Ngparts_written); } +/** + * @brief Copy every non-inhibited background DM #gpart into the gparts_written + * array. + * + * @param gparts The array of #gpart containing all particles. + * @param vr_data The array of gpart-related VELOCIraptor output. + * @param gparts_written The array of #gpart to fill with particles we want to + * write. + * @param vr_data_written The array of gpart-related VELOCIraptor with particles + * we want to write. + * @param Ngparts The total number of #part. + * @param Ngparts_written The total number of #part to write. + * @param with_stf Are we running with STF? i.e. do we want to collect vr data? + */ +void io_collect_gparts_background_to_write( + const struct gpart* restrict gparts, + const struct velociraptor_gpart_data* restrict vr_data, + struct gpart* restrict gparts_written, + struct velociraptor_gpart_data* restrict vr_data_written, + const size_t Ngparts, const size_t Ngparts_written, const int with_stf) { + + size_t count = 0; + + /* Loop over all parts */ + for (size_t i = 0; i < Ngparts; ++i) { + + /* And collect the ones that have not been removed */ + if ((gparts[i].time_bin != time_bin_inhibited) && + (gparts[i].time_bin != time_bin_not_created) && + (gparts[i].type == swift_type_dark_matter_background)) { + + if (with_stf) vr_data_written[count] = vr_data[i]; + + gparts_written[count] = gparts[i]; + count++; + } + } + + /* Check that everything is fine */ + if (count != Ngparts_written) + error("Collected the wrong number of g-particles (%zu vs. %zu expected)", + count, Ngparts_written); +} + /** * @brief Verify the io parameter file * diff --git a/src/common_io.h b/src/common_io.h index ca28d60901535422dd842f4554480d083c05e622..aa4b0c552f8a416a4aa42ecce5e9ae0b74dbc731 100644 --- a/src/common_io.h +++ b/src/common_io.h @@ -125,10 +125,18 @@ void io_collect_gparts_to_write(const struct gpart* restrict gparts, struct velociraptor_gpart_data* vr_data_written, const size_t Ngparts, const size_t Ngparts_written, int with_stf); +void io_collect_gparts_background_to_write( + const struct gpart* restrict gparts, + const struct velociraptor_gpart_data* vr_data, + struct gpart* restrict gparts_written, + struct velociraptor_gpart_data* vr_data_written, const size_t Ngparts, + 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); +size_t io_count_dm_background_gparts(const 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/single_io.c b/src/single_io.c index 6be604ddf215a5be527930f45c7dc483c0cd2b9f..6f1f4a0d8bbcb1c66c35e6efe6f2a5a13a9f5104 100644 --- a/src/single_io.c +++ b/src/single_io.c @@ -727,6 +727,7 @@ void write_output_single(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); @@ -742,7 +743,13 @@ void write_output_single(struct engine* e, const char* baseName, // const size_t Nbaryons = Ngas + Nstars; // const size_t Ndm = Ntot > 0 ? Ntot - Nbaryons : 0; - /* Number of particles that we will write */ + 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 + * Recall that background particles are never inhibited and have no extras */ const size_t Ntot_written = e->s->nr_gparts - e->s->nr_inhibited_gparts - e->s->nr_extra_gparts; const size_t Ngas_written = @@ -754,11 +761,12 @@ void write_output_single(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; /* Format things in a Gadget-friendly array */ long long N_total[swift_type_count] = { - (long long)Ngas_written, (long long)Ndm_written, 0, 0, + (long long)Ngas_written, (long long)Ndm_written, + (long long)Ndm_background, 0, (long long)Nstars_written, (long long)Nblackholes_written}; /* File name */ @@ -1058,7 +1066,7 @@ void write_output_single(struct engine* e, const char* baseName, case swift_type_dark_matter: { if (Ntot == Ndm_written) { - /* This is a DM-only run without inhibited particles */ + /* This is a DM-only run without background or inhibited particles */ N = Ntot; darkmatter_write_particles(gparts, list, &num_fields); if (with_fof) { @@ -1106,6 +1114,43 @@ void write_output_single(struct engine* e, const char* baseName, } } break; + case swift_type_dark_matter_background: { + + /* Ok, we need to fish out the particles we want */ + N = Ndm_background; + + /* Allocate temporary array */ + if (swift_memalign("gparts_written", (void**)&gparts_written, + gpart_align, + Ndm_background * sizeof(struct gpart)) != 0) + error("Error while allocating temporart memory for gparts"); + + if (with_stf) { + if (swift_memalign( + "gpart_group_written", (void**)&gpart_group_data_written, + gpart_align, + Ndm_background * sizeof(struct velociraptor_gpart_data)) != 0) + error( + "Error while allocating temporart memory for gparts STF " + "data"); + } + + /* Collect the non-inhibited DM particles from gpart */ + io_collect_gparts_background_to_write( + gparts, e->s->gpart_group_data, gparts_written, + gpart_group_data_written, Ntot, Ndm_background, with_stf); + + /* Select the fields to write */ + darkmatter_write_particles(gparts_written, list, &num_fields); + if (with_fof) { + num_fields += fof_write_gparts(gparts_written, list + num_fields); + } + if (with_stf) { + num_fields += velociraptor_write_gparts(gpart_group_data_written, + list + num_fields); + } + } break; + case swift_type_stars: { if (Nstars == Nstars_written) { diff --git a/src/space.c b/src/space.c index b25750b0d648a5123454a8d2551fb8999c6fa4cb..61ecf3fd81f4a57d274ad911e0471337625f944b 100644 --- a/src/space.c +++ b/src/space.c @@ -4532,6 +4532,7 @@ void space_convert_quantities(struct space *s, int verbose) { * @param hydro flag whether we are doing hydro or not? * @param self_gravity flag whether we are doing gravity or not? * @param star_formation flag whether we are doing star formation or not? + * @param DM_background Are we running with some DM background particles? * @param verbose Print messages to stdout or not. * @param dry_run If 1, just initialise stuff, don't do anything with the parts. * @@ -4546,7 +4547,8 @@ void space_init(struct space *s, struct swift_params *params, struct bpart *bparts, size_t Npart, size_t Ngpart, size_t Nspart, size_t Nbpart, int periodic, int replicate, int generate_gas_in_ics, int hydro, int self_gravity, - int star_formation, int verbose, int dry_run) { + int star_formation, int DM_background, int verbose, + int dry_run) { /* Clean-up everything */ bzero(s, sizeof(struct space)); @@ -4559,6 +4561,7 @@ void space_init(struct space *s, struct swift_params *params, s->with_self_gravity = self_gravity; s->with_hydro = hydro; s->with_star_formation = star_formation; + s->with_DM_background = DM_background; s->nr_parts = Npart; s->nr_gparts = Ngpart; s->nr_sparts = Nspart; diff --git a/src/space.h b/src/space.h index c294bfae3612c699345bd4a267c40b67cffc2bf1..94ba0654352c6c33f1115337073bb76d099a0354 100644 --- a/src/space.h +++ b/src/space.h @@ -103,6 +103,9 @@ struct space { /*! Are we doing star formation? */ int with_star_formation; + /*! Are we running with some DM background particles? */ + int with_DM_background; + /*! Width of the top-level cells. */ double width[3]; @@ -307,7 +310,8 @@ void space_init(struct space *s, struct swift_params *params, struct bpart *bparts, size_t Npart, size_t Ngpart, size_t Nspart, size_t Nbpart, int periodic, int replicate, int generate_gas_in_ics, int hydro, int gravity, - int star_formation, int verbose, int dry_run); + int star_formation, int DM_background, 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);