diff --git a/doc/RTD/source/FriendsOfFriends/stand_alone_fof.rst b/doc/RTD/source/FriendsOfFriends/stand_alone_fof.rst index 4d3a1b6eb838802d95df95e4c4c0dfe98adac2bc..d5b341409b1553dcfb4062765d83f7ba14e18570 100644 --- a/doc/RTD/source/FriendsOfFriends/stand_alone_fof.rst +++ b/doc/RTD/source/FriendsOfFriends/stand_alone_fof.rst @@ -21,7 +21,12 @@ command line: * ``--hydro``: Read and process the gas particles, * ``--stars``: Read and process the star particles, - * ``--black-holes``: Read and process the black hole particles. + * ``--black-holes``: Read and process the black hole particles, + * ``--cosmology``: Consider cosmological terms. + +Running with cosmology is necessary when using a linking length based +on the mean inter-particle separation. For analysis where an absolute +linking length is used this is not necessary. The group catalogues contain a unique identifier for each group as well as the number of particles in each group and their total mass (in diff --git a/examples/main.c b/examples/main.c index a6d3ab1a3823a0754561cd6cc4b71e8c2c9e7fba..a24714e9090f66c25bd55b8b48c670279989139e 100644 --- a/examples/main.c +++ b/examples/main.c @@ -991,7 +991,8 @@ int main(int argc, char *argv[]) { /* Initialise the FOF properties */ bzero(&fof_properties, sizeof(struct fof_props)); #ifdef WITH_FOF - if (with_fof) fof_init(&fof_properties, params, &prog_const, &us); + if (with_fof) + fof_init(&fof_properties, params, &prog_const, &us, /*stand-alone=*/0); #endif /* Be verbose about what happens next */ diff --git a/examples/main_fof.c b/examples/main_fof.c index 42a304abee2887f773aa54579479419ccf30628b..25eef7b6e5bd4795566732a58435705318f5c9c6 100644 --- a/examples/main_fof.c +++ b/examples/main_fof.c @@ -139,8 +139,8 @@ int main(int argc, char *argv[]) { int with_aff = 0; int dump_tasks = 0; int dump_threadpool = 0; - int with_fof = 1; int with_fp_exceptions = 0; + int with_cosmology = 0; int with_stars = 0; int with_black_holes = 0; int with_hydro = 0; @@ -160,6 +160,8 @@ int main(int argc, char *argv[]) { OPT_HELP(), OPT_GROUP(" Friends-of-Friends options:\n"), + OPT_BOOLEAN('c', "cosmology", &with_cosmology, + "Run with cosmological information.", NULL, 0, 0), OPT_BOOLEAN(0, "hydro", &with_hydro, "Read gas particles from the ICs.", NULL, 0, 0), OPT_BOOLEAN(0, "stars", &with_stars, "Read stars from the ICs.", NULL, 0, @@ -257,8 +259,7 @@ int main(int argc, char *argv[]) { /* Write output parameter file */ if (myrank == 0 && output_parameters_filename != NULL) { - io_write_output_field_parameter(output_parameters_filename, - /*with_cosmology=*/1); + io_write_output_field_parameter(output_parameters_filename, with_cosmology); printf("End of run.\n"); return 0; } @@ -389,6 +390,10 @@ int main(int argc, char *argv[]) { } #endif + /* Prepare and verify the selection of outputs */ + io_prepare_output_fields(output_options, with_cosmology, /*with_fof=*/1, + /*with_structure_finding=*/0); + /* Initialize unit system and constants */ units_init_from_params(&us, params, "InternalUnitSystem"); phys_const_init(&us, params, &prog_const); @@ -414,12 +419,23 @@ int main(int argc, char *argv[]) { params, "InitialConditions:cleanup_h_factors", 0); const int cleanup_sqrt_a = parser_get_opt_param_int( params, "InitialConditions:cleanup_velocity_factors", 0); + /* Initialise the cosmology */ - cosmology_init(params, &us, &prog_const, &cosmo); + if (with_cosmology) + cosmology_init(params, &us, &prog_const, &cosmo); + else + cosmology_init_no_cosmo(&cosmo); + if (myrank == 0 && with_cosmology) cosmology_print(&cosmo); + + /* Initialise the equation of state */ + if (with_hydro) + eos_init(&eos, &prog_const, &us, params); + else + bzero(&eos, sizeof(struct eos_parameters)); /* Initialise the FOF properties */ bzero(&fof_properties, sizeof(struct fof_props)); - if (with_fof) fof_init(&fof_properties, params, &prog_const, &us); + fof_init(&fof_properties, params, &prog_const, &us, /*stand-alone=*/1); /* Be verbose about what happens next */ if (myrank == 0) message("Reading ICs from file '%s'", ICfileName); @@ -438,27 +454,27 @@ 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, - &Ngpart_background, &Nspart, &Nbpart, &flag_entropy_ICs, with_hydro, - /*with_grav=*/1, with_stars, with_black_holes, /*with_cosmology=*/1, - cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank, nr_nodes, - MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads, - /*dry_run=*/0); + read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, + &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart, + &flag_entropy_ICs, with_hydro, + /*with_grav=*/1, with_stars, with_black_holes, + with_cosmology, cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, + myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads, + /*dry_run=*/0); #else - read_ic_serial( - ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas, &Ngpart, - &Ngpart_background, &Nspart, &Nbpart, &flag_entropy_ICs, with_hydro, - /*with_grav=*/1, with_stars, with_black_holes, /*with_cosmology=*/1, - cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank, nr_nodes, - MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads, /*dry_run=*/0); + read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas, + &Ngpart, &Ngpart_background, &Nspart, &Nbpart, + &flag_entropy_ICs, with_hydro, + /*with_grav=*/1, with_stars, with_black_holes, with_cosmology, + cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank, nr_nodes, + MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads, /*dry_run=*/0); #endif #else read_ic_single( ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart, &flag_entropy_ICs, with_hydro, - /*with_grav=*/1, with_stars, with_black_holes, /*with_cosmology=*/1, - cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads, /*dry_run=*/0); + /*with_grav=*/1, with_stars, with_black_holes, with_cosmology, cleanup_h, + cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads, /*dry_run=*/0); #endif #endif if (myrank == 0) { @@ -532,7 +548,7 @@ int main(int argc, char *argv[]) { /* Initialise the gravity scheme */ bzero(&gravity_properties, sizeof(struct gravity_props)); gravity_props_init(&gravity_properties, params, &prog_const, &cosmo, - /*with_cosmology=*/1, /*with_external_gravity=*/0, + with_cosmology, /*with_external_gravity=*/0, with_baryon_particles, with_DM_particles, with_DM_background_particles, periodic, s.dim); @@ -583,8 +599,8 @@ int main(int argc, char *argv[]) { /* Construct the engine policy */ int engine_policies = ENGINE_POLICY | engine_policy_steal; engine_policies |= engine_policy_self_gravity; - engine_policies |= engine_policy_cosmology; engine_policies |= engine_policy_fof; + if (with_cosmology) engine_policies |= engine_policy_cosmology; /* Initialize the engine with the space and policies. */ if (myrank == 0) clocks_gettime(&tic); diff --git a/src/common_io.c b/src/common_io.c index 194bf334870db682adddff86d42907fcc9aa2705..457016bb700bd0d31d6310e0e9a4bd80eeed1aac 100644 --- a/src/common_io.c +++ b/src/common_io.c @@ -206,11 +206,13 @@ void io_assert_valid_header_cosmology(hid_t h_grp, double a) { &redshift_from_snapshot); /* If the Header/Redshift value is not present, then we skip this check */ - if (redshift_from_snapshot == -1.0) { - return; - } + if (redshift_from_snapshot == -1.0) return; const double current_redshift = 1.0 / a - 1.0; + + /* Escape if non-cosmological */ + if (current_redshift == 0.) return; + const double redshift_fractional_difference = fabs(redshift_from_snapshot - current_redshift) / current_redshift; diff --git a/src/fof.c b/src/fof.c index ff15262bab336f336af47d3cc2e3710b1bee355a..58c9268d67ff5b20541e5c05999dfb01f3a4e17d 100644 --- a/src/fof.c +++ b/src/fof.c @@ -88,10 +88,12 @@ static integertime_t ti_current; * @param params the parameter file parser. * @param phys_const The physical constants in internal units. * @param us The internal unit system. + * @param stand_alone_fof Are we initialising for stand-alone? (1) or + * on-the-fly? (0) */ void fof_init(struct fof_props *props, struct swift_params *params, - const struct phys_const *phys_const, - const struct unit_system *us) { + const struct phys_const *phys_const, const struct unit_system *us, + const int stand_alone_fof) { /* Base name for the FOF output file */ parser_get_param_string(params, "FOF:basename", props->base_name); @@ -118,24 +120,27 @@ void fof_init(struct fof_props *props, struct swift_params *params, /* Read the linking length ratio to the mean inter-particle separation. */ props->l_x_ratio = - parser_get_param_double(params, "FOF:linking_length_ratio"); - - if (props->l_x_ratio <= 0.) - error("The FOF linking length can't be negative!"); + parser_get_opt_param_double(params, "FOF:linking_length_ratio", -1.); /* Read value of absolute linking length aksed by the user */ props->l_x_absolute = parser_get_opt_param_double(params, "FOF:absolute_linking_length", -1.); - if (props->l_x_ratio <= 0. && props->l_x_ratio != -1.) + if (props->l_x_ratio == -1. && props->l_x_absolute <= 0.) error("The FOF linking length can't be negative!"); - /* Read the minimal halo mass for black hole seeding */ - props->seed_halo_mass = - parser_get_param_double(params, "FOF:black_hole_seed_halo_mass_Msun"); + if (props->l_x_ratio <= 0. && props->l_x_absolute == -1.) + error("The FOF linking length ratio can't be negative!"); - /* Convert to internal units */ - props->seed_halo_mass *= phys_const->const_solar_mass; + if (!stand_alone_fof) { + + /* Read the minimal halo mass for black hole seeding */ + props->seed_halo_mass = + parser_get_param_double(params, "FOF:black_hole_seed_halo_mass_Msun"); + + /* Convert to internal units */ + props->seed_halo_mass *= phys_const->const_solar_mass; + } #if defined(WITH_MPI) && defined(UNION_BY_SIZE_OVER_MPI) if (engine_rank == 0) @@ -262,26 +267,6 @@ void fof_allocate(const struct space *s, const long long total_nr_DM_particles, } } - /* Calculate the mean inter-particle separation as if we were in - a scenario where the entire box was filled with high-resolution - particles */ - const double Omega_m = s->e->cosmology->Omega_m; - const double Omega_b = s->e->cosmology->Omega_b; - const double critical_density_0 = s->e->cosmology->critical_density_0; - double mean_matter_density; - if (s->with_hydro) - mean_matter_density = (Omega_m - Omega_b) * critical_density_0; - else - mean_matter_density = Omega_m * critical_density_0; - - /* Mean inter-particle separation of the DM particles */ - const double mean_inter_particle_sep = - cbrt(high_res_DM_mass / mean_matter_density); - - /* Calculate the particle linking length based upon the mean inter-particle - * spacing of the DM particles. */ - const double l_x = props->l_x_ratio * mean_inter_particle_sep; - /* Are we using the aboslute value or the one derived from the mean inter-particle sepration? */ if (props->l_x_absolute != -1.) { @@ -292,6 +277,33 @@ void fof_allocate(const struct space *s, const long long total_nr_DM_particles, props->l_x_absolute); } } else { + + /* Safety check */ + if (!(s->e->policy | engine_policy_cosmology)) + error( + "Attempting to run FoF on a simulation using cosmological " + "information but cosmology was not initialised"); + + /* Calculate the mean inter-particle separation as if we were in + a scenario where the entire box was filled with high-resolution + particles */ + const double Omega_m = s->e->cosmology->Omega_m; + const double Omega_b = s->e->cosmology->Omega_b; + const double critical_density_0 = s->e->cosmology->critical_density_0; + double mean_matter_density; + if (s->with_hydro) + mean_matter_density = (Omega_m - Omega_b) * critical_density_0; + else + mean_matter_density = Omega_m * critical_density_0; + + /* Mean inter-particle separation of the DM particles */ + const double mean_inter_particle_sep = + cbrt(high_res_DM_mass / mean_matter_density); + + /* Calculate the particle linking length based upon the mean inter-particle + * spacing of the DM particles. */ + const double l_x = props->l_x_ratio * mean_inter_particle_sep; + props->l_x2 = l_x * l_x; if (s->e->nodeID == 0) { diff --git a/src/fof.h b/src/fof.h index 160f7c0af8f7f6f7a589c6978cd107b7ffbff1fc..8b54534b1d8836cd0444b01a9d994f63b6ac438b 100644 --- a/src/fof.h +++ b/src/fof.h @@ -159,8 +159,8 @@ struct cell_pair_indices { /* Function prototypes. */ void fof_init(struct fof_props *props, struct swift_params *params, - const struct phys_const *phys_const, - const struct unit_system *us); + const struct phys_const *phys_const, const struct unit_system *us, + const int stand_alone_fof); void fof_create_mpi_types(void); void fof_allocate(const struct space *s, const long long total_nr_DM_particles, struct fof_props *props);