update cooling test

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
 

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;

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));
+        }
       }
     }
   }

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
+