diff --git a/examples/Cooling/CoolingBox/coolingBox.yml b/examples/Cooling/CoolingBox/coolingBox.yml
index ca97793aa9b4de43872db38f55c4060a2b27ddc9..f776823d1ba2b4703288913d4bf67f3e4851bebd 100644
--- a/examples/Cooling/CoolingBox/coolingBox.yml
+++ b/examples/Cooling/CoolingBox/coolingBox.yml
@@ -58,7 +58,7 @@ EAGLECooling:
dir_name: ./coolingtables/
H_reion_z: 11.5
H_reion_eV_p_H: 2.0
- He_reion_z_center: 3.5
+ He_reion_z_centre: 3.5
He_reion_z_sigma: 0.5
He_reion_eV_p_H: 2.0
diff --git a/src/cooling/EAGLE/cooling_rates.h b/src/cooling/EAGLE/cooling_rates.h
index 39978c9bfa8cf2fbea25e8856c2939a8d6ee69e4..d562f4adea9ff67f69da79669c25482807b70f57 100644
--- a/src/cooling/EAGLE/cooling_rates.h
+++ b/src/cooling/EAGLE/cooling_rates.h
@@ -503,8 +503,10 @@ INLINE static double eagle_metal_cooling_rate(
/* Sum up all the contributions */
double Lambda_net = Lambda_free + Lambda_Compton;
+ if (isnan(Lambda_net)) error("Lambda_net is nan free %.5e compton %.5e", Lambda_free, Lambda_Compton);
for (int elem = 2; elem < eagle_cooling_N_metal + 2; ++elem) {
Lambda_net += lambda_metal[elem];
+ if (isnan(Lambda_net)) error("Lambda_net is nan elem %d lambda_metal %.5e", elem, lambda_metal[elem]);
}
return Lambda_net;
diff --git a/tests/testCooling.c b/tests/testCooling.c
index 727a9638b09b871e866fe787438a5707fd43ec6b..dfcf22195a06235055f94719a5616b614d45b77b 100644
--- a/tests/testCooling.c
+++ b/tests/testCooling.c
@@ -20,8 +20,11 @@
/* Local headers. */
#include "swift.h"
+#include "../src/cooling/EAGLE/cooling_rates.h"
+#include "../src/cooling/EAGLE/interpolate.h"
+#include "../src/cooling/EAGLE/cooling_tables.h"
-#if 0
+#if 1
/*
* @brief Assign particle density and entropy corresponding to the
@@ -48,23 +51,21 @@ void set_quantities(struct part *restrict p, struct xpart *restrict xp,
cosmology_update(cosmo, phys_const, ti_current);
/* calculate density */
- double hydrogen_number_density = nh_cgs / cooling->number_density_scale;
+ double hydrogen_number_density = nh_cgs / cooling->number_density_to_cgs;
p->rho = hydrogen_number_density * phys_const->const_proton_mass /
p->chemistry_data.metal_mass_fraction[chemistry_element_H] *
(cosmo->a * cosmo->a * cosmo->a);
/* update entropy based on internal energy */
- float pressure = (u * cosmo->a * cosmo->a) / cooling->internal_energy_scale *
+ float pressure = (u * cosmo->a * cosmo->a) * cooling->internal_energy_from_cgs *
p->rho * (hydro_gamma_minus_one);
p->entropy = pressure * (pow(p->rho, -hydro_gamma));
xp->entropy_full = p->entropy;
}
/*
- * @brief Produces contributions to cooling rates for different
- * hydrogen number densities, from different metals,
- * tests 1d and 4d table interpolations produce
- * same results for cooling rate, dlambda/du and temperature.
+ * @brief Tests cooling integration scheme by comparing EAGLE
+ * integration to subcycled explicit equation.
*/
int main(int argc, char **argv) {
// Declare relevant structs
@@ -81,17 +82,23 @@ int main(int argc, char **argv) {
float nh; // hydrogen number density
double u; // internal energy
+ /* switch between checking comoving cooling and subcycling */
+ const int comoving_check = 1;
+
/* Number of values to test for in redshift,
* hydrogen number density and internal energy */
const int n_z = 50;
const int n_nh = 50;
const int n_u = 50;
+ //const int n_z = 2;
+ //const int n_nh = 2;
+ //const int n_u = 2;
/* Number of subcycles and tolerance used to compare
* subcycled and implicit solution. Note, high value
* of tolerance due to mismatch between explicit and
* implicit solution for large timesteps */
- const int n_subcycle = 1000;
+ const int n_subcycle = 10000;
const float integration_tolerance = 0.2;
/* Set dt */
@@ -110,17 +117,31 @@ int main(int argc, char **argv) {
/* Init chemistry */
chemistry_init(params, &us, &phys_const, &chem_data);
chemistry_first_init_part(&phys_const, &us, &cosmo, &chem_data, &p, &xp);
+ chemistry_part_has_no_neighbours(&p,&xp,&chem_data,&cosmo);
chemistry_print(&chem_data);
/* Init cosmology */
cosmology_init(params, &us, &phys_const, &cosmo);
cosmology_print(&cosmo);
+
+ /* Init hydro_props */
+ struct hydro_props hydro_properties;
+ hydro_props_init(&hydro_properties, &phys_const, &us, params);
/* Init cooling */
- cooling_init(params, &us, &phys_const, &cooling);
+ cooling_init(params, &us, &phys_const, &hydro_properties, &cooling);
cooling_print(&cooling);
cooling_update(&cosmo, &cooling, 0);
+ /* Init entropy floor */
+ struct entropy_floor_properties floor_props;
+ entropy_floor_init(&floor_props, &phys_const, &us, &hydro_properties, params);
+
+ /* Cooling function needs to know the minimal energy. Set it to the lowest
+ * internal energy in the cooling table. */
+ hydro_properties.minimal_internal_energy =
+ exp(M_LN10 * cooling.Therm[0]) * cooling.internal_energy_from_cgs;
+
/* Calculate abundance ratios */
float *abundance_ratio;
abundance_ratio = malloc((chemistry_element_count + 2) * sizeof(float));
@@ -133,61 +154,111 @@ int main(int argc, char **argv) {
(XH + p.chemistry_data.metal_mass_fraction[chemistry_element_He]);
int He_i;
float d_He;
- get_index_1d(cooling.HeFrac, cooling.N_He, HeFrac, &He_i, &d_He);
+ get_index_1d(cooling.HeFrac, eagle_cooling_N_He_frac, HeFrac, &He_i, &d_He);
- /* Cooling function needs to know the minimal energy. Set it to the lowest
- * internal energy in the cooling table. */
- struct hydro_props hydro_properties;
- hydro_properties.minimal_internal_energy =
- exp(M_LN10 * cooling.Therm[0]) / cooling.internal_energy_scale;
/* calculate spacing in nh and u */
- const float delta_nh = (cooling.nH[cooling.N_nH - 1] - cooling.nH[0]) / n_nh;
+ const float delta_nh = (cooling.nH[eagle_cooling_N_density - 1] - cooling.nH[0]) / n_nh;
const float delta_u =
- (cooling.Therm[cooling.N_Temp - 1] - cooling.Therm[0]) / n_u;
-
- for (int z_i = 0; z_i < n_z; z_i++) {
- integertime_t ti_current = max_nr_timesteps / n_z * z_i;
- for (int nh_i = 0; nh_i < n_nh; nh_i++) {
- nh = exp(M_LN10 * cooling.nH[0] + delta_nh * nh_i);
- for (int u_i = 0; u_i < n_u; u_i++) {
- u = exp(M_LN10 * cooling.Therm[0] + delta_u * u_i);
-
- /* update nh, u, z */
- set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const, nh, u,
- ti_current);
-
- /* calculate subcycled solution */
- for (int t_subcycle = 0; t_subcycle < n_subcycle; t_subcycle++) {
- p.entropy_dt = 0;
- cooling_cool_part(&phys_const, &us, &cosmo, &hydro_properties,
- &cooling, &p, &xp, dt_cool / n_subcycle,
- dt_therm / n_subcycle);
- xp.entropy_full += p.entropy_dt * dt_therm / n_subcycle;
- }
- double u_subcycled =
- hydro_get_physical_internal_energy(&p, &xp, &cosmo) *
- cooling.internal_energy_scale;
+ (cooling.Therm[eagle_cooling_N_temperature - 1] - cooling.Therm[0]) / n_u;
+
+ /* Declare variables we will be checking */
+ double u_implicit_cgs, u_check_cgs;
+ integertime_t ti_current;
- /* reset quantities to nh, u, and z that we want to test */
- set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const, nh, u,
- ti_current);
+ for (int nh_i = 0; nh_i < n_nh; nh_i++) {
+ nh = exp(M_LN10 * cooling.nH[0] + delta_nh * nh_i);
+ for (int u_i = 0; u_i < n_u; u_i++) {
+ u = exp(M_LN10 * cooling.Therm[0] + delta_u * u_i);
+ //u = exp(M_LN10 * cooling.Therm[0] + delta_u * u_i + 3.f);
- /* compute implicit solution */
- cooling_cool_part(&phys_const, &us, &cosmo, &hydro_properties, &cooling,
- &p, &xp, dt_cool, dt_therm);
- double u_implicit =
- hydro_get_physical_internal_energy(&p, &xp, &cosmo) *
- cooling.internal_energy_scale;
+ if (comoving_check) {
+ /* reset quantities to nh, u, and z that we want to test */
+ ti_current = 0;
+ set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const, nh, u,
+ ti_current);
+
+ /* compute implicit solution */
+ cooling_cool_part(&phys_const, &us, &cosmo, &hydro_properties, &floor_props,
+ &cooling, &p, &xp, dt_cool, dt_therm);
+ u_check_cgs =
+ hydro_get_physical_internal_energy(&p, &xp, &cosmo) *
+ cooling.internal_energy_to_cgs;
+ }
+ for (int z_i = 0; z_i < n_z; z_i++) {
+ ti_current = max_nr_timesteps / n_z * z_i;
+
+ if (!comoving_check) {
+ /* update nh, u, z */
+ set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const, nh, u,
+ ti_current);
+
+ /* calculate subcycled solution */
+ for (int t_subcycle = 0; t_subcycle < n_subcycle; t_subcycle++) {
+ p.entropy_dt = 0;
+ cooling_cool_part(&phys_const, &us, &cosmo, &hydro_properties,
+ &floor_props, &cooling, &p, &xp, dt_cool / n_subcycle,
+ dt_therm / n_subcycle);
+ xp.entropy_full += p.entropy_dt * dt_therm / n_subcycle;
+ }
+ u_check_cgs =
+ hydro_get_physical_internal_energy(&p, &xp, &cosmo) *
+ cooling.internal_energy_to_cgs;
+
+ /* reset quantities to nh, u, and z that we want to test */
+ set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const, nh, u,
+ ti_current);
+
+ /* compute implicit solution */
+ cooling_cool_part(&phys_const, &us, &cosmo, &hydro_properties, &floor_props,
+ &cooling, &p, &xp, dt_cool, dt_therm);
+ u_implicit_cgs =
+ hydro_get_physical_internal_energy(&p, &xp, &cosmo) *
+ cooling.internal_energy_to_cgs;
+
+ } else {
+ /* use set_quantities to set the redshift (and scalefactor) */
+ set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const, nh, u,
+ ti_current);
+ float density_1_physical = hydro_get_physical_density(&p, &cosmo);
+ float density_1_comoving = hydro_get_comoving_density(&p);
+ float internal_energy_1_physical = hydro_get_physical_internal_energy(&p, &xp, &cosmo);
+ float internal_energy_1_comoving = hydro_get_comoving_internal_energy(&p, &xp);
+
+ /* reset to get the comoving density */
+ set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const, nh * cosmo.a*cosmo.a*cosmo.a, u / cosmo.a2_inv,
+ ti_current);
+ float density_2_physical = hydro_get_physical_density(&p, &cosmo);
+ float density_2_comoving = hydro_get_comoving_density(&p);
+ float internal_energy_2_physical = hydro_get_physical_internal_energy(&p, &xp, &cosmo);
+ float internal_energy_2_comoving = hydro_get_comoving_internal_energy(&p, &xp);
+
+ message("scale factor %.5e density 1 2 physical comoving %.5e %.5e %.5e %.5e", cosmo.a, density_1_physical, density_1_comoving, density_2_physical, density_2_comoving);
+ message("scale factor %.5e internal_energy 1 2 physical comoving %.5e %.5e %.5e %.5e", cosmo.a, internal_energy_1_physical, internal_energy_1_comoving, internal_energy_2_physical, internal_energy_2_comoving);
+
+ /* compute implicit solution */
+ cooling_cool_part(&phys_const, &us, &cosmo, &hydro_properties, &floor_props,
+ &cooling, &p, &xp, dt_cool, dt_therm);
+ u_implicit_cgs =
+ hydro_get_physical_internal_energy(&p, &xp, &cosmo) *
+ cooling.internal_energy_to_cgs;
+ }
/* check if the two solutions are consistent */
- if (fabs((u_implicit - u_subcycled) / u_subcycled) >
- integration_tolerance)
+ if (fabs((u_implicit_cgs - u_check_cgs) / u_check_cgs) >
+ integration_tolerance) {
message(
- "implicit and subcycled solutions do not match. z_i %d nh_i %d "
- "u_i %d implicit %.5e subcycled %.5e error %.5e",
- z_i, nh_i, u_i, u_implicit, u_subcycled,
- fabs((u_implicit - u_subcycled) / u_subcycled));
+ "implicit and reference solutions do not match. z_i %d nh_i %d "
+ "u_i %d implicit %.5e reference %.5e error %.5e",
+ z_i, nh_i, u_i, u_implicit_cgs, u_check_cgs,
+ fabs((u_implicit_cgs - u_check_cgs) / u_check_cgs));
+ } else {
+ //message(
+ // "implicit and reference solutions match. z_i %d nh_i %d "
+ // "u_i %d implicit %.5e reference %.5e error %.5e",
+ // z_i, nh_i, u_i, u_implicit_cgs, u_check_cgs,
+ // fabs((u_implicit_cgs - u_check_cgs) / u_check_cgs));
+ }
}
}
}
diff --git a/tests/testCooling.yml b/tests/testCooling.yml
index faec32cdfec20b48af7341889c79b60bd2f6bb5b..9eea9b1041cdb7f8c65c1ec32a5ffe081a4b9a98 100644
--- a/tests/testCooling.yml
+++ b/tests/testCooling.yml
@@ -81,27 +81,51 @@ GrackleCooling:
MaxSteps: 1000
ConvergenceLimit: 1e-2
-EagleCooling:
- filename: /cosma5/data/Eagle/BG_Tables/CoolingTables/
- reionisation_redshift: 8.989
+EAGLECooling:
+ dir_name: ./coolingtables/
+ H_reion_z: 11.5
+ H_reion_eV_p_H: 2.0
He_reion_z_centre: 3.5
He_reion_z_sigma: 0.5
- He_reion_ev_pH: 2.0
+ He_reion_eV_p_H: 2.0
+
+#EAGLEChemistry:
+# InitMetallicity: 0.014
+# InitAbundance_Hydrogen: 0.70649785
+# InitAbundance_Helium: 0.28055534
+# InitAbundance_Carbon: 2.0665436e-3
+# InitAbundance_Nitrogen: 8.3562563e-4
+# InitAbundance_Oxygen: 5.4926244e-3
+# InitAbundance_Neon: 1.4144605e-3
+# InitAbundance_Magnesium: 5.907064e-4
+# InitAbundance_Silicon: 6.825874e-4
+# InitAbundance_Iron: 1.1032152e-3
+# CalciumOverSilicon: 0.0941736
+# SulphurOverSilicon: 0.6054160
+EAGLEChemistry: # Solar abundances
+ init_abundance_metal: 0.014
+ init_abundance_Hydrogen: 0.70649785
+ init_abundance_Helium: 0.28055534
+ init_abundance_Carbon: 2.0665436e-3
+ init_abundance_Nitrogen: 8.3562563e-4
+ init_abundance_Oxygen: 5.4926244e-3
+ init_abundance_Neon: 1.4144605e-3
+ init_abundance_Magnesium: 5.907064e-4
+ init_abundance_Silicon: 6.825874e-4
+ init_abundance_Iron: 1.1032152e-3
-EAGLEChemistry:
- InitMetallicity: 0.014
- InitAbundance_Hydrogen: 0.70649785
- InitAbundance_Helium: 0.28055534
- InitAbundance_Carbon: 2.0665436e-3
- InitAbundance_Nitrogen: 8.3562563e-4
- InitAbundance_Oxygen: 5.4926244e-3
- InitAbundance_Neon: 1.4144605e-3
- InitAbundance_Magnesium: 5.907064e-4
- InitAbundance_Silicon: 6.825874e-4
- InitAbundance_Iron: 1.1032152e-3
- CalciumOverSilicon: 0.0941736
- SulphurOverSilicon: 0.6054160
GearChemistry:
InitialMetallicity: 0.01295
+# Parameters for the EAGLE "equation of state"
+EAGLEEntropyFloor:
+ Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+ Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
+ Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
+ Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
+ Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+ Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
+ Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
+ Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
+