Redshift IC Check

doc/RTD/source/InitialConditions/index.rst

@@ -124,7 +124,13 @@ GADGET-2 based analysis programs:
   value different from 1, the code will return an error message.
 + ``Time``, time of the start of the simulation in internal units or expressed
   as a scale-factor for cosmological runs. **SWIFT ignores this and reads it
-  from the parameter file**, behaviour that matches the GADGET-2 code.
+  from the parameter file**, behaviour that matches the GADGET-2 code.  Note
+  that SWIFT writes the current time since the Big Bang, not scale-factor, to
+  this variable in snapshots.
++ ``Redshift``, the redshift at the start of the simulation. SWIFT checks this
+  (if present) against ``a_begin`` in the parameter file at the start of
+  cosmological runs.  Note that we explicitly do **not** compare the ``Time``
+  variable due to its ambiguous meaning.
 
 
 Particle Data

examples/main.c

@@ -869,25 +869,26 @@ int main(int argc, char *argv[]) {
     read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
                      &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
                      &flag_entropy_ICs, with_hydro, with_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);
+                     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);
 #else
     read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
                    &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
                    &flag_entropy_ICs, with_hydro, with_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);
+                   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);
 #endif
 #else
     read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
                    &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
                    &flag_entropy_ICs, with_hydro, with_gravity, with_stars,
-                   with_black_holes, cleanup_h, cleanup_sqrt_a, cosmo.h,
-                   cosmo.a, nr_threads, dry_run);
+                   with_black_holes, with_cosmology, cleanup_h, cleanup_sqrt_a,
+                   cosmo.h, cosmo.a, nr_threads, dry_run);
 #endif
 #endif
+
     if (myrank == 0) {
       clocks_gettime(&toc);
       message("Reading initial conditions took %.3f %s.",

examples/main_fof.c

@@ -430,27 +430,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, 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=*/1,
+      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, 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=*/1,
+      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, cleanup_h,
-                 cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads, /*dry_run=*/0);
+  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);
 #endif
 #endif
   if (myrank == 0) {

src/common_io.c

@@ -30,6 +30,7 @@
 /* Local includes. */
 #include "black_holes_io.h"
 #include "chemistry_io.h"
+#include "const.h"
 #include "cooling_io.h"
 #include "error.h"
 #include "fof_io.h"
@@ -138,6 +139,85 @@ void io_read_attribute(hid_t grp, const char* name, enum IO_DATA_TYPE type,
   H5Aclose(h_attr);
 }
 
+/**
+ * @brief Reads an attribute from a given HDF5 group.
+ *
+ * @param grp The group from which to read.
+ * @param name The name of the attribute to read.
+ * @param type The #IO_DATA_TYPE of the attribute.
+ * @param data (output) The attribute read from the HDF5 group.
+ *
+ * Exits gracefully (i.e. does not read the attribute at all) if
+ * it is not present, unless debugging checks are activated. If they are,
+ * and the read fails, we print a warning.
+ */
+void io_read_attribute_graceful(hid_t grp, const char* name,
+                                enum IO_DATA_TYPE type, void* data) {
+
+  /* First, we need to check if this attribute exists to avoid raising errors
+   * within the HDF5 library if we attempt to access an attribute that does
+   * not exist. */
+  const htri_t h_exists = H5Aexists(grp, name);
+
+  if (h_exists <= 0) {
+  /* Attribute either does not exist (0) or function failed (-ve) */
+#ifdef SWIFT_DEBUG_CHECKS
+    message("WARNING: attribute '%s' does not exist.", name);
+#endif
+  } else {
+    /* Ok, now we know that it exists we can read it. */
+    const hid_t h_attr = H5Aopen(grp, name, H5P_DEFAULT);
+
+    if (h_attr >= 0) {
+      const hid_t h_err = H5Aread(h_attr, io_hdf5_type(type), data);
+      if (h_err < 0) {
+      /* Explicitly do nothing unless debugging checks are activated */
+#ifdef SWIFT_DEBUG_CHECKS
+        message("WARNING: unable to read attribute '%s'", name);
+#endif
+      }
+    } else {
+#ifdef SWIFT_DEBUG_CHECKS
+      if (h_attr < 0) {
+        message("WARNING: was unable to open attribute '%s'", name);
+      }
+#endif
+    }
+
+    H5Aclose(h_attr);
+  }
+}
+
+/**
+ * @brief Asserts that the redshift in the initial conditions and the one
+ *        specified by the parameter file match.
+ *
+ * @param h_grp The Header group from the ICs
+ * @param a Current scale factor as specified by parameter file
+ */
+void io_assert_valid_header_cosmology(hid_t h_grp, double a) {
+
+  double redshift_from_snapshot = -1.0;
+  io_read_attribute_graceful(h_grp, "Redshift", DOUBLE,
+                             &redshift_from_snapshot);
+
+  /* If the Header/Redshift value is not present, then we skip this check */
+  if (redshift_from_snapshot == -1.0) {
+    return;
+  }
+
+  const double current_redshift = 1.0 / a - 1.0;
+  const double redshift_fractional_difference =
+      fabs(redshift_from_snapshot - current_redshift) / current_redshift;
+
+  if (redshift_fractional_difference >= io_redshift_tolerance) {
+    error(
+        "Initial redshift specified in parameter file (%lf) and redshift "
+        "read from initial conditions (%lf) are inconsistent.",
+        current_redshift, redshift_from_snapshot);
+  }
+}
+
 /**
  * @brief Write an attribute to a given HDF5 group.
  *

src/common_io.h

@@ -68,6 +68,9 @@ hid_t io_hdf5_type(enum IO_DATA_TYPE type);
 
 void io_read_attribute(hid_t grp, const char* name, enum IO_DATA_TYPE type,
                        void* data);
+void io_read_attribute_graceful(hid_t grp, const char* name,
+                                enum IO_DATA_TYPE type, void* data);
+void io_assert_valid_header_cosmology(hid_t h_grp, double a);
 
 void io_write_attribute(hid_t grp, const char* name, enum IO_DATA_TYPE type,
                         const void* data, int num);

src/const.h

@@ -26,6 +26,11 @@
 /* Time-step limiter maximal difference in signal velocity */
 #define const_limiter_max_v_sig_ratio 4.1f
 
+/* I/O Constant; this determines the relative tolerance between the value of
+ * redshift read from the snapshot, and the value from the parameter file. This
+ * current value asserts that they must match within 0.1%. */
+#define io_redshift_tolerance 1e-3f
+
 /* Type of gradients to use (GIZMO_SPH only) */
 /* If no option is chosen, no gradients are used (first order scheme) */
 //#define GRADIENTS_SPH

src/parallel_io.c

@@ -688,6 +688,7 @@ void writeArray(struct engine* e, hid_t grp, char* fileName,
  * @param with_gravity Are we running with gravity ?
  * @param with_stars Are we running with stars ?
  * @param with_black_holes Are we running with black holes ?
+ * @param with_cosmology Are we running with cosmology ?
  * @param cleanup_h Are we cleaning-up h-factors from the quantities we read?
  * @param cleanup_sqrt_a Are we cleaning-up the sqrt(a) factors in the Gadget
  * IC velocities?
@@ -707,9 +708,9 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
                       size_t* Ngas, 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 with_black_holes, int with_cosmology, 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;
@@ -773,6 +774,13 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
   io_read_attribute(h_grp, "NumPart_Total_HighWord", LONGLONG,
                     numParticles_highWord);
 
+  /* Check that the user is not doing something silly when they e.g. restart
+   * from a snapshot by asserting that the current scale-factor (from
+   * parameter file) and the redshift in the header are consistent */
+  if (with_cosmology) {
+    io_assert_valid_header_cosmology(h_grp, a);
+  }
+
   for (int ptype = 0; ptype < swift_type_count; ++ptype)
     N_total[ptype] =
         (numParticles[ptype]) + (numParticles_highWord[ptype] << 32);

src/parallel_io.h

@@ -39,10 +39,10 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
                       size_t* Ngas, size_t* Ngparts, size_t* Ngparts_background,
                       size_t* Nsparts, size_t* Nbparts, 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 nr_threads,
-                      int dry_run);
+                      int with_black_holes, int with_cosmology, int cleanup_h,
+                      int cleanup_sqrt_a, double h, double a, int mpi_rank,
+                      int mpi_size, MPI_Comm comm, MPI_Info info,
+                      int nr_threads, int dry_run);
 
 void write_output_parallel(struct engine* e, const char* baseName,
                            const struct unit_system* internal_units,

src/serial_io.c

@@ -471,6 +471,7 @@ void writeArray(const struct engine* e, hid_t grp, char* fileName,
  * @param with_gravity Are we reading/creating #gpart arrays ?
  * @param with_stars Are we reading star particles ?
  * @param with_black_holes Are we reading black hole particles ?
+ * @param with_cosmology Are we running with cosmology ?
  * @param cleanup_h Are we cleaning-up h-factors from the quantities we read?
  * @param cleanup_sqrt_a Are we cleaning-up the sqrt(a) factors in the Gadget
  * IC velocities?
@@ -497,9 +498,9 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
                     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) {
+                    int with_cosmology, 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) */
@@ -569,6 +570,13 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
     io_read_attribute(h_grp, "NumPart_Total_HighWord", LONGLONG,
                       numParticles_highWord);
 
+    /* Check that the user is not doing something silly when they e.g. restart
+     * from a snapshot by asserting that the current scale-factor (from
+     * parameter file) and the redshift in the header are consistent */
+    if (with_cosmology) {
+      io_assert_valid_header_cosmology(h_grp, a);
+    }
+
     for (int ptype = 0; ptype < swift_type_count; ++ptype)
       N_total[ptype] =
           (numParticles[ptype]) + (numParticles_highWord[ptype] << 32);

src/serial_io.h

@@ -41,9 +41,9 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
                     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);
+                    int with_cosmology, 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,

src/single_io.c

@@ -385,6 +385,7 @@ void writeArray(const struct engine* e, hid_t grp, char* fileName,
  * @param with_gravity Are we reading/creating #gpart arrays ?
  * @param with_stars Are we reading star particles ?
  * @param with_black_hole Are we reading black hole particles ?
+ * @param with_cosmology Are we running with cosmology ?
  * @param cleanup_h Are we cleaning-up h-factors from the quantities we read?
  * @param cleanup_sqrt_a Are we cleaning-up the sqrt(a) factors in the Gadget
  * IC velocities?
@@ -407,8 +408,8 @@ void read_ic_single(const char* fileName,
                     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 n_threads, int dry_run) {
+                    int with_cosmology, int cleanup_h, int cleanup_sqrt_a,
+                    double h, double a, int n_threads, int dry_run) {
 
   hid_t h_file = 0, h_grp = 0;
   /* GADGET has only cubic boxes (in cosmological mode) */
@@ -468,6 +469,13 @@ void read_ic_single(const char* fileName,
   io_read_attribute(h_grp, "NumPart_Total_HighWord", LONGLONG,
                     numParticles_highWord);
 
+  /* Check that the user is not doing something silly when they e.g. restart
+   * from a snapshot by asserting that the current scale-factor (from
+   * parameter file) and the redshift in the header are consistent */
+  if (with_cosmology) {
+    io_assert_valid_header_cosmology(h_grp, a);
+  }
+
   for (int ptype = 0; ptype < swift_type_count; ++ptype)
     N[ptype] = (numParticles[ptype]) + (numParticles_highWord[ptype] << 32);
 

src/single_io.h

@@ -37,8 +37,8 @@ void read_ic_single(const char* fileName,
                     size_t* Ndm, size_t* Ndm_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 nr_threads, int dry_run);
+                    int with_cosmology, int cleanup_h, int cleanup_sqrt_a,
+                    double h, double a, int nr_threads, int dry_run);
 
 void write_output_single(struct engine* e, const char* baseName,
                          const struct unit_system* internal_units,