diff --git a/examples/CosmoVolume/cosmoVolume.yml b/examples/CosmoVolume/cosmoVolume.yml
index 83d70f35c0d7ed5f0c81eb138b37eb0881125e85..548b85078952954f2bf97280be6feb25eb6ef444 100644
--- a/examples/CosmoVolume/cosmoVolume.yml
+++ b/examples/CosmoVolume/cosmoVolume.yml
@@ -23,11 +23,6 @@ Snapshots:
basename: cosmo # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.05 # Time difference between consecutive outputs (in internal units)
- UnitMass_in_cgs: 1 # Grams
- UnitLength_in_cgs: 1 # Centimeters
- UnitVelocity_in_cgs: 1 # Centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
# Parameters governing the conserved quantities statistics
Statistics:
diff --git a/examples/EAGLE_12/eagle_12.yml b/examples/EAGLE_12/eagle_12.yml
index 778d04fdd02004fd10ac468321568069e39dd217..bac3e4684d4436403bb2497afd34a865ea4bab87 100644
--- a/examples/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_12/eagle_12.yml
@@ -1,10 +1,10 @@
# Define the system of units to use internally.
InternalUnitSystem:
- UnitMass_in_cgs: 1 # Grams
- UnitLength_in_cgs: 1 # Centimeters
- UnitVelocity_in_cgs: 1 # Centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
+ UnitMass_in_cgs: 1.989e43 # 10^10 M_sun in grams
+ UnitLength_in_cgs: 3.085678e24 # Mpc in centimeters
+ UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
# Parameters for the task scheduling
Scheduler:
@@ -14,20 +14,15 @@ Scheduler:
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 1. # The end time of the simulation (in internal units).
- dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
- dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
+ time_end: 1e-2 # The end time of the simulation (in internal units).
+ dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: eagle # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
- delta_time: 0.05 # Time difference between consecutive outputs (in internal units)
- UnitMass_in_cgs: 1 # Grams
- UnitLength_in_cgs: 1 # Centimeters
- UnitVelocity_in_cgs: 1 # Centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
+ delta_time: 1e-3 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
diff --git a/examples/EAGLE_25/eagle_25.yml b/examples/EAGLE_25/eagle_25.yml
index 731d87ea360ced596c4190880c62fdfb1cd605f8..cd5a1211b68c4ec2bebcd565713cc52a7b8de6df 100644
--- a/examples/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_25/eagle_25.yml
@@ -1,10 +1,10 @@
# Define the system of units to use internally.
InternalUnitSystem:
- UnitMass_in_cgs: 1 # Grams
- UnitLength_in_cgs: 1 # Centimeters
- UnitVelocity_in_cgs: 1 # Centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
+ UnitMass_in_cgs: 1.989e43 # 10^10 M_sun in grams
+ UnitLength_in_cgs: 3.085678e24 # Mpc in centimeters
+ UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
# Parameters for the task scheduling
Scheduler:
@@ -14,20 +14,15 @@ Scheduler:
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 1. # The end time of the simulation (in internal units).
- dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
- dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
+ time_end: 1e-2 # The end time of the simulation (in internal units).
+ dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: eagle # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
- delta_time: 0.05 # Time difference between consecutive outputs (in internal units)
- UnitMass_in_cgs: 1 # Grams
- UnitLength_in_cgs: 1 # Centimeters
- UnitVelocity_in_cgs: 1 # Centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
+ delta_time: 1e-3 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
diff --git a/examples/EAGLE_50/eagle_50.yml b/examples/EAGLE_50/eagle_50.yml
index 1b1f4e4dfa0559f4fcf18d99c0fc649a3b8f5307..38a9c165796359dfdabef423154299ee04ae8c76 100644
--- a/examples/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_50/eagle_50.yml
@@ -1,10 +1,10 @@
# Define the system of units to use internally.
InternalUnitSystem:
- UnitMass_in_cgs: 1 # Grams
- UnitLength_in_cgs: 1 # Centimeters
- UnitVelocity_in_cgs: 1 # Centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
+ UnitMass_in_cgs: 1.989e43 # 10^10 M_sun in grams
+ UnitLength_in_cgs: 3.085678e24 # Mpc in centimeters
+ UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
# Parameters for the task scheduling
Scheduler:
@@ -14,20 +14,15 @@ Scheduler:
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 1. # The end time of the simulation (in internal units).
- dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
- dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
+ time_end: 1e-2 # The end time of the simulation (in internal units).
+ dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: eagle # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
- delta_time: 0.05 # Time difference between consecutive outputs (in internal units)
- UnitMass_in_cgs: 1 # Grams
- UnitLength_in_cgs: 1 # Centimeters
- UnitVelocity_in_cgs: 1 # Centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
+ delta_time: 1e-3 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
diff --git a/examples/ExternalPointMass/externalPointMass.yml b/examples/ExternalPointMass/externalPointMass.yml
index 52330163caa13609fd7674a5cdf2921743fbe227..c6ab8af09ccc4f659c241f2ec91ff86f16043ac6 100644
--- a/examples/ExternalPointMass/externalPointMass.yml
+++ b/examples/ExternalPointMass/externalPointMass.yml
@@ -18,11 +18,6 @@ Snapshots:
basename: pointMass # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.02 # Time difference between consecutive outputs (in internal units)
- UnitMass_in_cgs: 1.9885e33 # Grams
- UnitLength_in_cgs: 3.0856776e21 # Centimeters
- UnitVelocity_in_cgs: 1e5 # Centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
# Parameters governing the conserved quantities statistics
Statistics:
diff --git a/examples/MultiTypes/multiTypes.yml b/examples/MultiTypes/multiTypes.yml
index a1ff19f98e9fccab82d1749a45adf026b6675907..28d02fefa8168e35af696975a7c73a1bf767155e 100644
--- a/examples/MultiTypes/multiTypes.yml
+++ b/examples/MultiTypes/multiTypes.yml
@@ -18,11 +18,6 @@ Snapshots:
basename: multiTypes # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.01 # Time difference between consecutive outputs (in internal units)
- UnitMass_in_cgs: 1 # Grams
- UnitLength_in_cgs: 1 # Centimeters
- UnitVelocity_in_cgs: 1 # Centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
# Parameters governing the conserved quantities statistics
Statistics:
diff --git a/examples/SedovBlast/sedov.yml b/examples/SedovBlast/sedov.yml
index 9fbabb4969b6accdb7323d8270b735951ac0693a..eb95fd85d5599145b2ff037256dcde72fc306209 100644
--- a/examples/SedovBlast/sedov.yml
+++ b/examples/SedovBlast/sedov.yml
@@ -18,11 +18,6 @@ Snapshots:
basename: sedov # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.1 # Time difference between consecutive outputs (in internal units)
- UnitMass_in_cgs: 1 # Grams
- UnitLength_in_cgs: 1 # Centimeters
- UnitVelocity_in_cgs: 1 # Centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
# Parameters governing the conserved quantities statistics
Statistics:
diff --git a/examples/SodShock/sodShock.yml b/examples/SodShock/sodShock.yml
index 003b7286777230ceb3b84ee01c6a1f335aeb9476..19bacb57f3bfb173063dbac5d752929763dede29 100644
--- a/examples/SodShock/sodShock.yml
+++ b/examples/SodShock/sodShock.yml
@@ -18,11 +18,6 @@ Snapshots:
basename: sod # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.1 # Time difference between consecutive outputs (in internal units)
- UnitMass_in_cgs: 1 # Grams
- UnitLength_in_cgs: 1 # Centimeters
- UnitVelocity_in_cgs: 1 # Centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
# Parameters governing the conserved quantities statistics
Statistics:
diff --git a/examples/UniformBox/uniformBox.yml b/examples/UniformBox/uniformBox.yml
index 069a17fb549775acf918da88221f7b3ed7e80565..7c9c74e1342bffb939131a265188cae269cc773f 100644
--- a/examples/UniformBox/uniformBox.yml
+++ b/examples/UniformBox/uniformBox.yml
@@ -18,11 +18,6 @@ Snapshots:
basename: uniformBox # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.01 # Time difference between consecutive outputs (in internal units)
- UnitMass_in_cgs: 1 # Grams
- UnitLength_in_cgs: 1 # Centimeters
- UnitVelocity_in_cgs: 1 # Centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
# Parameters governing the conserved quantities statistics
Statistics:
diff --git a/examples/main.c b/examples/main.c
index a592c974feda86fd7c6537017a7b31480f4a1bc2..99f673307ab43b5e55991591177da30f4c212f63 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -358,8 +358,8 @@ int main(int argc, char *argv[]) {
&periodic, &flag_entropy_ICs, myrank, nr_nodes,
MPI_COMM_WORLD, MPI_INFO_NULL, dry_run);
#else
- read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic,
- &flag_entropy_ICs, myrank, nr_nodes, MPI_COMM_WORLD,
+ read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &Ngas, &Ngpart,
+ &periodic, &flag_entropy_ICs, myrank, nr_nodes, MPI_COMM_WORLD,
MPI_INFO_NULL, dry_run);
#endif
#else
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index b76b127f02f8253d5738ebdaae2a685dc74f96e0..b2000e978f48ced279e365f24bb717202584f3b1 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -25,11 +25,11 @@ Snapshots:
basename: output # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.01 # Time difference between consecutive outputs (in internal units)
- UnitMass_in_cgs: 1 # Unit system for the outputs (Grams)
- UnitLength_in_cgs: 1 # Unit system for the outputs (Centimeters)
- UnitVelocity_in_cgs: 1 # Unit system for the outputs (Centimeters per second)
- UnitCurrent_in_cgs: 1 # Unit system for the outputs (Amperes)
- UnitTemp_in_cgs: 1 # Unit system for the outputs (Kelvin)
+ UnitMass_in_cgs: 1 # (Optional) Unit system for the outputs (Grams)
+ UnitLength_in_cgs: 1 # (Optional) Unit system for the outputs (Centimeters)
+ UnitVelocity_in_cgs: 1 # (Optional) Unit system for the outputs (Centimeters per second)
+ UnitCurrent_in_cgs: 1 # (Optional) Unit system for the outputs (Amperes)
+ UnitTemp_in_cgs: 1 # (Optional) Unit system for the outputs (Kelvin)
# Parameters governing the conserved quantities statistics
Statistics:
diff --git a/src/adiabatic_index.h b/src/adiabatic_index.h
index 0df66d65a68bfaea63684e85e10c82940ded06b4..ba7a1a64d5c97f1f6d276e5969782351762be4b1 100644
--- a/src/adiabatic_index.h
+++ b/src/adiabatic_index.h
@@ -60,15 +60,12 @@ __attribute__((always_inline)) INLINE static float pow_gamma(float x) {
#if defined(HYDRO_GAMMA_5_3)
- const float x2 = x * x;
- const float x5 = x2 * x2 * x;
- return cbrtf(x5);
+ const float cbrt = cbrtf(x); /* x^(1/3) */
+ return cbrt * cbrt * x; /* x^(5/3) */
#elif defined(HYDRO_GAMMA_4_3)
- const float x2 = x * x;
- const float x4 = x2 * x2;
- return cbrtf(x4);
+ return cbrtf(x) * x; /* x^(4/3) */
#elif defined(HYDRO_GAMMA_2_1)
@@ -93,11 +90,12 @@ __attribute__((always_inline)) INLINE static float pow_gamma_minus_one(
#if defined(HYDRO_GAMMA_5_3)
- return cbrtf(x * x);
+ const float cbrt = cbrtf(x); /* x^(1/3) */
+ return cbrt * cbrt; /* x^(2/3) */
#elif defined(HYDRO_GAMMA_4_3)
- return cbrtf(x);
+ return cbrtf(x); /* x^(1/3) */
#elif defined(HYDRO_GAMMA_2_1)
@@ -122,11 +120,12 @@ __attribute__((always_inline)) INLINE static float pow_minus_gamma_minus_one(
#if defined(HYDRO_GAMMA_5_3)
- return 1.f / cbrtf(x * x);
+ const float inv_cbrt = 1.f / cbrtf(x); /* x^(-1/3) */
+ return inv_cbrt * inv_cbrt; /* x^(-2/3) */
#elif defined(HYDRO_GAMMA_4_3)
- return 1.f / cbrtf(x);
+ return 1.f / cbrtf(x); /* x^(-1/3) */
#elif defined(HYDRO_GAMMA_2_1)
diff --git a/src/common_io.c b/src/common_io.c
index 7108ff8d9209011cba5e131220d8e33fb1de00cf..d3590ae402487eb2f605a1cf3f5722f86704ad88 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -120,7 +120,7 @@ size_t sizeOfType(enum DATA_TYPE type) {
* Returns an error if the type is not FLOAT or DOUBLE
*/
int isDoublePrecision(enum DATA_TYPE type) {
-
+
switch (type) {
case FLOAT:
return 0;
diff --git a/src/drift.h b/src/drift.h
index 05b09bb7910e8cddaf9fb24bb248f120f9db9eea..932e956d8c9f9c8b1cb5afdbff33f99bc79c180d 100644
--- a/src/drift.h
+++ b/src/drift.h
@@ -37,7 +37,8 @@
* @param ti_current Integer end of time-step
*/
__attribute__((always_inline)) INLINE static void drift_gpart(
- struct gpart* gp, float dt, double timeBase, int ti_old, int ti_current) {
+ struct gpart *restrict gp, float dt, double timeBase, int ti_old,
+ int ti_current) {
/* Drift... */
gp->x[0] += gp->v_full[0] * dt;
gp->x[1] += gp->v_full[1] * dt;
@@ -60,8 +61,8 @@ __attribute__((always_inline)) INLINE static void drift_gpart(
* @param ti_current Integer end of time-step
*/
__attribute__((always_inline)) INLINE static void drift_part(
- struct part* p, struct xpart* xp, float dt, double timeBase, int ti_old,
- int ti_current) {
+ struct part *restrict p, struct xpart *restrict xp, float dt,
+ double timeBase, int ti_old, int ti_current) {
/* Useful quantity */
const float h_inv = 1.0f / p->h;
diff --git a/src/engine.c b/src/engine.c
index 96d72669ca6eb4282e3baa4276fab2e7987bd97a..1ae4dff2517b22aa4760a690c082b64743f54243 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2672,7 +2672,7 @@ void engine_dump_snapshot(struct engine *e) {
struct clocks_time time1, time2;
clocks_gettime(&time1);
- if (e->verbose) message("writing snapshot at t=%f.", e->time);
+ if (e->verbose) message("writing snapshot at t=%e.", e->time);
/* Dump... */
#if defined(WITH_MPI)
@@ -2808,7 +2808,7 @@ void engine_init(struct engine *e, struct space *s,
e->ti_nextSnapshot = 0;
parser_get_param_string(params, "Snapshots:basename", e->snapshotBaseName);
e->snapshotUnits = malloc(sizeof(struct UnitSystem));
- units_init(e->snapshotUnits, params, "Snapshots");
+ units_init_default(e->snapshotUnits, params, "Snapshots", internal_units);
e->dt_min = parser_get_param_double(params, "TimeIntegration:dt_min");
e->dt_max = parser_get_param_double(params, "TimeIntegration:dt_max");
e->file_stats = NULL;
@@ -3228,6 +3228,6 @@ void engine_compute_next_snapshot_time(struct engine *e) {
const float next_snapshot_time =
e->ti_nextSnapshot * e->timeBase + e->timeBegin;
if (e->verbose)
- message("Next output time set to t=%f.", next_snapshot_time);
+ message("Next output time set to t=%e.", next_snapshot_time);
}
}
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index a09920235eedcc6c997c1407d7f6c63d8f7a630f..ee693d60999779d5a3cd888e4737208868879ba8 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -28,8 +28,8 @@
*
*/
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
- const struct part* p, const struct xpart* xp,
- const struct hydro_props* hydro_properties) {
+ const struct part *restrict p, const struct xpart *restrict xp,
+ const struct hydro_props *restrict hydro_properties) {
const float CFL_condition = hydro_properties->CFL_condition;
@@ -55,7 +55,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
* @param xp The extended particle data to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
- struct part* p, struct xpart* xp) {}
+ struct part *restrict p, struct xpart *restrict xp) {}
/**
* @brief Prepares a particle for the density calculation.
@@ -66,15 +66,15 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
- struct part* p) {
+ struct part *restrict p) {
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
p->rho_dh = 0.f;
p->density.div_v = 0.f;
- p->density.curl_v[0] = 0.f;
- p->density.curl_v[1] = 0.f;
- p->density.curl_v[2] = 0.f;
+ p->density.rot_v[0] = 0.f;
+ p->density.rot_v[1] = 0.f;
+ p->density.rot_v[2] = 0.f;
}
/**
@@ -87,7 +87,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
* @param time The current time
*/
__attribute__((always_inline)) INLINE static void hydro_end_density(
- struct part* p, float time) {
+ struct part *restrict p, float time) {
/* Some smoothing length multiples. */
const float h = p->h;
@@ -110,6 +110,14 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* Compute the derivative term */
p->rho_dh = 1.f / (1.f + 0.33333333f * p->h * p->rho_dh * irho);
+
+ /* Finish calculation of the velocity curl components */
+ p->density.rot_v[0] *= ih4 * irho;
+ p->density.rot_v[1] *= ih4 * irho;
+ p->density.rot_v[2] *= ih4 * irho;
+
+ /* Finish calculation of the velocity divergence */
+ p->density.div_v *= ih4 * irho;
}
/**
@@ -122,7 +130,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
* @param time The current time
*/
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
- struct part* p, struct xpart* xp, int ti_current, double timeBase) {
+ struct part *restrict p, struct xpart *restrict xp, int ti_current,
+ double timeBase) {
/* Some smoothing length multiples. */
const float h = p->h;
@@ -132,10 +141,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Pre-compute some stuff for the balsara switch. */
const float normDiv_v = fabs(p->density.div_v / p->rho * ih4);
- const float normCurl_v = sqrtf(p->density.curl_v[0] * p->density.curl_v[0] +
- p->density.curl_v[1] * p->density.curl_v[1] +
- p->density.curl_v[2] * p->density.curl_v[2]) /
- p->rho * ih4;
+ const float normRot_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
+ p->density.rot_v[1] * p->density.rot_v[1] +
+ p->density.rot_v[2] * p->density.rot_v[2]) /
+ p->rho * ih4;
/* Compute this particle's sound speed. */
const float u = p->u;
@@ -146,7 +155,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
p->force.POrho2 = u * hydro_gamma_minus_one / (p->rho * xp->omega);
/* Balsara switch */
- p->force.balsara = normDiv_v / (normDiv_v + normCurl_v + 0.0001f * fc * ih);
+ p->force.balsara = normDiv_v / (normDiv_v + normRot_v + 0.0001f * fc * ih);
/* Viscosity parameter decay time */
const float tau = h / (2.f * const_viscosity_length * p->force.c);
@@ -171,7 +180,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
- struct part* p) {
+ struct part *restrict p) {
/* Reset the acceleration. */
p->a_hydro[0] = 0.0f;
@@ -194,7 +203,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
* @param timeBase The minimal time-step size
*/
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
- struct part* p, struct xpart* xp, int t0, int t1, double timeBase) {
+ struct part *restrict p, struct xpart *restrict xp, int t0, int t1,
+ double timeBase) {
float u, w;
const float dt = (t1 - t0) * timeBase;
@@ -218,7 +228,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_end_force(
- struct part* p) {}
+ struct part *restrict p) {}
/**
* @brief Kick the additional variables
@@ -229,17 +239,18 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param half_dt The half time-step for this kick
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
- struct part* p, struct xpart* xp, float dt, float half_dt) {}
+ struct part *restrict p, struct xpart *restrict xp, float dt,
+ float half_dt) {}
/**
- * @brief Converts hydro quantity of a particle
+ * @brief Converts hydro quantity of a particle at the start of a run
*
* Requires the density to be known
*
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
- struct part* p) {}
+ struct part *restrict p) {}
/**
* @brief Returns the internal energy of a particle
@@ -248,7 +259,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
- const struct part* p, float dt) {
+ const struct part *restrict p, float dt) {
return p->u;
}
diff --git a/src/hydro/Default/hydro_iact.h b/src/hydro/Default/hydro_iact.h
index e67d6422a8a3019beee7942f3381fabe5951640c..e4062f6a64821b06658d77c57254c91d97eb5f5e 100644
--- a/src/hydro/Default/hydro_iact.h
+++ b/src/hydro/Default/hydro_iact.h
@@ -79,8 +79,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
pi->density.wcount += wi;
pi->density.wcount_dh -= xi * wi_dx;
- pi->density.div_v += mj * dvdr * wi_dx;
- for (k = 0; k < 3; k++) pi->density.curl_v[k] += mj * curlvr[k] * wi_dx;
+ pi->density.div_v -= mj * dvdr * wi_dx;
+ for (k = 0; k < 3; k++) pi->density.rot_v[k] += mj * curlvr[k] * wi_dx;
/* Compute density of pj. */
h_inv = 1.0 / hj;
@@ -92,8 +92,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
pj->density.wcount += wj;
pj->density.wcount_dh -= xj * wj_dx;
- pj->density.div_v += mi * dvdr * wj_dx;
- for (k = 0; k < 3; k++) pj->density.curl_v[k] += mi * curlvr[k] * wj_dx;
+ pj->density.div_v -= mi * dvdr * wj_dx;
+ for (k = 0; k < 3; k++) pj->density.rot_v[k] += mi * curlvr[k] * wj_dx;
}
/**
@@ -111,7 +111,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
vector dx[3], dv[3];
vector vi[3], vj[3];
vector dvdr, div_vi, div_vj;
- vector curlvr[3], curl_vi[3], curl_vj[3];
+ vector curlvr[3], rot_vi[3], rot_vj[3];
int k, j;
#if VEC_SIZE == 8
@@ -178,28 +178,28 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
wcounti.v = wi.v;
wcounti_dh.v = xi.v * wi_dx.v;
div_vi.v = mj.v * dvdr.v * wi_dx.v;
- for (k = 0; k < 3; k++) curl_vi[k].v = mj.v * curlvr[k].v * wi_dx.v;
+ for (k = 0; k < 3; k++) rot_vi[k].v = mj.v * curlvr[k].v * wi_dx.v;
rhoj.v = mi.v * wj.v;
rhoj_dh.v = mi.v * (vec_set1(3.0f) * wj.v + xj.v * wj_dx.v);
wcountj.v = wj.v;
wcountj_dh.v = xj.v * wj_dx.v;
div_vj.v = mi.v * dvdr.v * wj_dx.v;
- for (k = 0; k < 3; k++) curl_vj[k].v = mi.v * curlvr[k].v * wj_dx.v;
+ for (k = 0; k < 3; k++) rot_vj[k].v = mi.v * curlvr[k].v * wj_dx.v;
for (k = 0; k < VEC_SIZE; k++) {
pi[k]->rho += rhoi.f[k];
pi[k]->rho_dh -= rhoi_dh.f[k];
pi[k]->density.wcount += wcounti.f[k];
pi[k]->density.wcount_dh -= wcounti_dh.f[k];
- pi[k]->density.div_v += div_vi.f[k];
- for (j = 0; j < 3; j++) pi[k]->density.curl_v[j] += curl_vi[j].f[k];
+ pi[k]->density.div_v -= div_vi.f[k];
+ for (j = 0; j < 3; j++) pi[k]->density.rot_v[j] += rot_vi[j].f[k];
pj[k]->rho += rhoj.f[k];
pj[k]->rho_dh -= rhoj_dh.f[k];
pj[k]->density.wcount += wcountj.f[k];
pj[k]->density.wcount_dh -= wcountj_dh.f[k];
- pj[k]->density.div_v += div_vj.f[k];
- for (j = 0; j < 3; j++) pj[k]->density.curl_v[j] += curl_vj[j].f[k];
+ pj[k]->density.div_v -= div_vj.f[k];
+ for (j = 0; j < 3; j++) pj[k]->density.rot_v[j] += rot_vj[j].f[k];
}
#else
@@ -254,8 +254,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
pi->density.wcount += wi;
pi->density.wcount_dh -= xi * wi_dx;
- pi->density.div_v += mj * dvdr * wi_dx;
- for (k = 0; k < 3; k++) pi->density.curl_v[k] += mj * curlvr[k] * wi_dx;
+ pi->density.div_v -= mj * dvdr * wi_dx;
+ for (k = 0; k < 3; k++) pi->density.rot_v[k] += mj * curlvr[k] * wi_dx;
}
/**
@@ -273,7 +273,7 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
vector dx[3], dv[3];
vector vi[3], vj[3];
vector dvdr;
- vector curlvr[3], curl_vi[3];
+ vector curlvr[3], rot_vi[3];
int k, j;
#if VEC_SIZE == 8
@@ -331,15 +331,15 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
wcounti.v = wi.v;
wcounti_dh.v = xi.v * wi_dx.v;
div_vi.v = mj.v * dvdr.v * wi_dx.v;
- for (k = 0; k < 3; k++) curl_vi[k].v = mj.v * curlvr[k].v * wi_dx.v;
+ for (k = 0; k < 3; k++) rot_vi[k].v = mj.v * curlvr[k].v * wi_dx.v;
for (k = 0; k < VEC_SIZE; k++) {
pi[k]->rho += rhoi.f[k];
pi[k]->rho_dh -= rhoi_dh.f[k];
pi[k]->density.wcount += wcounti.f[k];
pi[k]->density.wcount_dh -= wcounti_dh.f[k];
- pi[k]->density.div_v += div_vi.f[k];
- for (j = 0; j < 3; j++) pi[k]->density.curl_v[j] += curl_vi[j].f[k];
+ pi[k]->density.div_v -= div_vi.f[k];
+ for (j = 0; j < 3; j++) pi[k]->density.rot_v[j] += rot_vi[j].f[k];
}
#else
diff --git a/src/hydro/Default/hydro_part.h b/src/hydro/Default/hydro_part.h
index a4096d5c30d525307d5327559ef6b007c6931486..946145587262a1502f9e2ee847f10500c1a6b986 100644
--- a/src/hydro/Default/hydro_part.h
+++ b/src/hydro/Default/hydro_part.h
@@ -80,7 +80,7 @@ struct part {
float wcount_dh;
/* Particle velocity curl. */
- float curl_v[3];
+ float rot_v[3];
/* Particle number density. */
float wcount;
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 7789d4a80db6b37003fcbff92e5f8b42b1ca5ff3..a523fbb99f469004562fb7194d0c3825e9cd1857 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -27,8 +27,8 @@
*
*/
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
- const struct part* p, const struct xpart* xp,
- const struct hydro_props* hydro_properties) {
+ const struct part *restrict p, const struct xpart *restrict xp,
+ const struct hydro_props *restrict hydro_properties) {
const float CFL_condition = hydro_properties->CFL_condition;
@@ -49,7 +49,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
* @param xp The extended particle data to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
- struct part* p, struct xpart* xp) {}
+ struct part *restrict p, struct xpart *restrict xp) {}
/**
* @brief Prepares a particle for the density calculation.
@@ -60,7 +60,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
- struct part* p) {
+ struct part *restrict p) {
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
@@ -81,7 +81,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
* @param ti_current The current time (on the integer timeline)
*/
__attribute__((always_inline)) INLINE static void hydro_end_density(
- struct part* p, int ti_current) {
+ struct part *restrict p, int ti_current) {
/* Some smoothing length multiples. */
const float h = p->h;
@@ -125,7 +125,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
* @param timeBase The minimal time-step size
*/
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
- struct part* p, struct xpart* xp, int ti_current, double timeBase) {
+ struct part *restrict p, struct xpart *restrict xp, int ti_current,
+ double timeBase) {
/* Compute the norm of the curl */
p->force.curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
@@ -149,7 +150,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
- struct part* p) {
+ struct part *restrict p) {
/* Reset the acceleration. */
p->a_hydro[0] = 0.0f;
@@ -175,7 +176,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
* @param timeBase The minimal time-step size
*/
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
- struct part* p, struct xpart* xp, int t0, int t1, double timeBase) {
+ struct part *restrict p, const struct xpart *restrict xp, int t0, int t1,
+ double timeBase) {
/* Drift the pressure */
const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
@@ -194,7 +196,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_end_force(
- struct part* p) {
+ struct part *restrict p) {
p->entropy_dt *= hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho);
}
@@ -208,7 +210,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param half_dt The half time-step for this kick
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
- struct part* p, struct xpart* xp, float dt, float half_dt) {
+ struct part *restrict p, struct xpart *restrict xp, float dt,
+ float half_dt) {
/* Do not decrease the entropy (temperature) by more than a factor of 2*/
const float entropy_change = p->entropy_dt * dt;
@@ -223,14 +226,14 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
}
/**
- * @brief Converts hydro quantity of a particle
+ * @brief Converts hydro quantity of a particle at the start of a run
*
* Requires the density to be known
*
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
- struct part* p) {
+ struct part *restrict p) {
p->entropy =
hydro_gamma_minus_one * p->entropy * pow_minus_gamma_minus_one(p->rho);
@@ -243,7 +246,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
- const struct part* p, float dt) {
+ const struct part *restrict p, float dt) {
const float entropy = p->entropy + p->entropy_dt * dt;
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 848de3b5de8b916c35a3c6da7e414fd9a35d966b..87daae2bc43f671aa9549058bb8168c0af91ad5c 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -32,8 +32,8 @@
*
*/
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
- const struct part* p, const struct xpart* xp,
- const struct hydro_props* hydro_properties) {
+ const struct part *restrict p, const struct xpart *restrict xp,
+ const struct hydro_props *restrict hydro_properties) {
const float CFL_condition = hydro_properties->CFL_condition;
@@ -55,7 +55,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
* @param xp The extended particle data to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
- struct part* p, struct xpart* xp) {
+ struct part *restrict p, struct xpart *restrict xp) {
xp->u_full = p->u;
}
@@ -70,7 +70,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
- struct part* p) {
+ struct part *restrict p) {
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
@@ -90,7 +90,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
* @param time The current time
*/
__attribute__((always_inline)) INLINE static void hydro_end_density(
- struct part* p, float time) {
+ struct part *restrict p, float time) {
/* Some smoothing length multiples. */
const float h = p->h;
@@ -131,7 +131,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
* @param timeBase The minimal time-step size
*/
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
- struct part* p, struct xpart* xp, int ti_current, double timeBase) {
+ struct part *restrict p, struct xpart *restrict xp, int ti_current,
+ double timeBase) {
p->force.pressure = p->rho * p->u * hydro_gamma_minus_one;
}
@@ -145,7 +146,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
- struct part* p) {
+ struct part *restrict p) {
/* Reset the acceleration. */
p->a_hydro[0] = 0.0f;
@@ -171,7 +172,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
* @param timeBase The minimal time-step size
*/
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
- struct part* p, struct xpart* xp, int t0, int t1, double timeBase) {
+ struct part *restrict p, const struct xpart *restrict xp, int t0, int t1,
+ double timeBase) {
p->u = xp->u_full;
@@ -189,7 +191,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_end_force(
- struct part* p) {}
+ struct part *restrict p) {}
/**
* @brief Kick the additional variables
@@ -203,7 +205,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param half_dt The half time-step for this kick
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
- struct part* p, struct xpart* xp, float dt, float half_dt) {
+ struct part *restrict p, struct xpart *restrict xp, float dt,
+ float half_dt) {
/* Kick in momentum space */
xp->u_full += p->u_dt * dt;
@@ -223,7 +226,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
- struct part* p) {}
+ struct part *restrict p) {}
/**
* @brief Returns the internal energy of a particle
@@ -236,7 +239,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
- const struct part* p, float dt) {
+ const struct part *restrict p, float dt) {
return p->u;
}
diff --git a/src/kick.h b/src/kick.h
index df3bd4cc6ce9819d0b65640db51d2ed20a7d59fe..2cf50ca00ebe738a194d54fbd0142ed7e525b947 100644
--- a/src/kick.h
+++ b/src/kick.h
@@ -33,9 +33,8 @@
* @param new_dti The (integer) time-step for this kick.
* @param timeBase The minimal allowed time-step size.
*/
-__attribute__((always_inline)) INLINE static void kick_gpart(struct gpart* gp,
- int new_dti,
- double timeBase) {
+__attribute__((always_inline)) INLINE static void kick_gpart(
+ struct gpart *restrict gp, int new_dti, double timeBase) {
/* Compute the time step for this kick */
const int ti_start = (gp->ti_begin + gp->ti_end) / 2;
@@ -64,10 +63,9 @@ __attribute__((always_inline)) INLINE static void kick_gpart(struct gpart* gp,
* @param new_dti The (integer) time-step for this kick.
* @param timeBase The minimal allowed time-step size.
*/
-__attribute__((always_inline)) INLINE static void kick_part(struct part* p,
- struct xpart* xp,
- int new_dti,
- double timeBase) {
+__attribute__((always_inline)) INLINE static void kick_part(
+ struct part *restrict p, struct xpart *restrict xp, int new_dti,
+ double timeBase) {
/* Compute the time step for this kick */
const int ti_start = (p->ti_begin + p->ti_end) / 2;
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 47f2e893f9046dc9f517c1d362b57217a520a318..5e358c34f70065255e1eab61295b21487f07f77d 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -39,9 +39,9 @@
#include "common_io.h"
#include "engine.h"
#include "error.h"
-#include "hydro_properties.h"
#include "gravity_io.h"
#include "hydro_io.h"
+#include "hydro_properties.h"
#include "io_properties.h"
#include "kernel_hydro.h"
#include "part.h"
@@ -706,6 +706,39 @@ void write_output_parallel(struct engine* e, const char* baseName,
/* Print the system of Units used internally */
writeUnitSystem(h_file, internal_units, "InternalCodeUnits");
+ /* Tell the user if a conversion will be needed */
+ if (e->verbose && mpi_rank == 0) {
+ if (units_are_equal(snapshot_units, internal_units)) {
+
+ message("Snapshot and internal units match. No conversion needed.");
+
+ } else {
+
+ message("Conversion needed from:");
+ message("(Snapshot) Unit system: U_M = %e g.",
+ snapshot_units->UnitMass_in_cgs);
+ message("(Snapshot) Unit system: U_L = %e cm.",
+ snapshot_units->UnitLength_in_cgs);
+ message("(Snapshot) Unit system: U_t = %e s.",
+ snapshot_units->UnitTime_in_cgs);
+ message("(Snapshot) Unit system: U_I = %e A.",
+ snapshot_units->UnitCurrent_in_cgs);
+ message("(Snapshot) Unit system: U_T = %e K.",
+ snapshot_units->UnitTemperature_in_cgs);
+ message("to:");
+ message("(internal) Unit system: U_M = %e g.",
+ internal_units->UnitMass_in_cgs);
+ message("(internal) Unit system: U_L = %e cm.",
+ internal_units->UnitLength_in_cgs);
+ message("(internal) Unit system: U_t = %e s.",
+ internal_units->UnitTime_in_cgs);
+ message("(internal) Unit system: U_I = %e A.",
+ internal_units->UnitCurrent_in_cgs);
+ message("(internal) Unit system: U_T = %e K.",
+ internal_units->UnitTemperature_in_cgs);
+ }
+ }
+
/* Loop over all particle types */
for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ptype++) {
diff --git a/src/runner.c b/src/runner.c
index a0b8c51cfabf830afd685a210bd78d3f98af62f0..ae8a33f8b6f62ad0e89e78bc04f60c965c88cd93 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -119,7 +119,7 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
for (int i = 0; i < gcount; i++) {
/* Get a direct pointer on the part. */
- struct gpart *const g = &gparts[i];
+ struct gpart *restrict g = &gparts[i];
/* Is this part within the time step? */
if (g->ti_end <= ti_current) {
@@ -371,8 +371,8 @@ void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
*/
void runner_do_init(struct runner *r, struct cell *c, int timer) {
- struct part *const parts = c->parts;
- struct gpart *const gparts = c->gparts;
+ struct part *restrict parts = c->parts;
+ struct gpart *restrict gparts = c->gparts;
const int count = c->count;
const int gcount = c->gcount;
const int ti_current = r->e->ti_current;
@@ -390,7 +390,7 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
for (int i = 0; i < count; i++) {
/* Get a direct pointer on the part. */
- struct part *const p = &parts[i];
+ struct part *restrict p = &parts[i];
if (p->ti_end <= ti_current) {
@@ -403,7 +403,7 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
for (int i = 0; i < gcount; i++) {
/* Get a direct pointer on the part. */
- struct gpart *const gp = &gparts[i];
+ struct gpart *restrict gp = &gparts[i];
if (gp->ti_end <= ti_current) {
@@ -588,9 +588,9 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) {
const double dt = (r->e->ti_current - r->e->ti_old) * timeBase;
const int ti_old = r->e->ti_old;
const int ti_current = r->e->ti_current;
- struct part *const parts = c->parts;
- struct xpart *const xparts = c->xparts;
- struct gpart *const gparts = c->gparts;
+ struct part *restrict parts = c->parts;
+ struct xpart *restrict xparts = c->xparts;
+ struct gpart *restrict gparts = c->gparts;
float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, mass = 0.0;
@@ -611,7 +611,7 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) {
for (size_t k = 0; k < nr_gparts; k++) {
/* Get a handle on the gpart. */
- struct gpart *const gp = &gparts[k];
+ struct gpart *restrict gp = &gparts[k];
/* Drift... */
drift_gpart(gp, dt, timeBase, ti_old, ti_current);
@@ -628,8 +628,8 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) {
for (size_t k = 0; k < nr_parts; k++) {
/* Get a handle on the part. */
- struct part *const p = &parts[k];
- struct xpart *const xp = &xparts[k];
+ struct part *restrict p = &parts[k];
+ struct xpart *restrict xp = &xparts[k];
/* Drift... */
drift_part(p, xp, dt, timeBase, ti_old, ti_current);
@@ -684,7 +684,7 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) {
if (c->progeny[k] != NULL) {
/* Recurse */
- struct cell *cp = c->progeny[k];
+ struct cell *restrict cp = c->progeny[k];
runner_do_drift(r, cp, 0);
/* Collect */
@@ -734,9 +734,9 @@ void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) {
const double timeBase = r->e->timeBase;
const int count = c->count;
const int gcount = c->gcount;
- struct part *const parts = c->parts;
- struct xpart *const xparts = c->xparts;
- struct gpart *const gparts = c->gparts;
+ struct part *restrict parts = c->parts;
+ struct xpart *restrict xparts = c->xparts;
+ struct gpart *restrict gparts = c->gparts;
int updated = 0, g_updated = 0;
int ti_end_min = max_nr_timesteps, ti_end_max = 0;
@@ -757,7 +757,7 @@ void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) {
for (int k = 0; k < gcount; k++) {
/* Get a handle on the part. */
- struct gpart *const gp = &gparts[k];
+ struct gpart *restrict gp = &gparts[k];
/* If the g-particle has no counterpart */
if (gp->id_or_neg_offset > 0) {
@@ -783,8 +783,8 @@ void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) {
for (int k = 0; k < count; k++) {
/* Get a handle on the part. */
- struct part *const p = &parts[k];
- struct xpart *const xp = &xparts[k];
+ struct part *restrict p = &parts[k];
+ struct xpart *restrict xp = &xparts[k];
/* First, finish the force loop */
p->h_dt *= p->h * 0.333333333f;
@@ -812,7 +812,7 @@ void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) {
/* Loop over the progeny. */
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) {
- struct cell *const cp = c->progeny[k];
+ struct cell *restrict cp = c->progeny[k];
/* Recurse */
runner_do_kick_fixdt(r, cp, 0);
@@ -849,9 +849,9 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
const int ti_current = r->e->ti_current;
const int count = c->count;
const int gcount = c->gcount;
- struct part *const parts = c->parts;
- struct xpart *const xparts = c->xparts;
- struct gpart *const gparts = c->gparts;
+ struct part *restrict parts = c->parts;
+ struct xpart *restrict xparts = c->xparts;
+ struct gpart *restrict gparts = c->gparts;
int updated = 0, g_updated = 0;
int ti_end_min = max_nr_timesteps, ti_end_max = 0;
@@ -869,7 +869,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
for (int k = 0; k < gcount; k++) {
/* Get a handle on the part. */
- struct gpart *const gp = &gparts[k];
+ struct gpart *restrict gp = &gparts[k];
/* If the g-particle has no counterpart and needs to be kicked */
if (gp->id_or_neg_offset > 0) {
@@ -901,8 +901,8 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
for (int k = 0; k < count; k++) {
/* Get a handle on the part. */
- struct part *const p = &parts[k];
- struct xpart *const xp = &xparts[k];
+ struct part *restrict p = &parts[k];
+ struct xpart *restrict xp = &xparts[k];
/* If particle needs to be kicked */
if (p->ti_end <= ti_current) {
@@ -937,7 +937,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
/* Loop over the progeny. */
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) {
- struct cell *const cp = c->progeny[k];
+ struct cell *restrict cp = c->progeny[k];
/* Recurse */
runner_do_kick(r, cp, 0);
@@ -968,8 +968,8 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
*/
void runner_do_recv_cell(struct runner *r, struct cell *c, int timer) {
- const struct part *const parts = c->parts;
- const struct gpart *const gparts = c->gparts;
+ const struct part *restrict parts = c->parts;
+ const struct gpart *restrict gparts = c->gparts;
const size_t nr_parts = c->count;
const size_t nr_gparts = c->gcount;
// const int ti_current = r->e->ti_current;
diff --git a/src/serial_io.c b/src/serial_io.c
index 7bb829174c04545913a7ea560af1b4993526a0c5..0134fa3e52cd0df7341cbc88f11f4182b1c02949 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -39,9 +39,9 @@
#include "common_io.h"
#include "engine.h"
#include "error.h"
-#include "hydro_properties.h"
#include "gravity_io.h"
#include "hydro_io.h"
+#include "hydro_properties.h"
#include "io_properties.h"
#include "kernel_hydro.h"
#include "part.h"
@@ -764,6 +764,39 @@ void write_output_serial(struct engine* e, const char* baseName,
/* Print the system of Units used internally */
writeUnitSystem(h_file, internal_units, "InternalCodeUnits");
+ /* Tell the user if a conversion will be needed */
+ if (e->verbose) {
+ if (units_are_equal(snapshot_units, internal_units)) {
+
+ message("Snapshot and internal units match. No conversion needed.");
+
+ } else {
+
+ message("Conversion needed from:");
+ message("(Snapshot) Unit system: U_M = %e g.",
+ snapshot_units->UnitMass_in_cgs);
+ message("(Snapshot) Unit system: U_L = %e cm.",
+ snapshot_units->UnitLength_in_cgs);
+ message("(Snapshot) Unit system: U_t = %e s.",
+ snapshot_units->UnitTime_in_cgs);
+ message("(Snapshot) Unit system: U_I = %e A.",
+ snapshot_units->UnitCurrent_in_cgs);
+ message("(Snapshot) Unit system: U_T = %e K.",
+ snapshot_units->UnitTemperature_in_cgs);
+ message("to:");
+ message("(internal) Unit system: U_M = %e g.",
+ internal_units->UnitMass_in_cgs);
+ message("(internal) Unit system: U_L = %e cm.",
+ internal_units->UnitLength_in_cgs);
+ message("(internal) Unit system: U_t = %e s.",
+ internal_units->UnitTime_in_cgs);
+ message("(internal) Unit system: U_I = %e A.",
+ internal_units->UnitCurrent_in_cgs);
+ message("(internal) Unit system: U_T = %e K.",
+ internal_units->UnitTemperature_in_cgs);
+ }
+ }
+
/* Loop over all particle types */
for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ptype++) {
diff --git a/src/single_io.c b/src/single_io.c
index de7b96812365e136a078645511328338f8c4e8e3..3f634737194bf8bfdc780e662262b228ebc5c737 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -38,9 +38,9 @@
#include "common_io.h"
#include "engine.h"
#include "error.h"
-#include "hydro_properties.h"
#include "gravity_io.h"
#include "hydro_io.h"
+#include "hydro_properties.h"
#include "io_properties.h"
#include "kernel_hydro.h"
#include "part.h"
@@ -610,6 +610,39 @@ void write_output_single(struct engine* e, const char* baseName,
/* Print the system of Units used internally */
writeUnitSystem(h_file, internal_units, "InternalCodeUnits");
+ /* Tell the user if a conversion will be needed */
+ if (e->verbose) {
+ if (units_are_equal(snapshot_units, internal_units)) {
+
+ message("Snapshot and internal units match. No conversion needed.");
+
+ } else {
+
+ message("Conversion needed from:");
+ message("(Snapshot) Unit system: U_M = %e g.",
+ snapshot_units->UnitMass_in_cgs);
+ message("(Snapshot) Unit system: U_L = %e cm.",
+ snapshot_units->UnitLength_in_cgs);
+ message("(Snapshot) Unit system: U_t = %e s.",
+ snapshot_units->UnitTime_in_cgs);
+ message("(Snapshot) Unit system: U_I = %e A.",
+ snapshot_units->UnitCurrent_in_cgs);
+ message("(Snapshot) Unit system: U_T = %e K.",
+ snapshot_units->UnitTemperature_in_cgs);
+ message("to:");
+ message("(internal) Unit system: U_M = %e g.",
+ internal_units->UnitMass_in_cgs);
+ message("(internal) Unit system: U_L = %e cm.",
+ internal_units->UnitLength_in_cgs);
+ message("(internal) Unit system: U_t = %e s.",
+ internal_units->UnitTime_in_cgs);
+ message("(internal) Unit system: U_I = %e A.",
+ internal_units->UnitCurrent_in_cgs);
+ message("(internal) Unit system: U_T = %e K.",
+ internal_units->UnitTemperature_in_cgs);
+ }
+ }
+
/* Loop over all particle types */
for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ptype++) {
diff --git a/src/timestep.h b/src/timestep.h
index 1674958812d1f9e72eb618a9cb026befd610d06d..b4da4cdbaf5b082bd5c5890ca6305dd1dd952781 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -66,7 +66,7 @@ __attribute__((always_inline)) INLINE static int get_integer_timestep(
* @param e The #engine (used to get some constants).
*/
__attribute__((always_inline)) INLINE static int get_gpart_timestep(
- const struct gpart *gp, const struct engine *e) {
+ const struct gpart *restrict gp, const struct engine *restrict e) {
const float new_dt_external = gravity_compute_timestep_external(
e->external_potential, e->physical_constants, gp);
@@ -94,7 +94,8 @@ __attribute__((always_inline)) INLINE static int get_gpart_timestep(
* @param e The #engine (used to get some constants).
*/
__attribute__((always_inline)) INLINE static int get_part_timestep(
- const struct part *p, const struct xpart *xp, const struct engine *e) {
+ const struct part *restrict p, const struct xpart *restrict xp,
+ const struct engine *restrict e) {
/* Compute the next timestep (hydro condition) */
const float new_dt_hydro = hydro_compute_timestep(p, xp, e->hydro_properties);
diff --git a/src/units.c b/src/units.c
index 0284ebfb993c9c31d399a6ae45ec2b8d543e2a34..8fa4abfc618dae669bbc7be566a6db085d7be047 100644
--- a/src/units.c
+++ b/src/units.c
@@ -41,6 +41,20 @@
#include "const.h"
#include "error.h"
+/**
+ * @brief Initialises the UnitSystem structure with CGS system
+ *
+ * @param us The UnitSystem to initialize
+ */
+void units_init_cgs(struct UnitSystem* us) {
+
+ us->UnitMass_in_cgs = 1.;
+ us->UnitLength_in_cgs = 1.;
+ us->UnitTime_in_cgs = 1.;
+ us->UnitCurrent_in_cgs = 1.;
+ us->UnitTemperature_in_cgs = 1.;
+}
+
/**
* @brief Initialises the UnitSystem structure with the constants given in
* rhe parameter file.
@@ -66,6 +80,41 @@ void units_init(struct UnitSystem* us, const struct swift_params* params,
us->UnitTemperature_in_cgs = parser_get_param_double(params, buffer);
}
+/**
+ * @brief Initialises the UnitSystem structure with the constants given in
+ * rhe parameter file. Uses a default if the values are not present in the file.
+ *
+ * @param us The UnitSystem to initialize.
+ * @param params The parsed parameter file.
+ * @param category The section of the parameter file to read from.
+ * @param def The default unit system to copy from if required.
+ */
+void units_init_default(struct UnitSystem* us,
+ const struct swift_params* params, const char* category,
+ const struct UnitSystem* def) {
+
+ if (!def) error("Default UnitSystem not allocated");
+
+ char buffer[200];
+ sprintf(buffer, "%s:UnitMass_in_cgs", category);
+ us->UnitMass_in_cgs =
+ parser_get_opt_param_double(params, buffer, def->UnitMass_in_cgs);
+ sprintf(buffer, "%s:UnitLength_in_cgs", category);
+ us->UnitLength_in_cgs =
+ parser_get_opt_param_double(params, buffer, def->UnitLength_in_cgs);
+ sprintf(buffer, "%s:UnitVelocity_in_cgs", category);
+ const double defaultVelocity = def->UnitLength_in_cgs / def->UnitTime_in_cgs;
+ const double unitVelocity =
+ parser_get_opt_param_double(params, buffer, defaultVelocity);
+ us->UnitTime_in_cgs = us->UnitLength_in_cgs / unitVelocity;
+ sprintf(buffer, "%s:UnitCurrent_in_cgs", category);
+ us->UnitCurrent_in_cgs =
+ parser_get_opt_param_double(params, buffer, def->UnitCurrent_in_cgs);
+ sprintf(buffer, "%s:UnitTemp_in_cgs", category);
+ us->UnitTemperature_in_cgs =
+ parser_get_opt_param_double(params, buffer, def->UnitTemperature_in_cgs);
+}
+
/**
* @brief Returns the base unit conversion factor for a given unit system
* @param us The UnitSystem used
diff --git a/src/units.h b/src/units.h
index 29b4563163b4fe224c881b1c1055f2cbfbcc95a1..c67f6ebbab324e3c90ae86eb1ea11dcf013c5dc7 100644
--- a/src/units.h
+++ b/src/units.h
@@ -92,8 +92,13 @@ enum UnitConversionFactor {
UNIT_CONV_TEMPERATURE
};
+void units_init_cgs(struct UnitSystem*);
void units_init(struct UnitSystem*, const struct swift_params*,
const char* category);
+void units_init_default(struct UnitSystem* us,
+ const struct swift_params* params, const char* category,
+ const struct UnitSystem* def);
+
int units_are_equal(const struct UnitSystem* a, const struct UnitSystem* b);
/* Base units */
diff --git a/tests/makeInput.py b/tests/makeInput.py
index 8d8c07faf1e528869616bf93fbb7ffe8bf196af0..93c173b75aab3e14639372df5a77474acd3fc4e5 100644
--- a/tests/makeInput.py
+++ b/tests/makeInput.py
@@ -87,6 +87,15 @@ grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+
#Particle group
grp = file.create_group("/PartType0")
ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
diff --git a/tests/test27cells.c b/tests/test27cells.c
index b3bfdd18bd9165d68ecb380e31b9dee094c80f12..6f5de58d395c40e239b7da80fec22cd2bffb097f 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -183,10 +183,15 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
main_cell->parts[pid].v[2], main_cell->parts[pid].rho,
main_cell->parts[pid].rho_dh, main_cell->parts[pid].density.wcount,
main_cell->parts[pid].density.wcount_dh,
-#ifdef GADGET2_SPH
+#if defined(GADGET2_SPH)
main_cell->parts[pid].div_v, main_cell->parts[pid].density.rot_v[0],
main_cell->parts[pid].density.rot_v[1],
main_cell->parts[pid].density.rot_v[2]
+#elif defined(DEFAULT_SPH)
+ main_cell->parts[pid].density.div_v,
+ main_cell->parts[pid].density.rot_v[0],
+ main_cell->parts[pid].density.rot_v[1],
+ main_cell->parts[pid].density.rot_v[2]
#else
0., 0., 0., 0.
#endif
@@ -213,9 +218,12 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
cj->parts[pjd].x[2], cj->parts[pjd].v[0], cj->parts[pjd].v[1],
cj->parts[pjd].v[2], cj->parts[pjd].rho, cj->parts[pjd].rho_dh,
cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh,
-#ifdef GADGET2_SPH
+#if defined(GADGET2_SPH)
cj->parts[pjd].div_v, cj->parts[pjd].density.rot_v[0],
cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2]
+#elif defined(DEFAULT_SPH)
+ cj->parts[pjd].density.div_v, cj->parts[pjd].density.rot_v[0],
+ cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2]
#else
0., 0., 0., 0.
#endif
diff --git a/tests/testReading.c b/tests/testReading.c
index 600c69c47b615e9b3f3a4aaa7a85056e74976d4f..24b52371707cc42a34bcf6ae43cb82441ed3ee00 100644
--- a/tests/testReading.c
+++ b/tests/testReading.c
@@ -33,14 +33,18 @@ int main() {
struct part *parts = NULL;
struct gpart *gparts = NULL;
+ /* Default unit system */
+ struct UnitSystem us;
+ units_init_cgs(&us);
+
/* Properties of the ICs */
const double boxSize = 1.;
const int L = 4;
const double rho = 2.;
/* Read data */
- read_ic_single("input.hdf5", dim, &parts, &gparts, &Ngas, &Ngpart, &periodic,
- &flag_entropy_ICs, 0);
+ read_ic_single("input.hdf5", &us, dim, &parts, &gparts, &Ngas, &Ngpart,
+ &periodic, &flag_entropy_ICs, 0);
/* Check global properties read are correct */
assert(dim[0] == boxSize);
diff --git a/tests/testSPHStep.c b/tests/testSPHStep.c
index 12882347c2a7db1f8d89fd4b339722a8d0d43485..bca40fedb749e5d01cd7ef8a1fd2c77a6bb93b2a 100644
--- a/tests/testSPHStep.c
+++ b/tests/testSPHStep.c
@@ -96,9 +96,7 @@ void runner_dopair2_force(struct runner *r, struct cell *ci, struct cell *cj);
/* Run a full time step integration for one cell */
int main() {
-#ifndef DEFAULT_SPH
return 0;
-#else
int i, j, k, offset[3];
struct part *p;
@@ -131,7 +129,7 @@ int main() {
for (j = 0; j < 27; ++j)
for (i = 0; i < cells[j]->count; ++i) {
cells[j]->parts[i].mass = dim * dim * dim * rho / (N * N * N);
- cells[j]->parts[i].u = P / ((const_hydro_gamma - 1.) * rho);
+ cells[j]->parts[i].u = P / (hydro_gamma_minus_one * rho);
}
message("m=%f", dim * dim * dim * rho / (N * N * N));
@@ -199,6 +197,4 @@ int main() {
free(ci->xparts);
return 0;
-
-#endif
}
diff --git a/tests/testSingle.c b/tests/testSingle.c
index d37f20908a28fb2e1098043dd7fbbcf76bd8c247..6f6a8c8ad6cffbc8594d774c4a230220b8039890 100644
--- a/tests/testSingle.c
+++ b/tests/testSingle.c
@@ -85,10 +85,10 @@ int main(int argc, char *argv[]) {
p1.force.balsara = 0.0f;
p2.force.c = 58.8972740361f;
p2.force.balsara = 0.0f;
- p1.u = 1.e-5 / ((const_hydro_gamma - 1.) * p1.rho);
- p2.u = 1.e-5 / ((const_hydro_gamma - 1.) * p2.rho) + 100.0f / (33 * p2.mass);
- p1.force.POrho2 = p1.u * (const_hydro_gamma - 1.0f) / p1.rho;
- p2.force.POrho2 = p2.u * (const_hydro_gamma - 1.0f) / p2.rho;
+ p1.u = 1.e-5 / (hydro_gamma_minus_one * p1.rho);
+ p2.u = 1.e-5 / (hydro_gamma_minus_one * p2.rho) + 100.0f / (33 * p2.mass);
+ p1.force.POrho2 = p1.u * hydro_gamma_minus_one / p1.rho;
+ p2.force.POrho2 = p2.u * hydro_gamma_minus_one / p2.rho;
/* Dump a header. */
// printParticle_single(&p1, NULL);