Commit 21144edc authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'fof_io_fix' into 'master'

Initialise the output_options also in the stand-alone FoF code

Closes #689

See merge request !1117
parents c577fae8 1ac1d920
...@@ -21,7 +21,12 @@ command line: ...@@ -21,7 +21,12 @@ command line:
* ``--hydro``: Read and process the gas particles, * ``--hydro``: Read and process the gas particles,
* ``--stars``: Read and process the star 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 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 well as the number of particles in each group and their total mass (in
......
...@@ -991,7 +991,8 @@ int main(int argc, char *argv[]) { ...@@ -991,7 +991,8 @@ int main(int argc, char *argv[]) {
/* Initialise the FOF properties */ /* Initialise the FOF properties */
bzero(&fof_properties, sizeof(struct fof_props)); bzero(&fof_properties, sizeof(struct fof_props));
#ifdef WITH_FOF #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 #endif
/* Be verbose about what happens next */ /* Be verbose about what happens next */
......
...@@ -139,8 +139,8 @@ int main(int argc, char *argv[]) { ...@@ -139,8 +139,8 @@ int main(int argc, char *argv[]) {
int with_aff = 0; int with_aff = 0;
int dump_tasks = 0; int dump_tasks = 0;
int dump_threadpool = 0; int dump_threadpool = 0;
int with_fof = 1;
int with_fp_exceptions = 0; int with_fp_exceptions = 0;
int with_cosmology = 0;
int with_stars = 0; int with_stars = 0;
int with_black_holes = 0; int with_black_holes = 0;
int with_hydro = 0; int with_hydro = 0;
...@@ -160,6 +160,8 @@ int main(int argc, char *argv[]) { ...@@ -160,6 +160,8 @@ int main(int argc, char *argv[]) {
OPT_HELP(), OPT_HELP(),
OPT_GROUP(" Friends-of-Friends options:\n"), 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.", OPT_BOOLEAN(0, "hydro", &with_hydro, "Read gas particles from the ICs.",
NULL, 0, 0), NULL, 0, 0),
OPT_BOOLEAN(0, "stars", &with_stars, "Read stars from the ICs.", NULL, 0, OPT_BOOLEAN(0, "stars", &with_stars, "Read stars from the ICs.", NULL, 0,
...@@ -257,8 +259,7 @@ int main(int argc, char *argv[]) { ...@@ -257,8 +259,7 @@ int main(int argc, char *argv[]) {
/* Write output parameter file */ /* Write output parameter file */
if (myrank == 0 && output_parameters_filename != NULL) { if (myrank == 0 && output_parameters_filename != NULL) {
io_write_output_field_parameter(output_parameters_filename, io_write_output_field_parameter(output_parameters_filename, with_cosmology);
/*with_cosmology=*/1);
printf("End of run.\n"); printf("End of run.\n");
return 0; return 0;
} }
...@@ -389,6 +390,10 @@ int main(int argc, char *argv[]) { ...@@ -389,6 +390,10 @@ int main(int argc, char *argv[]) {
} }
#endif #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 */ /* Initialize unit system and constants */
units_init_from_params(&us, params, "InternalUnitSystem"); units_init_from_params(&us, params, "InternalUnitSystem");
phys_const_init(&us, params, &prog_const); phys_const_init(&us, params, &prog_const);
...@@ -414,12 +419,23 @@ int main(int argc, char *argv[]) { ...@@ -414,12 +419,23 @@ int main(int argc, char *argv[]) {
params, "InitialConditions:cleanup_h_factors", 0); params, "InitialConditions:cleanup_h_factors", 0);
const int cleanup_sqrt_a = parser_get_opt_param_int( const int cleanup_sqrt_a = parser_get_opt_param_int(
params, "InitialConditions:cleanup_velocity_factors", 0); params, "InitialConditions:cleanup_velocity_factors", 0);
/* Initialise the cosmology */ /* Initialise the cosmology */
if (with_cosmology)
cosmology_init(params, &us, &prog_const, &cosmo); 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 */ /* Initialise the FOF properties */
bzero(&fof_properties, sizeof(struct fof_props)); 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 */ /* Be verbose about what happens next */
if (myrank == 0) message("Reading ICs from file '%s'", ICfileName); if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
...@@ -438,18 +454,18 @@ int main(int argc, char *argv[]) { ...@@ -438,18 +454,18 @@ int main(int argc, char *argv[]) {
#if defined(HAVE_HDF5) #if defined(HAVE_HDF5)
#if defined(WITH_MPI) #if defined(WITH_MPI)
#if defined(HAVE_PARALLEL_HDF5) #if defined(HAVE_PARALLEL_HDF5)
read_ic_parallel( read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas, &Ngpart, &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
&Ngpart_background, &Nspart, &Nbpart, &flag_entropy_ICs, with_hydro, &flag_entropy_ICs, with_hydro,
/*with_grav=*/1, with_stars, with_black_holes, /*with_cosmology=*/1, /*with_grav=*/1, with_stars, with_black_holes,
cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank, nr_nodes, with_cosmology, cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a,
MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
/*dry_run=*/0); /*dry_run=*/0);
#else #else
read_ic_serial( read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas,
ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas, &Ngpart, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
&Ngpart_background, &Nspart, &Nbpart, &flag_entropy_ICs, with_hydro, &flag_entropy_ICs, with_hydro,
/*with_grav=*/1, with_stars, with_black_holes, /*with_cosmology=*/1, /*with_grav=*/1, with_stars, with_black_holes, with_cosmology,
cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank, nr_nodes, cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank, nr_nodes,
MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads, /*dry_run=*/0); MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads, /*dry_run=*/0);
#endif #endif
...@@ -457,8 +473,8 @@ int main(int argc, char *argv[]) { ...@@ -457,8 +473,8 @@ int main(int argc, char *argv[]) {
read_ic_single( read_ic_single(
ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas, &Ngpart, ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas, &Ngpart,
&Ngpart_background, &Nspart, &Nbpart, &flag_entropy_ICs, with_hydro, &Ngpart_background, &Nspart, &Nbpart, &flag_entropy_ICs, with_hydro,
/*with_grav=*/1, with_stars, with_black_holes, /*with_cosmology=*/1, /*with_grav=*/1, with_stars, with_black_holes, with_cosmology, cleanup_h,
cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads, /*dry_run=*/0); cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads, /*dry_run=*/0);
#endif #endif
#endif #endif
if (myrank == 0) { if (myrank == 0) {
...@@ -532,7 +548,7 @@ int main(int argc, char *argv[]) { ...@@ -532,7 +548,7 @@ int main(int argc, char *argv[]) {
/* Initialise the gravity scheme */ /* Initialise the gravity scheme */
bzero(&gravity_properties, sizeof(struct gravity_props)); bzero(&gravity_properties, sizeof(struct gravity_props));
gravity_props_init(&gravity_properties, params, &prog_const, &cosmo, 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_baryon_particles, with_DM_particles,
with_DM_background_particles, periodic, s.dim); with_DM_background_particles, periodic, s.dim);
...@@ -583,8 +599,8 @@ int main(int argc, char *argv[]) { ...@@ -583,8 +599,8 @@ int main(int argc, char *argv[]) {
/* Construct the engine policy */ /* Construct the engine policy */
int engine_policies = ENGINE_POLICY | engine_policy_steal; int engine_policies = ENGINE_POLICY | engine_policy_steal;
engine_policies |= engine_policy_self_gravity; engine_policies |= engine_policy_self_gravity;
engine_policies |= engine_policy_cosmology;
engine_policies |= engine_policy_fof; engine_policies |= engine_policy_fof;
if (with_cosmology) engine_policies |= engine_policy_cosmology;
/* Initialize the engine with the space and policies. */ /* Initialize the engine with the space and policies. */
if (myrank == 0) clocks_gettime(&tic); if (myrank == 0) clocks_gettime(&tic);
......
...@@ -206,11 +206,13 @@ void io_assert_valid_header_cosmology(hid_t h_grp, double a) { ...@@ -206,11 +206,13 @@ void io_assert_valid_header_cosmology(hid_t h_grp, double a) {
&redshift_from_snapshot); &redshift_from_snapshot);
/* If the Header/Redshift value is not present, then we skip this check */ /* If the Header/Redshift value is not present, then we skip this check */
if (redshift_from_snapshot == -1.0) { if (redshift_from_snapshot == -1.0) return;
return;
}
const double current_redshift = 1.0 / a - 1.0; const double current_redshift = 1.0 / a - 1.0;
/* Escape if non-cosmological */
if (current_redshift == 0.) return;
const double redshift_fractional_difference = const double redshift_fractional_difference =
fabs(redshift_from_snapshot - current_redshift) / current_redshift; fabs(redshift_from_snapshot - current_redshift) / current_redshift;
......
...@@ -88,10 +88,12 @@ static integertime_t ti_current; ...@@ -88,10 +88,12 @@ static integertime_t ti_current;
* @param params the parameter file parser. * @param params the parameter file parser.
* @param phys_const The physical constants in internal units. * @param phys_const The physical constants in internal units.
* @param us The internal unit system. * @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, void fof_init(struct fof_props *props, struct swift_params *params,
const struct phys_const *phys_const, const struct phys_const *phys_const, const struct unit_system *us,
const struct unit_system *us) { const int stand_alone_fof) {
/* Base name for the FOF output file */ /* Base name for the FOF output file */
parser_get_param_string(params, "FOF:basename", props->base_name); 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, ...@@ -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. */ /* Read the linking length ratio to the mean inter-particle separation. */
props->l_x_ratio = props->l_x_ratio =
parser_get_param_double(params, "FOF:linking_length_ratio"); parser_get_opt_param_double(params, "FOF:linking_length_ratio", -1.);
if (props->l_x_ratio <= 0.)
error("The FOF linking length can't be negative!");
/* Read value of absolute linking length aksed by the user */ /* Read value of absolute linking length aksed by the user */
props->l_x_absolute = props->l_x_absolute =
parser_get_opt_param_double(params, "FOF:absolute_linking_length", -1.); 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!"); error("The FOF linking length can't be negative!");
if (props->l_x_ratio <= 0. && props->l_x_absolute == -1.)
error("The FOF linking length ratio can't be negative!");
if (!stand_alone_fof) {
/* Read the minimal halo mass for black hole seeding */ /* Read the minimal halo mass for black hole seeding */
props->seed_halo_mass = props->seed_halo_mass =
parser_get_param_double(params, "FOF:black_hole_seed_halo_mass_Msun"); parser_get_param_double(params, "FOF:black_hole_seed_halo_mass_Msun");
/* Convert to internal units */ /* Convert to internal units */
props->seed_halo_mass *= phys_const->const_solar_mass; props->seed_halo_mass *= phys_const->const_solar_mass;
}
#if defined(WITH_MPI) && defined(UNION_BY_SIZE_OVER_MPI) #if defined(WITH_MPI) && defined(UNION_BY_SIZE_OVER_MPI)
if (engine_rank == 0) if (engine_rank == 0)
...@@ -262,6 +267,23 @@ void fof_allocate(const struct space *s, const long long total_nr_DM_particles, ...@@ -262,6 +267,23 @@ void fof_allocate(const struct space *s, const long long total_nr_DM_particles,
} }
} }
/* Are we using the aboslute value or the one derived from the mean
inter-particle sepration? */
if (props->l_x_absolute != -1.) {
props->l_x2 = props->l_x_absolute * props->l_x_absolute;
if (s->e->nodeID == 0) {
message("Linking length is set to %e [internal units].",
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 /* Calculate the mean inter-particle separation as if we were in
a scenario where the entire box was filled with high-resolution a scenario where the entire box was filled with high-resolution
particles */ particles */
...@@ -282,16 +304,6 @@ void fof_allocate(const struct space *s, const long long total_nr_DM_particles, ...@@ -282,16 +304,6 @@ void fof_allocate(const struct space *s, const long long total_nr_DM_particles,
* spacing of the DM particles. */ * spacing of the DM particles. */
const double l_x = props->l_x_ratio * mean_inter_particle_sep; 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.) {
props->l_x2 = props->l_x_absolute * props->l_x_absolute;
if (s->e->nodeID == 0) {
message("Linking length is set to %e [internal units].",
props->l_x_absolute);
}
} else {
props->l_x2 = l_x * l_x; props->l_x2 = l_x * l_x;
if (s->e->nodeID == 0) { if (s->e->nodeID == 0) {
......
...@@ -159,8 +159,8 @@ struct cell_pair_indices { ...@@ -159,8 +159,8 @@ struct cell_pair_indices {
/* Function prototypes. */ /* Function prototypes. */
void fof_init(struct fof_props *props, struct swift_params *params, void fof_init(struct fof_props *props, struct swift_params *params,
const struct phys_const *phys_const, const struct phys_const *phys_const, const struct unit_system *us,
const struct unit_system *us); const int stand_alone_fof);
void fof_create_mpi_types(void); void fof_create_mpi_types(void);
void fof_allocate(const struct space *s, const long long total_nr_DM_particles, void fof_allocate(const struct space *s, const long long total_nr_DM_particles,
struct fof_props *props); struct fof_props *props);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment