No cooling limit

configure.ac

@@ -1371,7 +1371,7 @@ esac
 # Hydro scheme.
 AC_ARG_WITH([hydro],
    [AS_HELP_STRING([--with-hydro=<scheme>],
-      [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, pressure-energy-monaghan, default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, anarchy-pu debug default: gadget2@:>@]
+      [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, pressure-energy-monaghan, default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, anarchy-pu default: gadget2@:>@]
    )],
    [with_hydro="$withval"],
    [with_hydro="gadget2"]

doc/Doxyfile.in

@@ -811,7 +811,7 @@ RECURSIVE              = NO
 # Note that relative paths are relative to the directory from which doxygen is
 # run.
 
-EXCLUDE                =
+EXCLUDE                = @top_srcdir@/src/cooling/EAGLE/newton_cooling.c
 
 # The EXCLUDE_SYMLINKS tag can be used to select whether or not files or
 # directories that are symbolic links (a Unix file system feature) are excluded

examples/Cooling/CoolingRates/cooling_rates.c

@@ -104,9 +104,9 @@ INLINE static double eagle_print_metal_cooling_rate(
 
   /* calculate cooling rates */
   for (int j = 0; j < eagle_cooling_N_metal + 2; j++) element_lambda[j] = 0.0;
-  lambda_net = eagle_metal_cooling_rate(
-      log10(u), cosmo->z, n_h, abundance_ratio, n_h_i, d_n_h, He_i, d_He,
-      cooling, /*dLambdaNet_du=*/NULL, element_lambda);
+  lambda_net =
+      eagle_metal_cooling_rate(log10(u), cosmo->z, n_h, abundance_ratio, n_h_i,
+                               d_n_h, He_i, d_He, cooling, element_lambda);
 
   /* write cooling rate contributions to their own files. */
   for (int j = 0; j < eagle_cooling_N_metal + 2; j++) {
@@ -235,18 +235,26 @@ int main(int argc, char **argv) {
 
   // Init cooling
   cooling_init(params, &us, &internal_const, &hydro_properties, &cooling);
+  cooling.H_reion_done = 1;
   cooling_print(&cooling);
   cooling_update(&cosmo, &cooling, &s);
 
+  // Copy over the raw metals into the smoothed metals
+  memcpy(&p.chemistry_data.smoothed_metal_mass_fraction,
+         &p.chemistry_data.metal_mass_fraction,
+         chemistry_element_count * sizeof(float));
+  p.chemistry_data.smoothed_metal_mass_fraction_total =
+      p.chemistry_data.metal_mass_fraction_total;
+
   // Calculate abundance ratios
   float abundance_ratio[(chemistry_element_count + 2)];
   abundance_ratio_to_solar(&p, &cooling, abundance_ratio);
 
   // extract mass fractions, calculate table indices and offsets
-  float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
-  float HeFrac =
-      p.chemistry_data.metal_mass_fraction[chemistry_element_He] /
-      (XH + p.chemistry_data.metal_mass_fraction[chemistry_element_He]);
+  float XH = p.chemistry_data.smoothed_metal_mass_fraction[chemistry_element_H];
+  float XHe =
+      p.chemistry_data.smoothed_metal_mass_fraction[chemistry_element_He];
+  float HeFrac = XHe / (XH + XHe);
   int He_i, n_h_i;
   float d_He, d_n_h;
   get_index_1d(cooling.HeFrac, eagle_cooling_N_He_frac, HeFrac, &He_i, &d_He);
@@ -286,7 +294,7 @@ int main(int argc, char **argv) {
 
     // calculate cooling rates
     const double temperature = eagle_convert_u_to_temp(
-        log10(u), cosmo.z, 0, NULL, n_h_i, He_i, d_n_h, d_He, &cooling);
+        log10(u), cosmo.z, n_h_i, He_i, d_n_h, d_He, &cooling);
 
     const double cooling_du_dt = eagle_print_metal_cooling_rate(
         n_h_i, d_n_h, He_i, d_He, &p, &xp, &cooling, &cosmo, &internal_const,

examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_high_z.yml

@@ -44,7 +44,20 @@ InitialConditions:
   periodic:   1
 
 # Parameters for the EAGLE chemistry
-EAGLEChemistry: 	     # Solar abundances
+# 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
+
+# Parameters for the EAGLE chemistry
+EAGLEChemistry: 	     # Primoridal abundances
   init_abundance_metal:      0.
   init_abundance_Hydrogen:   0.752
   init_abundance_Helium:     0.248
@@ -82,4 +95,4 @@ LambdaCooling:
 ConstCooling:
   cooling_rate: 1.          # Cooling rate (du/dt) (internal units)
   min_energy:   1.          # Minimal internal energy per unit mass (internal units)
-  cooling_tstep_mult: 1.    # Dimensionless pre-factor for the time-step condition
\ No newline at end of file
+  cooling_tstep_mult: 1.    # Dimensionless pre-factor for the time-step condition

examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml

@@ -44,7 +44,19 @@ InitialConditions:
   periodic:   1
 
 # Parameters for the EAGLE chemistry
-EAGLEChemistry: 	     # Primordial
+# 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: 	     # Primordial abundances
   init_abundance_metal:      0.
   init_abundance_Hydrogen:   0.752
   init_abundance_Helium:     0.248

examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml

@@ -35,7 +35,20 @@ InitialConditions:
   periodic:   1
 
 # Parameters for the EAGLE chemistry
-EAGLEChemistry: 	     # Solar abundances
+# 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
+
+# Parameters for the EAGLE chemistry
+EAGLEChemistry: 	     # Primoridal abundances
   init_abundance_metal:      0.
   init_abundance_Hydrogen:   0.752
   init_abundance_Helium:     0.248
@@ -73,4 +86,4 @@ LambdaCooling:
 ConstCooling:
   cooling_rate: 1.          # Cooling rate (du/dt) (internal units)
   min_energy:   1.          # Minimal internal energy per unit mass (internal units)
-  cooling_tstep_mult: 1.    # Dimensionless pre-factor for the time-step condition
\ No newline at end of file
+  cooling_tstep_mult: 1.    # Dimensionless pre-factor for the time-step condition

examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml

@@ -10,8 +10,8 @@ InternalUnitSystem:
 Gravity:
   eta:          0.025               # Constant dimensionless multiplier for time integration.
   theta:        0.7                 # Opening angle (Multipole acceptance criterion).
-  comoving_softening:     0.01      # Comoving softening length (in internal units).
-  max_physical_softening: 0.01      # Physical softening length (in internal units).
+  comoving_softening:     0.1       # Comoving softening length (in internal units).
+  max_physical_softening: 0.1       # Physical softening length (in internal units).
 
 # Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
 TimeIntegration:
@@ -69,16 +69,16 @@ EAGLEChemistry:
 
 # Hernquist potential parameters
 HernquistPotential:
-  useabspos:       0        # 0 -> positions based on centre, 1 -> absolute positions 
+  useabspos:       0       # 0 -> positions based on centre, 1 -> absolute positions 
   position:        [0.,0.,0.]    # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
-  idealizeddisk:   1        # Run with an idealized galaxy disk
+  idealizeddisk:   1       # Run with an idealized galaxy disk
   M200:            137.0   # M200 of the galaxy disk
-  h:               0.704    # reduced Hubble constant (value does not specify the used units!)
-  concentration:   9.0      # concentration of the Halo
-  diskfraction:              0.040   # Disk mass fraction
-  bulgefraction:              0.014   # Bulge mass fraction
-  timestep_mult:   0.01     # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
-  epsilon:         0.01      # Softening size (internal units)
+  h:               0.704   # reduced Hubble constant (value does not specify the used units!)
+  concentration:   9.0     # concentration of the Halo
+  diskfraction:    0.040   # Disk mass fraction
+  bulgefraction:   0.014   # Bulge mass fraction
+  timestep_mult:   0.01    # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
+  epsilon:         0.1     # Softening size (internal units)
  
 # EAGLE star formation parameters
 EAGLEStarFormation:

src/cell.c

@@ -53,6 +53,7 @@
 #include "chemistry.h"
 #include "drift.h"
 #include "engine.h"
+#include "entropy_floor.h"
 #include "error.h"
 #include "feedback.h"
 #include "gravity.h"
@@ -3595,7 +3596,8 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
 
       /* Drift... */
       drift_part(p, xp, dt_drift, dt_kick_hydro, dt_kick_grav, dt_therm,
-                 ti_old_part, ti_current);
+                 ti_old_part, ti_current, e->cosmology, e->hydro_properties,
+                 e->entropy_floor);
 
       /* Update the tracers properties */
       tracers_after_drift(p, xp, e->internal_units, e->physical_constants,

src/chemistry/EAGLE/chemistry.h

@@ -90,8 +90,9 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density(
 
   /* Some smoothing length multiples. */
   const float h = p->h;
-  const float h_inv = 1.0f / h;                       /* 1/h */
-  const float factor = pow_dimension(h_inv) / p->rho; /* 1 / h^d * rho */
+  const float h_inv = 1.0f / h; /* 1/h */
+  const float rho = hydro_get_comoving_density(p);
+  const float factor = pow_dimension(h_inv) / rho; /* 1 / (h^d * rho) */
   const float m = hydro_get_mass(p);
 
   struct chemistry_part_data* cpd = &p->chemistry_data;

src/cooling/EAGLE/cooling.c

@@ -48,18 +48,13 @@
 #include "space.h"
 #include "units.h"
 
-/* Maximum number of iterations for newton
- * and bisection integration schemes */
-static const int newton_max_iterations = 15;
+/* Maximum number of iterations for bisection scheme */
 static const int bisection_max_iterations = 150;
 
 /* Tolerances for termination criteria. */
 static const float explicit_tolerance = 0.05;
-static const float newton_tolerance = 1.0e-4;
 static const float bisection_tolerance = 1.0e-6;
-static const float rounding_tolerance = 1.0e-4;
 static const double bracket_factor = 1.5;
-static const double newton_log_u_guess_cgs = 12;
 
 /**
  * @brief Find the index of the current redshift along the redshift dimension
@@ -189,102 +184,6 @@ void cooling_update(const struct cosmology *cosmo,
   cooling->z_index = z_index;
 }
 
-/**
- * @brief Newton Raphson integration scheme to calculate particle cooling over
- * timestep. This replaces bisection scheme used in EAGLE to minimize the
- * number of array accesses. Integration defaults to bisection scheme (see
- * function bisection_iter) if this function does not converge within a
- * specified number of steps
- *
- * @param logu_init Initial guess for log(internal energy)
- * @param u_ini Internal energy at beginning of hydro step
- * @param n_H_index Particle hydrogen number density index
- * @param d_n_H Particle hydrogen number density offset
- * @param He_index Particle helium fraction index
- * @param d_He Particle helium fraction offset
- * @param He_reion_heat Heating due to helium reionization
- * (only depends on redshift, so passed as parameter)
- * @param p #part structure
- * @param cosmo #cosmology structure
- * @param cooling #cooling_function_data structure
- * @param phys_const #phys_const data structure
- * @param abundance_ratio Array of ratios of metal abundance to solar
- * @param dt timestep
- * @param bisection_flag Flag to identify if scheme failed to converge
- */
-INLINE static float newton_iter(
-    float logu_init, double u_ini, int n_H_index, float d_n_H, int He_index,
-    float d_He, float He_reion_heat, struct part *restrict p,
-    const struct cosmology *restrict cosmo,
-    const struct cooling_function_data *restrict cooling,
-    const struct phys_const *restrict phys_const,
-    const float abundance_ratio[chemistry_element_count + 2], float dt,
-    int *bisection_flag) {
-
-  double logu, logu_old;
-  double dLambdaNet_du = 0.0, LambdaNet;
-
-  /* table bounds */
-  const float log_table_bound_high =
-      (cooling->Therm[eagle_cooling_N_temperature - 1] - 0.05) / M_LOG10E;
-  const float log_table_bound_low = (cooling->Therm[0] + 0.05) / M_LOG10E;
-
-  /* convert Hydrogen mass fraction in Hydrogen number density */
-  const float XH =
-      p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_H];
-  const double n_H =
-      hydro_get_physical_density(p, cosmo) * XH / phys_const->const_proton_mass;
-  const double n_H_cgs = n_H * cooling->number_density_to_cgs;
-
-  /* compute ratefact = n_H * n_H / rho; Might lead to round-off error:
-   * replaced by equivalent expression below */
-  const double ratefact_cgs = n_H_cgs * XH * cooling->inv_proton_mass_cgs;
-
-  logu_old = logu_init;
-  logu = logu_old;
-  int i = 0;
-
-  float LambdaNet_old = 0;
-  LambdaNet = 0;
-  do /* iterate to convergence */
-  {
-    logu_old = logu;
-    LambdaNet_old = LambdaNet;
-    LambdaNet = (He_reion_heat / (dt * ratefact_cgs)) +
-                eagle_cooling_rate(logu_old, cosmo->z, n_H_cgs, abundance_ratio,
-                                   n_H_index, d_n_H, He_index, d_He, cooling,
-                                   &dLambdaNet_du);
-
-    /* Newton iteration. For details on how the cooling equation is integrated
-     * see documentation in theory/Cooling/ */
-    logu = logu_old - (1.0 - u_ini * exp(-logu_old) -
-                       LambdaNet * ratefact_cgs * dt * exp(-logu_old)) /
-                          (1.0 - dLambdaNet_du * ratefact_cgs * dt);
-    /* Check if first step passes over equilibrium solution, if it does adjust
-     * next guess */
-    if (i == 1 && LambdaNet_old * LambdaNet < 0) logu = newton_log_u_guess_cgs;
-
-    /* check whether iterations go within about 10% of the table bounds,
-     * if they do default to bisection method */
-    if (logu > log_table_bound_high) {
-      i = newton_max_iterations;
-      break;
-    } else if (logu < log_table_bound_low) {
-      i = newton_max_iterations;
-      break;
-    }
-
-    i++;
-  } while (fabs(logu - logu_old) > newton_tolerance &&
-           i < newton_max_iterations);
-  if (i >= newton_max_iterations) {
-    /* flag to trigger bisection scheme */
-    *bisection_flag = 1;
-  }
-
-  return logu;
-}
-
 /**
  * @brief Bisection integration scheme
  *
@@ -304,11 +203,12 @@ INLINE static float newton_iter(
  */
 INLINE static double bisection_iter(
     const double u_ini_cgs, const double n_H_cgs, const double redshift,
-    int n_H_index, float d_n_H, int He_index, float d_He,
-    double Lambda_He_reion_cgs, double ratefact_cgs,
+    const int n_H_index, const float d_n_H, const int He_index,
+    const float d_He, const double Lambda_He_reion_cgs,
+    const double ratefact_cgs,
     const struct cooling_function_data *restrict cooling,
-    const float abundance_ratio[chemistry_element_count + 2], double dt_cgs,
-    long long ID) {
+    const float abundance_ratio[chemistry_element_count + 2],
+    const double dt_cgs, const long long ID) {
 
   /* Bracketing */
   double u_lower_cgs = u_ini_cgs;
@@ -320,9 +220,8 @@ INLINE static double bisection_iter(
 
   double LambdaNet_cgs =
       Lambda_He_reion_cgs +
-      eagle_cooling_rate(log(u_ini_cgs), redshift, n_H_cgs, abundance_ratio,
-                         n_H_index, d_n_H, He_index, d_He, cooling,
-                         /*dLambdaNet_du=*/NULL);
+      eagle_cooling_rate(log10(u_ini_cgs), redshift, n_H_cgs, abundance_ratio,
+                         n_H_index, d_n_H, He_index, d_He, cooling);
 
   /*************************************/
   /* Let's try to bracket the solution */
@@ -335,11 +234,10 @@ INLINE static double bisection_iter(
     u_upper_cgs *= bracket_factor;
 
     /* Compute a new rate */
-    LambdaNet_cgs =
-        Lambda_He_reion_cgs +
-        eagle_cooling_rate(log(u_lower_cgs), redshift, n_H_cgs, abundance_ratio,
-                           n_H_index, d_n_H, He_index, d_He, cooling,
-                           /*dLambdaNet_du=*/NULL);
+    LambdaNet_cgs = Lambda_He_reion_cgs +
+                    eagle_cooling_rate(log10(u_lower_cgs), redshift, n_H_cgs,
+                                       abundance_ratio, n_H_index, d_n_H,
+                                       He_index, d_He, cooling);
 
     int i = 0;
     while (u_lower_cgs - u_ini_cgs - LambdaNet_cgs * ratefact_cgs * dt_cgs >
@@ -351,10 +249,9 @@ INLINE static double bisection_iter(
 
       /* Compute a new rate */
       LambdaNet_cgs = Lambda_He_reion_cgs +
-                      eagle_cooling_rate(log(u_lower_cgs), redshift, n_H_cgs,
+                      eagle_cooling_rate(log10(u_lower_cgs), redshift, n_H_cgs,
                                          abundance_ratio, n_H_index, d_n_H,
-                                         He_index, d_He, cooling,
-                                         /*dLambdaNet_du=*/NULL);
+                                         He_index, d_He, cooling);
       i++;
     }
 
@@ -371,11 +268,10 @@ INLINE static double bisection_iter(
     u_upper_cgs *= bracket_factor;
 
     /* Compute a new rate */
-    LambdaNet_cgs =
-        Lambda_He_reion_cgs +
-        eagle_cooling_rate(log(u_upper_cgs), redshift, n_H_cgs, abundance_ratio,
-                           n_H_index, d_n_H, He_index, d_He, cooling,
-                           /*dLambdaNet_du=*/NULL);
+    LambdaNet_cgs = Lambda_He_reion_cgs +
+                    eagle_cooling_rate(log10(u_upper_cgs), redshift, n_H_cgs,
+                                       abundance_ratio, n_H_index, d_n_H,
+                                       He_index, d_He, cooling);
 
     int i = 0;
     while (u_upper_cgs - u_ini_cgs - LambdaNet_cgs * ratefact_cgs * dt_cgs <
@@ -387,10 +283,9 @@ INLINE static double bisection_iter(
 
       /* Compute a new rate */
       LambdaNet_cgs = Lambda_He_reion_cgs +
-                      eagle_cooling_rate(log(u_upper_cgs), redshift, n_H_cgs,
+                      eagle_cooling_rate(log10(u_upper_cgs), redshift, n_H_cgs,
                                          abundance_ratio, n_H_index, d_n_H,
-                                         He_index, d_He, cooling,
-                                         /*dLambdaNet_du=*/NULL);
+                                         He_index, d_He, cooling);
       i++;
     }
 
@@ -417,11 +312,10 @@ INLINE static double bisection_iter(
     u_next_cgs = 0.5 * (u_lower_cgs + u_upper_cgs);
 
     /* New rate */
-    LambdaNet_cgs =
-        Lambda_He_reion_cgs +
-        eagle_cooling_rate(log(u_next_cgs), redshift, n_H_cgs, abundance_ratio,
-                           n_H_index, d_n_H, He_index, d_He, cooling,
-                           /*dLambdaNet_du=*/NULL);
+    LambdaNet_cgs = Lambda_He_reion_cgs +
+                    eagle_cooling_rate(log10(u_next_cgs), redshift, n_H_cgs,
+                                       abundance_ratio, n_H_index, d_n_H,
+                                       He_index, d_He, cooling);
 
     // Debugging
     if (u_next_cgs <= 0)
@@ -503,7 +397,7 @@ void cooling_cool_part(const struct phys_const *phys_const,
   /* Get the change in internal energy due to hydro forces */
   const float hydro_du_dt = hydro_get_physical_internal_energy_dt(p, cosmo);
 
-  /* Get internal energy at the end of the next kick step (assuming dt does not
+  /* Get internal energy at the end of the step (assuming dt does not
    * increase) */
   double u_0 = (u_start + hydro_du_dt * dt_therm);
 
@@ -511,7 +405,6 @@ void cooling_cool_part(const struct phys_const *phys_const,
   u_0 = max(u_0, hydro_properties->minimal_internal_energy);
 
   /* Convert to CGS units */
-  const double u_start_cgs = u_start * cooling->internal_energy_to_cgs;
   const double u_0_cgs = u_0 * cooling->internal_energy_to_cgs;
   const double dt_cgs = dt * units_cgs_conversion_factor(us, UNIT_CONV_TIME);
 
@@ -536,7 +429,7 @@ void cooling_cool_part(const struct phys_const *phys_const,
    * metal-free Helium mass fraction as per the Wiersma+08 definition */
   const float HeFrac = XHe / (XH + XHe);
 
-  /* convert Hydrogen mass fraction into Hydrogen number density */
+  /* convert Hydrogen mass fraction into physical Hydrogen number density */
   const double n_H =
       hydro_get_physical_density(p, cosmo) * XH / phys_const->const_proton_mass;
   const double n_H_cgs = n_H * cooling->number_density_to_cgs;
@@ -573,10 +466,9 @@ void cooling_cool_part(const struct phys_const *phys_const,
 
   /* First try an explicit integration (note we ignore the derivative) */
   const double LambdaNet_cgs =
-      Lambda_He_reion_cgs + eagle_cooling_rate(log(u_0_cgs), cosmo->z, n_H_cgs,
-                                               abundance_ratio, n_H_index,
-                                               d_n_H, He_index, d_He, cooling,
-                                               /*dLambdaNet_du=*/NULL);
+      Lambda_He_reion_cgs +
+      eagle_cooling_rate(log10(u_0_cgs), cosmo->z, n_H_cgs, abundance_ratio,
+                         n_H_index, d_n_H, He_index, d_He, cooling);
 
   /* if cooling rate is small, take the explicit solution */
   if (fabs(ratefact_cgs * LambdaNet_cgs * dt_cgs) <
@@ -586,76 +478,32 @@ void cooling_cool_part(const struct phys_const *phys_const,
 
   } else {
 
-    int bisection_flag = 1;
-
-    // MATTHIEU: TO DO restore the Newton-Raphson scheme
-    if (0 && cooling->newton_flag) {
-
-      /* Ok, try a Newton-Raphson scheme instead */
-      double log_u_final_cgs =
-          newton_iter(log(u_0_cgs), u_0_cgs, n_H_index, d_n_H, He_index, d_He,
-                      Lambda_He_reion_cgs, p, cosmo, cooling, phys_const,
-                      abundance_ratio, dt_cgs, &bisection_flag);
-
-      /* Check if newton scheme sent us to a higher energy despite being in
-         a  cooling regime If it did try newton scheme with a better guess.
-         (Guess internal energy near equilibrium solution).  */
-      if (LambdaNet_cgs < 0 && log_u_final_cgs > log(u_0_cgs)) {
-        bisection_flag = 0;
-        log_u_final_cgs =
-            newton_iter(newton_log_u_guess_cgs, u_0_cgs, n_H_index, d_n_H,
-                        He_index, d_He, Lambda_He_reion_cgs, p, cosmo, cooling,
-                        phys_const, abundance_ratio, dt_cgs, &bisection_flag);
-      }
-
-      u_final_cgs = exp(log_u_final_cgs);
-    }
-
-    /* Alright, all else failed, let's bisect */
-    if (bisection_flag || !(cooling->newton_flag)) {
-      u_final_cgs =
-          bisection_iter(u_0_cgs, n_H_cgs, cosmo->z, n_H_index, d_n_H, He_index,
-                         d_He, Lambda_He_reion_cgs, ratefact_cgs, cooling,
-                         abundance_ratio, dt_cgs, p->id);
-    }
+    /* Otherwise, go the bisection route. */
+    u_final_cgs =
+        bisection_iter(u_0_cgs, n_H_cgs, cosmo->z, n_H_index, d_n_H, He_index,
+                       d_He, Lambda_He_reion_cgs, ratefact_cgs, cooling,
+                       abundance_ratio, dt_cgs, p->id);
   }
 
-  /* Expected change in energy over the next kick step
-     (assuming no change in dt) */
-  const double delta_u_cgs = u_final_cgs - u_start_cgs;
-
   /* Convert back to internal units */
-  double delta_u = delta_u_cgs * cooling->internal_energy_from_cgs;
+  double u_final = u_final_cgs * cooling->internal_energy_from_cgs;
 
   /* We now need to check that we are not going to go below any of the limits */
 
-  /* Limit imposed by the entropy floor */
-  const double A_floor = entropy_floor(p, cosmo, floor_props);
-  const double rho = hydro_get_physical_density(p, cosmo);
-  const double u_floor = gas_internal_energy_from_entropy(rho, A_floor);
-
   /* Absolute minimum */
   const double u_minimal = hydro_properties->minimal_internal_energy;
+  u_final = max(u_final, u_minimal);
 
-  /* Largest of both limits */
-  const double u_limit = max(u_minimal, u_floor);
-
-  /* First, check whether we may end up below the minimal energy after
-   * this step 1/2 kick + another 1/2 kick that could potentially be for
-   * a time-step twice as big. We hence check for 1.5 delta_u. */
-  if (u_start + 1.5 * delta_u < u_limit) {
-    delta_u = (u_limit - u_start) / 1.5;
-  }
+  /* Limit imposed by the entropy floor */
+  const double A_floor = entropy_floor(p, cosmo, floor_props);
+  const double rho_physical = hydro_get_physical_density(p, cosmo);
+  const double u_floor =
+      gas_internal_energy_from_entropy(rho_physical, A_floor);
+  u_final = max(u_final, u_floor);
 
-  /* Second, check whether the energy used in the prediction could get negative.
-   * We need to check for the 1/2 dt kick followed by a full time-step drift
-   * that could potentially be for a time-step twice as big. We hence check
-   * for 2.5 delta_u but this time against 0 energy not the minimum.
-   * To avoid numerical rounding bringing us below 0., we add a tiny tolerance.
-   */
-  if (u_start + 2.5 * delta_u < 0.) {
-    delta_u = -u_start / (2.5 + rounding_tolerance);
-  }
+  /* Expected change in energy over the next kick step
+     (assuming no change in dt) */
+  const double delta_u = u_final - max(u_start, u_floor);
 
   /* Turn this into a rate of change (including cosmology term) */
   const float cooling_du_dt = delta_u / dt_therm;
@@ -774,8 +622,7 @@ float cooling_get_temperature(
 
   /* Compute the log10 of the temperature by interpolating the table */
   const double log_10_T = eagle_convert_u_to_temp(
-      log10(u_cgs), cosmo->z, /*compute_dT_du=*/0, /*dT_du=*/NULL, n_H_index,
-      He_index, d_n_H, d_He, cooling);
+      log10(u_cgs), cosmo->z, n_H_index, He_index, d_n_H, d_He, cooling);
 
   /* Undo the log! */
   return exp10(log_10_T);
@@ -940,10 +787,6 @@ void cooling_init_backend(struct swift_params *parameter_file,
 
   /* set previous_z_index and to last value of redshift table*/
   cooling->previous_z_index = eagle_cooling_N_redshifts - 2;
-
-  /* Check if we are running with the newton scheme */
-  cooling->newton_flag = parser_get_opt_param_int(
-      parameter_file, "EAGLECooling:newton_integration", 0);
 }
 
 /**

src/cooling/EAGLE/cooling_rates.h

@@ -99,7 +99,7 @@ __attribute__((always_inline)) INLINE void abundance_ratio_to_solar(
       cooling->Ca_over_Si_ratio_in_solar;
 
   ratio_solar[10] =
-      p->chemistry_data.metal_mass_fraction[chemistry_element_Fe] *
+      p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_Fe] *
       cooling->SolarAbundances_inv[10 /* Fe */];
 }
 
@@ -118,7 +118,8 @@ __attribute__((always_inline)) INLINE void abundance_ratio_to_solar(
  */
 __attribute__((always_inline)) INLINE double
 eagle_helium_reionization_extraheat(