Initialise the output_options also in the stand-alone FoF code

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

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 */

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);

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;
 

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) {

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);