diff --git a/examples/CoolingBox/coolingBox.yml b/examples/CoolingBox/coolingBox.yml
index c251edf54885e9490b9e07f7adb80fb51475be9f..781919a498e1cac7626c42d33352e0c507133a0a 100644
--- a/examples/CoolingBox/coolingBox.yml
+++ b/examples/CoolingBox/coolingBox.yml
@@ -9,8 +9,8 @@ InternalUnitSystem:
# Cosmological parameters
Cosmology:
h: 0.6777 # Reduced Hubble constant
- a_begin: 0.142857142857 # Initial scale-factor of the simulation
- a_end: 0.144285714286 # Final scale factor of the simulation
+ a_begin: 0.04 # Initial scale-factor of the simulation
+ a_end: 1.0 # Final scale factor of the simulation
Omega_m: 0.307 # Matter density parameter
Omega_lambda: 0.693 # Dark-energy density parameter
Omega_b: 0.0455 # Baryon density parameter
diff --git a/examples/CoolingBox/makeIC.py b/examples/CoolingBox/makeIC.py
index bd0311ec9c0b6cf23ad80df0804c69bf14361caf..a8e888526fb1f86dd1341af377c91a480f124ca5 100644
--- a/examples/CoolingBox/makeIC.py
+++ b/examples/CoolingBox/makeIC.py
@@ -27,7 +27,7 @@ from numpy import *
# Parameters
periodic= 1 # 1 For periodic box
boxSize = 3.0857e21 # 1 kiloparsec
-rho = 1.e-25 # Density in cgs
+rho = 2.36748e-25 # Density in cgs
P = 6.68e-12 # Pressure in cgs
gamma = 5./3. # Gas adiabatic index
eta = 1.2349 # 48 ngbs with cubic spline kernel
diff --git a/examples/CoolingBox/testCooling.c b/examples/CoolingBox/testCooling.c
index 2a5969e6464befa3681427175e94f356f8575dfd..37db87f2bcbbb493128c3b7c7923c70502c4ed63 100644
--- a/examples/CoolingBox/testCooling.c
+++ b/examples/CoolingBox/testCooling.c
@@ -94,7 +94,7 @@ int main(int argc, char **argv) {
/* Number of subcycles and tolerance used to compare
* subcycled and implicit solution. Note, high value
* of tolerance due to mismatch between explicit and
- * implicit solutioin for large timesteps */
+ * implicit solution for large timesteps */
const int n_subcycle = 1000;
const float integration_tolerance = 0.2;
diff --git a/examples/CoolingRates/Makefile.am b/examples/CoolingRates/Makefile.am
index ca9ddd8670957273b22cd8d4c6d91db782009f0e..058cdaf2efa3df3647af6f6e0263f65a0e515a15 100644
--- a/examples/CoolingRates/Makefile.am
+++ b/examples/CoolingRates/Makefile.am
@@ -23,10 +23,10 @@ AM_LDFLAGS = $(HDF5_LDFLAGS) $(HDF5_LIBS) $(FFTW_LIBS) $(TCMALLOC_LIBS) $(JEMALL
EXTRA_LIBS = $(HDF5_LIBS) $(FFTW_LIBS) $(PROFILER_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS) $(TBBMALLOC_LIBS) $(GRACKLE_LIBS) $(VELOCIRAPTOR_LIBS) $(GSL_LIBS)
# Programs.
-bin_PROGRAMS = testCooling
+bin_PROGRAMS = cooling_rates
# Sources
-testCooling_SOURCES = testCooling.c
-testCooling_CFLAGS = $(MYFLAGS) $(AM_CFLAGS)
-testCooling_LDADD = ../../src/.libs/libswiftsim.a $(EXTRA_LIBS)
+cooling_rates_SOURCES = cooling_rates.c
+cooling_rates_CFLAGS = $(MYFLAGS) $(AM_CFLAGS)
+cooling_rates_LDADD = ../../src/.libs/libswiftsim.a $(EXTRA_LIBS)
diff --git a/examples/CoolingRates/README b/examples/CoolingRates/README
index 414a4dc9ebb453cd0a1c59160a36139b437c025f..1beea8fafe1786862f77eff0c435414da7020d2f 100644
--- a/examples/CoolingRates/README
+++ b/examples/CoolingRates/README
@@ -1,6 +1,6 @@
This is a test that produces a plot of the contribution to the cooling rate from each of the elements depending on internal energy, density and redshift based on the EAGLE tables. To do so, the function in src/cooling/EAGLE returning the cooling rate is run for multiple values of the internal energy. The resulting cooling rates are written to files and plotted with a python script (cooling_rates_plot.py).
-The test may be run by using the shell script, test_script.sh, or by running:
+The test may be run by:
./testCooling -z X -d Y
python cooling_rates_plot.py
diff --git a/examples/CoolingRates/testCooling.c b/examples/CoolingRates/cooling_rates.c
similarity index 86%
rename from examples/CoolingRates/testCooling.c
rename to examples/CoolingRates/cooling_rates.c
index 87b07b1520bb2eb683820c9bee86f00554e9a188..31546c6a00f8e19bcf175361fe65ad81ad3cf620 100644
--- a/examples/CoolingRates/testCooling.c
+++ b/examples/CoolingRates/cooling_rates.c
@@ -44,23 +44,16 @@ void set_quantities(struct part *restrict p,
const struct phys_const *restrict internal_const, float nh,
double u) {
- float scale_factor = 1.0 / (1.0 + cosmo->z);
double hydrogen_number_density =
nh * pow(units_cgs_conversion_factor(us, UNIT_CONV_LENGTH), 3);
p->rho = hydrogen_number_density * internal_const->const_proton_mass /
p->chemistry_data.metal_mass_fraction[chemistry_element_H];
- float pressure = (u * scale_factor * scale_factor) / cooling->internal_energy_scale *
+ float pressure = (u * cosmo->a * cosmo->a) / cooling->internal_energy_scale *
p->rho * (hydro_gamma_minus_one);
p->entropy = pressure * (pow(p->rho, -hydro_gamma));
xp->entropy_full = p->entropy;
- // Using hydro_set_init_internal_energy seems to work better for higher z for
- // setting the internal energy correctly However, with Gadget2 this just sets
- // the entropy to the internal energy, which needs to be converted somehow
- if (cosmo->z >= 1)
- hydro_set_init_internal_energy(
- p, (u * scale_factor * scale_factor) / cooling->internal_energy_scale);
}
/*
@@ -79,18 +72,18 @@ int main(int argc, char **argv) {
struct phys_const internal_const;
struct cooling_function_data cooling;
struct cosmology cosmo;
- char *parametersFileName = "./testCooling.yml";
+ char *parametersFileName = "./cooling_rates.yml";
float nh; // hydrogen number density
double u; // internal energy
const int npts = 250; // number of values for the internal energy at which
// cooling rate is evaluated
+ // Set some default values
+ float redshift = 0.0, log_10_nh = -1;
+
// Read options
int param;
- float redshift = -1.0,
- log_10_nh =
- 100; // unreasonable values will be changed if not set in options
while ((param = getopt(argc, argv, "z:d:t")) != -1) switch (param) {
case 'z':
// read redshift
@@ -105,7 +98,7 @@ int main(int argc, char **argv) {
printf("Option -%c requires an argument.\n", optopt);
else
printf("Unknown option character `\\x%x'.\n", optopt);
- error("invalid option(s) to testCooling");
+ error("invalid option(s) to cooling_rates");
}
// Read the parameter file
@@ -125,17 +118,16 @@ int main(int argc, char **argv) {
// Init cosmology
cosmology_init(params, &us, &internal_const, &cosmo);
cosmology_print(&cosmo);
- if (redshift == -1.0) {
- cosmo.z = 3.0;
- } else {
- cosmo.z = redshift;
- }
+
+ // Set redshift and associated quantities
+ const float scale_factor = 1.0/(1.0+redshift);
+ integertime_t ti_current = log(scale_factor/cosmo.a_begin)/cosmo.time_base;
+ cosmology_update(&cosmo,&internal_const,ti_current);
// Init cooling
cooling_init(params, &us, &internal_const, &cooling);
cooling_print(&cooling);
cooling_update(&cosmo, &cooling, 0);
- message("redshift %.5e", cosmo.z);
// Calculate abundance ratios
float *abundance_ratio;
@@ -160,12 +152,8 @@ int main(int argc, char **argv) {
}
// set hydrogen number density
- if (log_10_nh == 100) {
- // hydrogen number density not specified in options
- nh = 1.0e-1;
- } else {
- nh = pow(10.0, log_10_nh);
- }
+ nh = exp(M_LN10 * log_10_nh);
+
// set internal energy to dummy value, will get reset when looping over
// internal energies
diff --git a/examples/CoolingRates/testCooling.yml b/examples/CoolingRates/cooling_rates.yml
similarity index 98%
rename from examples/CoolingRates/testCooling.yml
rename to examples/CoolingRates/cooling_rates.yml
index 63f3e381d73f706a2dcd97cf7882f3136eb2e35a..d90ed72764b67ee4e2f10b0993575a0c7ad8c8fe 100644
--- a/examples/CoolingRates/testCooling.yml
+++ b/examples/CoolingRates/cooling_rates.yml
@@ -9,7 +9,7 @@ InternalUnitSystem:
# Cosmological parameters
Cosmology:
h: 0.6777 # Reduced Hubble constant
- a_begin: 0.9090909 # Initial scale-factor of the simulation
+ a_begin: 0.04 # Initial scale-factor of the simulation
a_end: 1.0 # Final scale factor of the simulation
Omega_m: 0.307 # Matter density parameter
Omega_lambda: 0.693 # Dark-energy density parameter
diff --git a/examples/CoolingRates/test_script.sh b/examples/CoolingRates/test_script.sh
deleted file mode 100644
index 880dc02b8ee139c131d1bff6ed6f9d1d0a7e3015..0000000000000000000000000000000000000000
--- a/examples/CoolingRates/test_script.sh
+++ /dev/null
@@ -1,19 +0,0 @@
-#!/bin/bash
-
-#redshift_array=( 0 0.049 0.101 0.155 0.211 0.271 0.333 0.399 0.468 0.540)
-redshift_array=( 1 2 3 4 5 6 7 8 9 10 )
-
-for i in $(seq 10 -1 0); do
- redshift=${redshift_array[$i]}
- ./testCooling -z $redshift -d 3
- file='cooling_output_'$i'.dat'
- cp cooling_output.dat $file
-done
-
-#for nh in $(seq -1 -1 -6); do
-# ./testCooling -z 1 -d $nh
-# file='cooling_output_'$(expr 0 - $nh)'.dat'
-# cp cooling_output.dat $file
-#done
-
-python cooling_rates_plot.py
diff --git a/src/chemistry/EAGLE/chemistry.h b/src/chemistry/EAGLE/chemistry.h
index 1429b8887c4a8ce17fdb62ebcbaa94cbb3b1f3c3..6c14737b53ef7ecf7a38fc77981a1c3f0298560c 100644
--- a/src/chemistry/EAGLE/chemistry.h
+++ b/src/chemistry/EAGLE/chemistry.h
@@ -105,13 +105,13 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density(
cpd->smoothed_metal_mass_fraction[i] *= factor;
}
- // Smooth mass fraction of all metals
+ /* Smooth mass fraction of all metals */
cpd->smoothed_metal_mass_fraction_total +=
m * cpd->metal_mass_fraction_total * kernel_root;
cpd->smoothed_metal_mass_fraction_total *= factor;
- // Smooth iron mass fraction from SNIa
+ /* Smooth iron mass fraction from SNIa */
cpd->smoothed_iron_mass_fraction_from_SNIa +=
m * cpd->iron_mass_fraction_from_SNIa * kernel_root;
cpd->smoothed_iron_mass_fraction_from_SNIa *= factor;