diff --git a/examples/EAGLE_6/eagle_6.yml b/examples/EAGLE_6/eagle_6.yml
index d0055607d3cd7c1453e3b93ed8b21cd5ff2673d5..0f457ba5120d17d81273435f9752fb2e046523b2 100644
--- a/examples/EAGLE_6/eagle_6.yml
+++ b/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
diff --git a/examples/main.c b/examples/main.c
index ab1dfad6da036c781d2e07f6a0672905d1220574..1948110c1bb68a8cec641f7c23be90d49ef9da04 100644
--- a/examples/main.c
+++ b/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);
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 2dfc112435e797d9c024eff3dfbaa2eb67445e56..07c02b76bfd36cd9021e3420d3fce9a1383b4fb4 100644
--- a/examples/parameter_example.yml
+++ b/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).
diff --git a/src/cosmology.c b/src/cosmology.c
index b1bf71ad25451d2185d911320fb757193e990945..15a2d09f909f5de2675ec8de12ab400c92d18beb 100644
--- a/src/cosmology.c
+++ b/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.;
diff --git a/src/engine.c b/src/engine.c
index 64f557b5d6e3c76623a8247af873b3f359625b9f..9abb75936052c1f5fad6c4efe6170467a3111a49 100644
--- a/src/engine.c
+++ b/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) {
diff --git a/src/single_io.c b/src/single_io.c
index 0b9057dfb24aab52f70ae6f375d383e2a4994697..6950a1c7b75fc6f5db240a7a612de09601fccb4b 100644
--- a/src/single_io.c
+++ b/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);
diff --git a/src/single_io.h b/src/single_io.h
index efed4ef27ade61f7726777d71a2fdc586ea26030..26b849716e3e018d9a10c5c5c513ad26c7ccb274 100644
--- a/src/single_io.h
+++ b/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 */