Also apply the h-factor correction in the serial-io and parallel-io routines.

examples/main.c

@@ -639,14 +639,14 @@ int main(int argc, char *argv[]) {
     read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
                      &Ngpart, &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
                      (with_external_gravity || with_self_gravity), with_stars,
-                     myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL,
-                     nr_threads, dry_run);
+                     cleanup_h, cosmo.h, myrank, nr_nodes, MPI_COMM_WORLD,
+                     MPI_INFO_NULL, nr_threads, dry_run);
 #else
     read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
                    &Ngpart, &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
                    (with_external_gravity || with_self_gravity), with_stars,
-                   myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
-                   dry_run);
+                   cleanup_h, cosmo.h, myrank, nr_nodes, MPI_COMM_WORLD,
+                   MPI_INFO_NULL, nr_threads, dry_run);
 #endif
 #else
     read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,

src/parallel_io.c

@@ -69,11 +69,14 @@
  * @param offset Offset in the array where this mpi task starts writing.
  * @param internal_units The #unit_system used internally.
  * @param ic_units The #unit_system used in the snapshots.
+ * @param cleanup_h Are we removing h-factors from the ICs?
+ * @param h The value of the reduced Hubble constant to use for cleaning.
  */
 void readArray_chunk(hid_t h_data, hid_t h_plist_id,
                      const struct io_props props, size_t N, long long offset,
                      const struct unit_system* internal_units,
-                     const struct unit_system* ic_units) {
+                     const struct unit_system* ic_units, int cleanup_h,
+                     double h) {
 
   const size_t typeSize = io_sizeof_type(props.type);
   const size_t copySize = typeSize * props.dimension;
@@ -136,6 +139,23 @@ void readArray_chunk(hid_t h_data, hid_t h_plist_id,
     }
   }
 
+  /* Clean-up h if necessary */
+  const float h_factor_exp = units_h_factor(internal_units, props.units);
+  if (h_factor_exp != 0.f && exist != 0) {
+    const double h_factor = pow(h, h_factor_exp);
+
+    /* message("Multipltying '%s' by h^%f=%f", props.name, h_factor_exp,
+     * h_factor); */
+
+    if (io_is_double_precision(props.type)) {
+      double* temp_d = (double*)temp;
+      for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= h_factor;
+    } else {
+      float* temp_f = (float*)temp;
+      for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= h_factor;
+    }
+  }
+
   /* Copy temporary buffer to particle data */
   char* temp_c = (char*)temp;
   for (size_t i = 0; i < N; ++i)
@@ -158,11 +178,13 @@ void readArray_chunk(hid_t h_data, hid_t h_plist_id,
  * @param offset The offset in the array on disk for this rank.
  * @param internal_units The #unit_system used internally.
  * @param ic_units The #unit_system used in the ICs.
+ * @param cleanup_h Are we removing h-factors from the ICs?
+ * @param h The value of the reduced Hubble constant to use for cleaning.
  */
 void readArray(hid_t grp, struct io_props props, size_t N, long long N_total,
                int mpi_rank, long long offset,
                const struct unit_system* internal_units,
-               const struct unit_system* ic_units) {
+               const struct unit_system* ic_units, int cleanup_h, double h) {
 
   const size_t typeSize = io_sizeof_type(props.type);
   const size_t copySize = typeSize * props.dimension;
@@ -541,6 +563,8 @@ void writeArray(struct engine* e, hid_t grp, char* fileName,
  * @param with_hydro Are we running with hydro ?
  * @param with_gravity Are we running with gravity ?
  * @param with_stars Are we running with stars ?
+ * @param cleanup_h Are we cleaning-up h-factors from the quantities we read?
+ * @param h The value of the reduced Hubble constant to use for correction.
  * @param mpi_rank The MPI rank of this node
  * @param mpi_size The number of MPI ranks
  * @param comm The MPI communicator
@@ -554,8 +578,9 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
                       struct spart** sparts, size_t* Ngas, size_t* Ngparts,
                       size_t* Nstars, int* periodic, int* flag_entropy,
                       int with_hydro, int with_gravity, int with_stars,
-                      int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
-                      int n_threads, int dry_run) {
+                      int cleanup_h, double h, 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) */
@@ -626,6 +651,13 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
   else if (hydro_dimension == 1)
     dim[2] = dim[1] = dim[0];
 
+  /* Convert the box size if we want to clean-up h-factors */
+  if (cleanup_h) {
+    dim[0] /= h;
+    dim[1] /= h;
+    dim[2] /= h;
+  }
+
   /* message("Found %lld particles in a %speriodic box of size [%f %f %f].", */
   /* 	  N_total[0], (periodic ? "": "non-"), dim[0], */
   /* 	  dim[1], dim[2]); */
@@ -771,7 +803,7 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
     if (!dry_run)
       for (int i = 0; i < num_fields; ++i)
         readArray(h_grp, list[i], Nparticles, N_total[ptype], mpi_rank,
-                  offset[ptype], internal_units, ic_units);
+                  offset[ptype], internal_units, ic_units, cleanup_h, h);
 
     /* Close particle group */
     H5Gclose(h_grp);

src/parallel_io.h

@@ -22,6 +22,8 @@
 /* Config parameters. */
 #include "../config.h"
 
+#if defined(HAVE_HDF5) && defined(WITH_MPI) && defined(HAVE_PARALLEL_HDF5)
+
 /* MPI headers. */
 #ifdef WITH_MPI
 #include <mpi.h>
@@ -32,15 +34,14 @@
 #include "part.h"
 #include "units.h"
 
-#if defined(HAVE_HDF5) && defined(WITH_MPI) && defined(HAVE_PARALLEL_HDF5)
-
 void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
                       double dim[3], struct part** parts, struct gpart** gparts,
                       struct spart** sparts, size_t* Ngas, size_t* Ngparts,
                       size_t* Nsparts, int* periodic, int* flag_entropy,
                       int with_hydro, int with_gravity, int with_stars,
-                      int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
-                      int nr_threads, int dry_run);
+                      int cleanup_h, double h, 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

@@ -67,6 +67,8 @@
  * @param offset The offset position where this rank starts reading.
  * @param internal_units The #unit_system used internally
  * @param ic_units The #unit_system used in the ICs
+ * @param cleanup_h Are we removing h-factors from the ICs?
+ * @param h The value of the reduced Hubble constant to use for cleaning.
  *
  * @todo A better version using HDF5 hyper-slabs to read the file directly into
  * the part array will be written once the structures have been stabilized.
@@ -74,7 +76,7 @@
 void readArray(hid_t grp, const struct io_props props, size_t N,
                long long N_total, long long offset,
                const struct unit_system* internal_units,
-               const struct unit_system* ic_units) {
+               const struct unit_system* ic_units, int cleanup_h, double h) {
 
   const size_t typeSize = io_sizeof_type(props.type);
   const size_t copySize = typeSize * props.dimension;
@@ -161,6 +163,23 @@ void readArray(hid_t grp, const struct io_props props, size_t N,
     }
   }
 
+  /* Clean-up h if necessary */
+  const float h_factor_exp = units_h_factor(internal_units, props.units);
+  if (h_factor_exp != 0.f && exist != 0) {
+    const double h_factor = pow(h, h_factor_exp);
+
+    /* message("Multipltying '%s' by h^%f=%f", props.name, h_factor_exp,
+     * h_factor); */
+
+    if (io_is_double_precision(props.type)) {
+      double* temp_d = (double*)temp;
+      for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= h_factor;
+    } else {
+      float* temp_f = (float*)temp;
+      for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= h_factor;
+    }
+  }
+
   /* Copy temporary buffer to particle data */
   char* temp_c = temp;
   for (size_t i = 0; i < N; ++i)
@@ -382,6 +401,8 @@ void writeArray(const struct engine* e, hid_t grp, char* fileName,
  * @param with_hydro Are we reading gas particles ?
  * @param with_gravity Are we reading/creating #gpart arrays ?
  * @param with_stars Are we reading star particles ?
+ * @param cleanup_h Are we cleaning-up h-factors from the quantities we read?
+ * @param h The value of the reduced Hubble constant to use for correction.
  * @param mpi_rank The MPI rank of this node
  * @param mpi_size The number of MPI ranks
  * @param comm The MPI communicator
@@ -402,8 +423,8 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
                     struct spart** sparts, size_t* Ngas, size_t* Ngparts,
                     size_t* Nstars, int* periodic, int* flag_entropy,
                     int with_hydro, int with_gravity, int with_stars,
-                    int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
-                    int n_threads, int dry_run) {
+                    int cleanup_h, double h, 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) */
@@ -476,6 +497,13 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
     else if (hydro_dimension == 1)
       dim[2] = dim[1] = dim[0];
 
+    /* Convert the box size if we want to clean-up h-factors */
+    if (cleanup_h) {
+      dim[0] /= h;
+      dim[1] /= h;
+      dim[2] /= h;
+    }
+
     /* message("Found %lld particles in a %speriodic box of size [%f %f %f].",
      */
     /* 	    N_total, (periodic ? "": "non-"), dim[0], dim[1], dim[2]); */
@@ -641,7 +669,7 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
         if (!dry_run)
           for (int i = 0; i < num_fields; ++i)
             readArray(h_grp, list[i], Nparticles, N_total[ptype], offset[ptype],
-                      internal_units, ic_units);
+                      internal_units, ic_units, cleanup_h, h);
 
         /* Close particle group */
         H5Gclose(h_grp);

src/serial_io.h

@@ -22,6 +22,8 @@
 /* Config parameters. */
 #include "../config.h"
 
+#if defined(HAVE_HDF5) && defined(WITH_MPI) && !defined(HAVE_PARALLEL_HDF5)
+
 /* MPI headers. */
 #ifdef WITH_MPI
 #include <mpi.h>
@@ -32,15 +34,13 @@
 #include "part.h"
 #include "units.h"
 
-#if defined(HAVE_HDF5) && defined(WITH_MPI) && !defined(HAVE_PARALLEL_HDF5)
-
 void read_ic_serial(char* fileName, const struct unit_system* internal_units,
                     double dim[3], struct part** parts, struct gpart** gparts,
                     struct spart** sparts, size_t* Ngas, size_t* Ngparts,
                     size_t* Nstars, int* periodic, int* flag_entropy,
                     int with_hydro, int with_gravity, int with_stars,
-                    int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
-                    int nr_threads, int dry_run);
+                    int cleanup_h, double h, int mpi_rank, int mpi_size,
+                    MPI_Comm comm, MPI_Info info, int nr_threads, int dry_run);
 
 void write_output_serial(struct engine* e, const char* baseName,
                          const struct unit_system* internal_units,