Added facility to clean-up h-factors out of ICs when using the single_io mode.

examples/EAGLE_6/eagle_6.yml

@@ -50,4 +50,5 @@ SPH:
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  ./EAGLE_ICs_6.hdf5     # The file to read
+  cleanup_h_factors: 1               # Remove the h-factors inherited from Gadget
 

examples/main.c

@@ -622,7 +622,11 @@ int main(int argc, char *argv[]) {
         parser_get_opt_param_int(params, "InitialConditions:replicate", 1);
     clean_smoothing_length_values = parser_get_opt_param_int(
         params, "InitialConditions:cleanup_smoothing_lengths", 0);
+    const int cleanup_h = parser_get_opt_param_int(
+        params, "InitialConditions:cleanup_h_factors", 0);
     if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
+    if (myrank == 0 && cleanup_h)
+      message("Cleaning up h-factors (h=%f)", cosmo.h);
     fflush(stdout);
 
     /* Get ready to read particles of all kinds */
@@ -648,7 +652,7 @@ int main(int argc, char *argv[]) {
     read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
                    &Ngpart, &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
                    (with_external_gravity || with_self_gravity), with_stars,
-                   nr_threads, dry_run);
+                   cleanup_h, cosmo.h, nr_threads, dry_run);
 #endif
     if (myrank == 0) {
       clocks_gettime(&toc);

examples/parameter_example.yml

@@ -81,6 +81,7 @@ Statistics:
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  SedovBlast/sedov.hdf5 # The file to read
+  cleanup_h_factors:           0    # (Optional) Clean up the h-factors used in the ICs (e.g. in Gadget files).
   cleanup_smoothing_lengths:   0    # (Optional) Clean the values of the smoothing lengths that are read in to remove stupid values. Set to 1 to activate.
   smoothing_length_scaling:    1.   # (Optional) A scaling factor to apply to all smoothing lengths in the ICs.
   shift_x:    0.                    # (Optional) A shift to apply to all particles read from the ICs (in internal units).

src/cosmology.c

@@ -454,7 +454,7 @@ void cosmology_init_no_cosmo(struct cosmology *c) {
   c->Omega_b = 0.;
   c->w_0 = 0.;
   c->w_a = 0.;
-  c->h = 0.;
+  c->h = 1.;
   c->w = 0.;
 
   c->a_begin = 1.;

src/engine.c

@@ -3665,7 +3665,7 @@ int engine_estimate_nr_tasks(struct engine *e) {
  * @brief Rebuild the space and tasks.
  *
  * @param e The #engine.
- * @param clean_smoothing_length_values Are we cleaning up the values of 
+ * @param clean_smoothing_length_values Are we cleaning up the values of
  * the smoothing lengths before building the tasks ?
  */
 void engine_rebuild(struct engine *e, int clean_smoothing_length_values) {

src/single_io.c

@@ -64,13 +64,15 @@
  * @param N The number of particles.
  * @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 The value of the reduced Hubble constant.
  *
  * @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.
  */
 void readArray(hid_t h_grp, const struct io_props prop, size_t N,
                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(prop.type);
   const size_t copySize = typeSize * prop.dimension;
@@ -124,18 +126,35 @@ void readArray(hid_t h_grp, const struct io_props prop, size_t N,
   }
 
   /* Unit conversion if necessary */
-  const double factor =
+  const double unit_factor =
       units_conversion_factor(ic_units, internal_units, prop.units);
-  if (factor != 1. && exist != 0) {
+  if (unit_factor != 1. && exist != 0) {
 
     /* message("Converting ! factor=%e", factor); */
 
     if (io_is_double_precision(prop.type)) {
       double* temp_d = (double*)temp;
-      for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= factor;
+      for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= unit_factor;
     } else {
       float* temp_f = (float*)temp;
-      for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= factor;
+      for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= unit_factor;
+    }
+  }
+
+  /* Clean-up h if necessary */
+  const float h_factor_exp = units_h_factor(internal_units, prop.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", prop.name, h_factor_exp,
+     * h_factor); */
+
+    if (io_is_double_precision(prop.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;
     }
   }
 
@@ -301,6 +320,9 @@ 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.
+ * @prarm n_threads The number of threads to use for the temporary threadpool.
  * @param dry_run If 1, don't read the particle. Only allocates the arrays.
  *
  * Opens the HDF5 file fileName and reads the particles contained
@@ -316,7 +338,7 @@ void read_ic_single(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 n_threads, int dry_run) {
+                    int cleanup_h, double h, int n_threads, int dry_run) {
 
   hid_t h_file = 0, h_grp = 0;
   /* GADGET has only cubic boxes (in cosmological mode) */
@@ -384,6 +406,13 @@ void read_ic_single(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 %d particles in a %speriodic box of size [%f %f %f].",  */
   /* 	  *N, (periodic ? "": "non-"), dim[0], dim[1], dim[2]);  */
 
@@ -517,7 +546,8 @@ void read_ic_single(char* fileName, const struct unit_system* internal_units,
     /* Read everything */
     if (!dry_run)
       for (int i = 0; i < num_fields; ++i)
-        readArray(h_grp, list[i], Nparticles, internal_units, ic_units);
+        readArray(h_grp, list[i], Nparticles, internal_units, ic_units,
+                  cleanup_h, h);
 
     /* Close particle group */
     H5Gclose(h_grp);

src/single_io.h

@@ -22,24 +22,24 @@
 /* Config parameters. */
 #include "../config.h"
 
+#if defined(HAVE_HDF5) && !defined(WITH_MPI)
+
 /* Includes. */
 #include "engine.h"
 #include "part.h"
 #include "units.h"
 
-#if defined(HAVE_HDF5) && !defined(WITH_MPI)
-
 void read_ic_single(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* Ndm,
                     size_t* Nstars, int* periodic, int* flag_entropy,
                     int with_hydro, int with_gravity, int with_stars,
-                    int nr_threads, int dry_run);
+                    int cleanup_h, double h, int nr_threads, int dry_run);
 
 void write_output_single(struct engine* e, const char* baseName,
                          const struct unit_system* internal_units,
                          const struct unit_system* snapshot_units);
 
-#endif
+#endif /* HAVE_HDF5 && !WITH_MPI */
 
 #endif /* SWIFT_SINGLE_IO_H */