diff --git a/doc/RTD/source/InitialConditions/index.rst b/doc/RTD/source/InitialConditions/index.rst
index 52d79ff8ddf12cf21cc08a9e0fa3eba69a8ca78d..66fb9f01f1189c8d9f96ed2f1399e0c5001a8bc3 100644
--- a/doc/RTD/source/InitialConditions/index.rst
+++ b/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
diff --git a/examples/main.c b/examples/main.c
index 62e5e85ba13527f1119b73234e2ceb1122f9f397..35f20ca1331b02c5e3e93796446a2f8d75b985cc 100644
--- a/examples/main.c
+++ b/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.",
diff --git a/examples/main_fof.c b/examples/main_fof.c
index c96cb1267e8c3ab6890de121a9d7fb3fb16443df..e4b13187435e1eb600cae4ce05e9332650c4fddb 100644
--- a/examples/main_fof.c
+++ b/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) {
diff --git a/src/common_io.c b/src/common_io.c
index 6f090ee01d1c870408223f90dba0770b5f76f567..63faace19fbb388dc85d00b603499518039bd66f 100644
--- a/src/common_io.c
+++ b/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.
*
diff --git a/src/common_io.h b/src/common_io.h
index aa4b0c552f8a416a4aa42ecce5e9ae0b74dbc731..b4a846dc4f683cb4e18e3bba03c4cc95235a8854 100644
--- a/src/common_io.h
+++ b/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);
diff --git a/src/const.h b/src/const.h
index d7017699cff56b1de545e6e0ae78013f000b76d0..521d094f51e61fde1e1f983caa7194bc0242214b 100644
--- a/src/const.h
+++ b/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
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 8ce21085e0dc0d1d126ed510262e14842ddf4019..ccba33d07500f4e22e365942622f6392cfbb0166 100644
--- a/src/parallel_io.c
+++ b/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);
diff --git a/src/parallel_io.h b/src/parallel_io.h
index 465835013847347a907b3f81fd3611992070d0d7..6f1b397ef337c05da99a8017cc562212067e32a5 100644
--- a/src/parallel_io.h
+++ b/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,
diff --git a/src/serial_io.c b/src/serial_io.c
index ca832a4beb16f8f39760c3c0a2950ecd003cc02b..0035ba89efb73e0e6d585bbc4b1e180aae507530 100644
--- a/src/serial_io.c
+++ b/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);
diff --git a/src/serial_io.h b/src/serial_io.h
index c6ca7ad66df038791fbd62277f406a7f34340e04..5b51b9b41870846c40b8d3d358f31ebfe7202203 100644
--- a/src/serial_io.h
+++ b/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,
diff --git a/src/single_io.c b/src/single_io.c
index 017f9498ffc7d11e1a03d27eca733d2d417b9777..ad2560db96f7bf2b4c60141bf9ff945b2bb745b0 100644
--- a/src/single_io.c
+++ b/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);
diff --git a/src/single_io.h b/src/single_io.h
index 47b575272c2fbe5b170eb4370743cd6830b99ae0..591f0ef8be6ed831d5cf2a6d465498cc259286d5 100644
--- a/src/single_io.h
+++ b/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,