From 222c41275ba8fbe8abdbc31639a33c8df739e66b Mon Sep 17 00:00:00 2001
From: Alexei Borissov <dc-bori1@cosma-b.pri.cosma7.alces.network>
Date: Mon, 26 Nov 2018 16:21:29 +0000
Subject: [PATCH] Tidy CoolingRates example
---
examples/CoolingBox/coolingBox.yml | 4 +-
examples/CoolingBox/makeIC.py | 2 +-
examples/CoolingBox/testCooling.c | 2 +-
examples/CoolingRates/Makefile.am | 8 ++--
examples/CoolingRates/README | 2 +-
.../{testCooling.c => cooling_rates.c} | 38 +++++++------------
.../{testCooling.yml => cooling_rates.yml} | 2 +-
examples/CoolingRates/test_script.sh | 19 ----------
src/chemistry/EAGLE/chemistry.h | 4 +-
9 files changed, 25 insertions(+), 56 deletions(-)
rename examples/CoolingRates/{testCooling.c => cooling_rates.c} (86%)
rename examples/CoolingRates/{testCooling.yml => cooling_rates.yml} (98%)
delete mode 100644 examples/CoolingRates/test_script.sh
diff --git a/examples/CoolingBox/coolingBox.yml b/examples/CoolingBox/coolingBox.yml
index c251edf548..781919a498 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 bd0311ec9c..a8e888526f 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 2a5969e646..37db87f2bc 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 ca9ddd8670..058cdaf2ef 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 414a4dc9eb..1beea8fafe 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 87b07b1520..31546c6a00 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 63f3e381d7..d90ed72764 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 880dc02b8e..0000000000
--- 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 1429b8887c..6c14737b53 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;
--
GitLab