Commit a7bfbf6f authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Add an option to cancel the sqrt(a) factors in the velocities of Gadget ICs.

parent a5d29e70
......@@ -644,6 +644,8 @@ int main(int argc, char *argv[]) {
params, "InitialConditions:cleanup_smoothing_lengths", 0);
const int cleanup_h = parser_get_opt_param_int(
params, "InitialConditions:cleanup_h_factors", 0);
const int cleanup_sqrt_a = parser_get_opt_param_int(
params, "InitialConditions:cleanup_velocity_factors", 0);
const int generate_gas_in_ics = parser_get_opt_param_int(
params, "InitialConditions:generate_gas_in_ics", 0);
if (generate_gas_in_ics && flag_entropy_ICs)
......@@ -653,6 +655,8 @@ int main(int argc, char *argv[]) {
if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
if (myrank == 0 && cleanup_h)
message("Cleaning up h-factors (h=%f)", cosmo.h);
if (myrank == 0 && cleanup_sqrt_a)
message("Cleaning up a-factors from velocity (a=%f)", cosmo.a);
fflush(stdout);
/* Get ready to read particles of all kinds */
......@@ -666,20 +670,23 @@ 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,
cleanup_h, cosmo.h, myrank, nr_nodes, MPI_COMM_WORLD,
MPI_INFO_NULL, nr_threads, dry_run);
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, &Ngas,
&Ngpart, &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
(with_external_gravity || with_self_gravity), with_stars,
cleanup_h, cosmo.h, myrank, nr_nodes, MPI_COMM_WORLD,
MPI_INFO_NULL, nr_threads, dry_run);
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, &Ngas,
&Ngpart, &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
(with_external_gravity || with_self_gravity), with_stars,
cleanup_h, cosmo.h, nr_threads, dry_run);
cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads,
dry_run);
#endif
#endif
if (myrank == 0) {
......
......@@ -94,6 +94,7 @@ InitialConditions:
file_name: SedovBlast/sedov.hdf5 # The file to read
generate_gas_in_ics: 0 # (Optional) Generate gas particles from the DM-only ICs (e.g. from panphasia).
cleanup_h_factors: 0 # (Optional) Clean up the h-factors used in the ICs (e.g. in Gadget files).
cleanup_velocity_factors: 0 # (Optional) Clean up the scale-factors used in the definition of the velocity variable 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).
......
......@@ -77,7 +77,7 @@ 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, int cleanup_h,
double h) {
int cleanup_sqrt_a, double h, double a) {
const size_t typeSize = io_sizeof_type(props.type);
const size_t copySize = typeSize * props.dimension;
......@@ -141,20 +141,35 @@ 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 (cleanup_h && h_factor_exp != 0.f) {
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;
const double h_factor = pow(h, h_factor_exp);
for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= h_factor;
} else {
float* temp_f = (float*)temp;
const float h_factor = pow(h, h_factor_exp);
for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= h_factor;
}
}
/* Clean-up a if necessary */
if (cleanup_sqrt_a && a != 1. && (strcmp(props.name, "Velocities") == 0)) {
if (io_is_double_precision(props.type)) {
double* temp_d = (double*)temp;
const double vel_factor = sqrt(a);
for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= vel_factor;
} else {
float* temp_f = (float*)temp;
const float vel_factor = sqrt(a);
for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= vel_factor;
}
}
/* Copy temporary buffer to particle data */
char* temp_c = (char*)temp;
for (size_t i = 0; i < N; ++i)
......@@ -183,7 +198,8 @@ void readArray_chunk(hid_t h_data, hid_t h_plist_id,
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, int cleanup_h, double h) {
const struct unit_system* ic_units, int cleanup_h,
int cleanup_sqrt_a, double h, double a) {
const size_t typeSize = io_sizeof_type(props.type);
const size_t copySize = typeSize * props.dimension;
......@@ -273,7 +289,7 @@ void readArray(hid_t grp, struct io_props props, size_t N, long long N_total,
/* Write the first chunk */
const size_t this_chunk = (N > max_chunk_size) ? max_chunk_size : N;
readArray_chunk(h_data, h_plist_id, props, this_chunk, offset,
internal_units, ic_units, cleanup_h, h);
internal_units, ic_units, cleanup_h, cleanup_sqrt_a, h, a);
/* Compute how many items are left */
if (N > max_chunk_size) {
......@@ -601,7 +617,10 @@ 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 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?
* @param h The value of the reduced Hubble constant to use for correction.
* @param a The current value of the scale-factor.
* @param mpi_rank The MPI rank of this node
* @param mpi_size The number of MPI ranks
* @param comm The MPI communicator
......@@ -615,9 +634,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 cleanup_h, double h, int mpi_rank, int mpi_size,
MPI_Comm comm, MPI_Info info, int n_threads,
int dry_run) {
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) */
......@@ -836,7 +855,8 @@ 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, cleanup_h, h);
offset[ptype], internal_units, ic_units, cleanup_h,
cleanup_sqrt_a, h, a);
/* Close particle group */
H5Gclose(h_grp);
......
......@@ -39,9 +39,9 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
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 cleanup_h, double h, int mpi_rank, int mpi_size,
MPI_Comm comm, MPI_Info info, int nr_threads,
int dry_run);
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,
......
......@@ -77,7 +77,8 @@
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, int cleanup_h, double h) {
const struct unit_system* ic_units, int cleanup_h,
int cleanup_sqrt_a, double h, double a) {
const size_t typeSize = io_sizeof_type(props.type);
const size_t copySize = typeSize * props.dimension;
......@@ -159,20 +160,35 @@ 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 (cleanup_h && 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;
const double h_factor = pow(h, h_factor_exp);
for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= h_factor;
} else {
float* temp_f = (float*)temp;
const float h_factor = pow(h, h_factor_exp);
for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= h_factor;
}
}
/* Clean-up a if necessary */
if (cleanup_sqrt_a && a != 1. && (strcmp(props.name, "Velocities") == 0)) {
if (io_is_double_precision(props.type)) {
double* temp_d = (double*)temp;
const double vel_factor = sqrt(a);
for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= vel_factor;
} else {
float* temp_f = (float*)temp;
const float vel_factor = sqrt(a);
for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= vel_factor;
}
}
/* Copy temporary buffer to particle data */
char* temp_c = temp;
for (size_t i = 0; i < N; ++i)
......@@ -393,7 +409,10 @@ 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 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?
* @param h The value of the reduced Hubble constant to use for correction.
* @param a The current value of the scale-factor.
* @param mpi_rank The MPI rank of this node
* @param mpi_size The number of MPI ranks
* @param comm The MPI communicator
......@@ -414,8 +433,9 @@ 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 cleanup_h, double h, int mpi_rank, int mpi_size,
MPI_Comm comm, MPI_Info info, int n_threads, int dry_run) {
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) */
......@@ -656,7 +676,8 @@ 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, cleanup_h, h);
internal_units, ic_units, cleanup_h, cleanup_sqrt_a, h,
a);
/* Close particle group */
H5Gclose(h_grp);
......
......@@ -39,8 +39,9 @@ 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 cleanup_h, double h, int mpi_rank, int mpi_size,
MPI_Comm comm, MPI_Info info, int nr_threads, int dry_run);
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_serial(struct engine* e, const char* baseName,
const struct unit_system* internal_units,
......
......@@ -73,7 +73,8 @@
*/
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, int cleanup_h, double h) {
const struct unit_system* ic_units, int cleanup_h, double h,
int cleanup_sqrt_a, double a) {
const size_t typeSize = io_sizeof_type(prop.type);
const size_t copySize = typeSize * prop.dimension;
......@@ -135,20 +136,35 @@ void readArray(hid_t h_grp, const struct io_props prop, size_t N,
/* Clean-up h if necessary */
const float h_factor_exp = units_h_factor(internal_units, prop.units);
if (cleanup_h && 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;
const double h_factor = pow(h, h_factor_exp);
for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= h_factor;
} else {
float* temp_f = (float*)temp;
const float h_factor = pow(h, h_factor_exp);
for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= h_factor;
}
}
/* Clean-up a if necessary */
if (cleanup_sqrt_a && a != 1. && (strcmp(prop.name, "Velocities") == 0)) {
if (io_is_double_precision(prop.type)) {
double* temp_d = (double*)temp;
const double vel_factor = sqrt(a);
for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= vel_factor;
} else {
float* temp_f = (float*)temp;
const float vel_factor = sqrt(a);
for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= vel_factor;
}
}
/* Copy temporary buffer to particle data */
char* temp_c = (char*)temp;
for (size_t i = 0; i < N; ++i)
......@@ -309,7 +325,10 @@ 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 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?
* @param h The value of the reduced Hubble constant to use for correction.
* @param a The current value of the scale-factor.
* @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.
*
......@@ -319,7 +338,6 @@ void writeArray(const struct engine* e, hid_t grp, char* fileName,
*
* @warning Can not read snapshot distributed over more than 1 file !!!
* @todo Read snapshots distributed in more than one file.
*
*/
void read_ic_single(const char* fileName,
const struct unit_system* internal_units, double dim[3],
......@@ -327,7 +345,8 @@ void read_ic_single(const char* fileName,
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 cleanup_h, double h, int n_threads, int dry_run) {
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) */
......@@ -533,7 +552,7 @@ void read_ic_single(const char* fileName,
if (!dry_run)
for (int i = 0; i < num_fields; ++i)
readArray(h_grp, list[i], Nparticles, internal_units, ic_units,
cleanup_h, h);
cleanup_h, cleanup_sqrt_a, h, a);
/* Close particle group */
H5Gclose(h_grp);
......
......@@ -35,7 +35,8 @@ void read_ic_single(const char* fileName,
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 cleanup_h, double h, int nr_threads, int dry_run);
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,
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment