Tidy CoolingRates example

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

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

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;
 

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)
 

examples/CoolingRates/README

 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
 

examples/CoolingRates/testCooling.c → 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

examples/CoolingRates/testCooling.yml → 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

examples/CoolingRates/test_script.sh

-#!/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

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;