diff --git a/examples/main.c b/examples/main.c index 5127ab9ef169be30153ea2c033c17a23c8b5461f..e795e210ca84b87410e1503c063d30c9a114ae83 100644 --- a/examples/main.c +++ b/examples/main.c @@ -867,19 +867,19 @@ int main(int argc, char *argv[]) { #if defined(HAVE_HDF5) #if defined(WITH_MPI) #if defined(HAVE_PARALLEL_HDF5) - read_ic_parallel(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, myrank, nr_nodes, MPI_COMM_WORLD, - MPI_INFO_NULL, nr_threads, dry_run); + read_ic_parallel( + ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &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, myrank, + nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads, dry_run); #else - read_ic_serial(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, myrank, nr_nodes, MPI_COMM_WORLD, - MPI_INFO_NULL, nr_threads, dry_run); + read_ic_serial( + ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &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, myrank, + nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads, dry_run); #endif #else read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, diff --git a/src/serial_io.c b/src/serial_io.c index 0d367b1cf423d0471a0f6235fa3a1c6877e9f279..3b3872762f237c4922bbc30fd866599547b23557 100644 --- a/src/serial_io.c +++ b/src/serial_io.c @@ -492,12 +492,12 @@ void writeArray(const struct engine* e, hid_t grp, char* fileName, void read_ic_serial(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 mpi_rank, - int mpi_size, MPI_Comm comm, MPI_Info info, 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 mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info, + int n_threads, int dry_run) { hid_t h_file = 0, h_grp = 0; /* GADGET has only cubic boxes (in cosmological mode) */ @@ -510,11 +510,13 @@ void read_ic_serial(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; struct unit_system* ic_units = (struct unit_system*)malloc(sizeof(struct unit_system)); /* Initialise counters */ - *Ngas = 0, *Ngparts = 0, *Nstars = 0, *Nblackholes = 0; + *Ngas = 0, *Ngparts = 0, *Ngparts_background = 0, *Nstars = 0, + *Nblackholes = 0; /* First read some information about the content */ if (mpi_rank == 0) { @@ -676,10 +678,13 @@ 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_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"); @@ -740,6 +745,13 @@ void read_ic_serial(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; @@ -789,17 +801,23 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units, /* 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 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 black holes 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); } @@ -847,6 +865,7 @@ void write_output_serial(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); @@ -864,7 +883,13 @@ void write_output_serial(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 = @@ -876,7 +901,7 @@ void write_output_serial(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; /* File name */ char fileName[FILENAME_BUFFER_SIZE]; @@ -888,7 +913,8 @@ void write_output_serial(struct engine* e, const char* baseName, e->snapshot_output_count); /* 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}; @@ -1242,7 +1268,8 @@ void write_output_serial(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 + */ Nparticles = Ntot; darkmatter_write_particles(gparts, list, &num_fields); if (with_fof) { @@ -1293,6 +1320,46 @@ void write_output_serial(struct engine* e, const char* baseName, } } break; + case swift_type_dark_matter_background: { + + /* Ok, we need to fish out the particles we want */ + Nparticles = 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/serial_io.h b/src/serial_io.h index 7d994181b135e9758a9e3a52971366c2d871f627..6badf076b8373ce2b6d5b2cd3d91a5b2f57ed8ae 100644 --- a/src/serial_io.h +++ b/src/serial_io.h @@ -38,12 +38,12 @@ void read_ic_serial(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 mpi_rank, - int mpi_size, MPI_Comm comm, MPI_Info info, 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 mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info, + int n_threads, int dry_run); void write_output_serial(struct engine* e, const char* baseName, const struct unit_system* internal_units,