Update the lambda-cooling model with correct cosmological terms.

examples/SmallCosmoVolume/small_cosmo_volume.yml

@@ -56,3 +56,8 @@ InitialConditions:
   cleanup_velocity_factors:    1  
   generate_gas_in_ics: 1            # Generate gas particles from the DM-only ICs
   cleanup_smoothing_lengths: 1      # Since we generate gas, make use of the (expensive) cleaning-up procedure.
+
+# Constant lambda cooling function
+LambdaCooling:
+  lambda_cgs:                  1e-22 # Cooling rate (in cgs units)
+  cooling_tstep_mult:          1.0   # Dimensionless pre-factor for the time-step condition

examples/parameter_example.yml

@@ -198,7 +198,7 @@ ConstCooling:
 
 # Constant lambda cooling function
 LambdaCooling:
-  lambda:                      2.0   # Cooling rate (in cgs units)
+  lambda_cgs:                  2.0   # Cooling rate (in cgs units)
   minimum_temperature:         1.0e4 # Minimal temperature (Kelvin)
   mean_molecular_weight:       0.59  # Mean molecular weight
   hydrogen_mass_abundance:     0.75  # Hydrogen mass abundance (dimensionless)

src/cooling/const_lambda/cooling.h

@@ -40,30 +40,33 @@
 #include "units.h"
 
 /**
- * @brief Calculates du/dt in code units for a particle.
+ * @brief Calculates du/dt in CGS units for a particle.
  *
  * @param phys_const The physical constants in internal units.
  * @param us The internal system of units.
  * @param cosmo The current cosmological model.
+ * @param hydro_props The properties of the hydro scheme.
  * @param cooling The #cooling_function_data used in the run.
  * @param p Pointer to the particle data..
  */
-__attribute__((always_inline)) INLINE static float cooling_rate(
-    const struct phys_const* const phys_const, const struct unit_system* us,
+__attribute__((always_inline)) INLINE static double cooling_rate_cgs(
     const struct cosmology* restrict cosmo,
+    const struct hydro_props* hydro_props,
     const struct cooling_function_data* cooling, const struct part* p) {
 
   /* Get particle density */
-  const float rho = hydro_get_physical_density(p, cosmo);
+  const double rho = hydro_get_physical_density(p, cosmo);
+  const double rho_cgs = rho * cooling->conv_factor_cgs_density;
 
   /* Get cooling function properties */
-  const float X_H = cooling->hydrogen_mass_abundance;
+  const double X_H = hydro_props->hydrogen_mass_fraction;
 
   /* Calculate du_dt */
-  const float du_dt = -cooling->lambda *
-                      (X_H * rho / phys_const->const_proton_mass) *
-                      (X_H * rho / phys_const->const_proton_mass) / rho;
-  return du_dt;
+  const double du_dt_cgs = -cooling->lambda_cgs *
+                           (X_H * rho_cgs / cooling->proton_mass_cgs) *
+                           (X_H * rho_cgs / cooling->proton_mass_cgs) / rho_cgs;
+
+  return du_dt_cgs;
 }
 
 /**
@@ -80,29 +83,42 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
     const struct phys_const* restrict phys_const,
     const struct unit_system* restrict us,
     const struct cosmology* restrict cosmo,
+    const struct hydro_props* hydro_props,
     const struct cooling_function_data* restrict cooling,
-    struct part* restrict p, struct xpart* restrict xp, float dt) {
+    struct part* restrict p, struct xpart* restrict xp, const float dt,
+    const float dt_therm) {
+
+  if (dt == 0.) return;
 
   /* Internal energy floor */
-  const float u_floor = cooling->min_energy;
+  const double u_floor = hydro_props->minimal_internal_energy;
 
   /* Current energy */
-  const float u_old = hydro_get_physical_internal_energy(p, cosmo);
+  const double u_old = hydro_get_physical_internal_energy2(p, xp, cosmo);
 
-  /* Current du_dt */
-  const float hydro_du_dt = hydro_get_internal_energy_dt(p);
+  /* Current du_dt in physical coordinates */
+  const double hydro_du_dt = hydro_get_physical_internal_energy_dt(p, cosmo);
 
   /* Calculate cooling du_dt */
-  float cooling_du_dt = cooling_rate(phys_const, us, cosmo, cooling, p);
+  const double cooling_du_dt_cgs =
+      cooling_rate_cgs(cosmo, hydro_props, cooling, p);
+
+  /* Convert to internal units */
+  double cooling_du_dt =
+      cooling_du_dt_cgs / cooling->conv_factor_cgs_energy_rate;
 
   /* Integrate cooling equation to enforce energy floor */
   /* Factor of 1.5 included since timestep could potentially double */
-  if (u_old + (hydro_du_dt + cooling_du_dt) * 1.5f * dt < u_floor) {
-    cooling_du_dt = -(u_old + 1.5f * dt * hydro_du_dt - u_floor) / (1.5f * dt);
+  if (u_old + (hydro_du_dt * dt_therm + cooling_du_dt * dt) * 2.5f < u_floor) {
+    cooling_du_dt =
+        0.;  //-(u_old + 2.5f * dt * hydro_du_dt - u_floor) / (2.5f * dt);
   }
 
+  /* Add cosmological term */
+  cooling_du_dt *= cosmo->a * cosmo->a;
+
   /* Update the internal energy time derivative */
-  hydro_set_internal_energy_dt(p, hydro_du_dt + cooling_du_dt);
+  hydro_set_physical_internal_energy_dt(p, cosmo, hydro_du_dt + cooling_du_dt);
 
   /* Store the radiated energy */
   xp->cooling_data.radiated_energy += -hydro_get_mass(p) * cooling_du_dt * dt;
@@ -121,17 +137,24 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
     const struct cooling_function_data* restrict cooling,
     const struct phys_const* restrict phys_const,
     const struct cosmology* restrict cosmo,
-    const struct unit_system* restrict us, const struct part* restrict p) {
-
-  /* Get current internal energy */
-  const float u = hydro_get_physical_internal_energy(p, cosmo);
-  const float du_dt = cooling_rate(phys_const, us, cosmo, cooling, p);
-
-  /* If we are close to (or below) the energy floor, we ignore the condition */
-  if (u < 1.01f * cooling->min_energy)
-    return FLT_MAX;
-  else
-    return cooling->cooling_tstep_mult * u / fabsf(du_dt);
+    const struct unit_system* restrict us,
+    const struct hydro_props* hydro_props, const struct part* restrict p) {
+
+  /* /\* Get current internal energy *\/ */
+  /* const float u = hydro_get_physical_internal_energy(p, cosmo); */
+  /* float cooling_du_dt = */
+  /*     cooling_rate(phys_const, us, cosmo, hydro_props, cooling, p); */
+
+  /* /\* Add cosmological term *\/ */
+  /* cooling_du_dt *= cosmo->a * cosmo->a; */
+
+  /* /\* If we are close to (or below) the energy floor, we ignore the condition
+   * *\/ */
+  /* if (u < 1.01f * hydro_props->minimal_internal_energy) */
+  /*   return FLT_MAX; */
+  /* else */
+  /*   return cooling->cooling_tstep_mult * u / fabsf(cooling_du_dt); */
+  return FLT_MAX;
 }
 
 /**
@@ -176,30 +199,24 @@ static INLINE void cooling_init_backend(struct swift_params* parameter_file,
                                         const struct phys_const* phys_const,
                                         struct cooling_function_data* cooling) {
 
-  const double lambda_cgs =
+  /* Read in the cooling parameters */
+  cooling->lambda_cgs =
       parser_get_param_double(parameter_file, "LambdaCooling:lambda_cgs");
-  const float min_temperature = parser_get_param_double(
-      parameter_file, "LambdaCooling:minimum_temperature");
-  cooling->hydrogen_mass_abundance = parser_get_param_double(
-      parameter_file, "LambdaCooling:hydrogen_mass_abundance");
-  cooling->mean_molecular_weight = parser_get_param_double(
-      parameter_file, "LambdaCooling:mean_molecular_weight");
   cooling->cooling_tstep_mult = parser_get_param_double(
       parameter_file, "LambdaCooling:cooling_tstep_mult");
 
-  /* convert minimum temperature into minimum internal energy */
-  const float u_floor =
-      phys_const->const_boltzmann_k * min_temperature /
-      (hydro_gamma_minus_one * cooling->mean_molecular_weight *
-       phys_const->const_proton_mass);
-
-  cooling->min_energy = u_floor;
-
-  /* convert lambda to code units */
-  cooling->lambda = lambda_cgs *
-                    units_cgs_conversion_factor(us, UNIT_CONV_TIME) /
-                    (units_cgs_conversion_factor(us, UNIT_CONV_ENERGY) *
-                     units_cgs_conversion_factor(us, UNIT_CONV_VOLUME));
+  /* Some useful conversion values */
+  cooling->conv_factor_cgs_density =
+      units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
+  cooling->conv_factor_cgs_energy =
+      units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+  cooling->conv_factor_cgs_energy_rate =
+      units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS) /
+      units_cgs_conversion_factor(us, UNIT_CONV_TIME);
+
+  /* Useful constants */
+  cooling->proton_mass_cgs = phys_const->const_proton_mass *
+                             units_cgs_conversion_factor(us, UNIT_CONV_MASS);
 }
 
 /**
@@ -210,12 +227,10 @@ static INLINE void cooling_init_backend(struct swift_params* parameter_file,
 static INLINE void cooling_print_backend(
     const struct cooling_function_data* cooling) {
 
-  message(
-      "Cooling function is 'Constant lambda' with "
-      "(lambda,min_energy,hydrogen_mass_abundance,mean_molecular_weight) "
-      "=  (%g,%g,%g,%g)",
-      cooling->lambda, cooling->min_energy, cooling->hydrogen_mass_abundance,
-      cooling->mean_molecular_weight);
+  message("Cooling function is 'Constant lambda' with Lambda=%g [cgs]",
+          cooling->lambda_cgs);
+
+  message("%e", cooling->proton_mass_cgs);
 }
 
 #endif /* SWIFT_COOLING_CONST_LAMBDA_H */

src/cooling/const_lambda/cooling_struct.h

@@ -28,17 +28,21 @@
  */
 struct cooling_function_data {
 
-  /*! Cooling rate in internal units */
-  double lambda;
+  /*! Cooling rate in cgs units */
+  double lambda_cgs;
 
-  /*! Fraction of gas mass that is Hydrogen. Used to calculate n_H */
-  float hydrogen_mass_abundance;
+  /*! Conversion factor from internal units to cgs for density */
+  double conv_factor_cgs_density;
 
-  /*! 'mu', used to convert min_temperature to min_internal energy */
-  float mean_molecular_weight;
+  /*! Conversion factor from internal units to cgs for internal energy */
+  double conv_factor_cgs_energy;
 
-  /*! Minimally allowed internal energy of all the particles */
-  float min_energy;
+  /*! Conversion factor from internal units to cgs for internal energy
+   * derivative */
+  double conv_factor_cgs_energy_rate;
+
+  /*! Proton mass in cgs units */
+  double proton_mass_cgs;
 
   /*! Constant multiplication factor for time-step criterion */
   float cooling_tstep_mult;

src/hydro/Gadget2/hydro.h

@@ -65,6 +65,21 @@ hydro_get_physical_internal_energy(const struct part *restrict p,
   return gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, p->entropy);
 }
 
+/**
+ * @brief Returns the physical internal energy of a particle
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_internal_energy2(const struct part *restrict p,
+                                    const struct xpart *restrict xp,
+                                    const struct cosmology *cosmo) {
+
+  return gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv,
+                                          xp->entropy_full);
+}
+
 /**
  * @brief Returns the comoving pressure of a particle
  *
@@ -204,32 +219,66 @@ __attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
 }
 
 /**
- * @brief Returns the time derivative of internal energy of a particle
+ * @brief Returns the time derivative of co-moving internal energy of a particle
  *
  * We assume a constant density.
  *
  * @param p The particle of interest
  */
-__attribute__((always_inline)) INLINE static float hydro_get_internal_energy_dt(
-    const struct part *restrict p) {
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_internal_energy_dt(const struct part *restrict p) {
 
   return gas_internal_energy_from_entropy(p->rho, p->entropy_dt);
 }
 
 /**
- * @brief Returns the time derivative of internal energy of a particle
+ * @brief Returns the time derivative of physical internal energy of a particle
  *
  * We assume a constant density.
  *
  * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_internal_energy_dt(const struct part *restrict p,
+                                      const struct cosmology *cosmo) {
+
+  return gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv,
+                                          p->entropy_dt);
+}
+
+/**
+ * @brief Sets the time derivative of the co-moving internal energy of a
+ * particle
+ *
+ * We assume a constant density for the conversion to entropy.
+ *
+ * @param p The particle of interest.
  * @param du_dt The new time derivative of the internal energy.
  */
-__attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
-    struct part *restrict p, float du_dt) {
+__attribute__((always_inline)) INLINE static void
+hydro_set_comoving_internal_energy_dt(struct part *restrict p, float du_dt) {
 
   p->entropy_dt = gas_entropy_from_internal_energy(p->rho, du_dt);
 }
 
+/**
+ * @brief Sets the time derivative of the physical internal energy of a particle
+ *
+ * We assume a constant density for the conversion to entropy.
+ *
+ * @param p The particle of interest.
+ * @param cosmo Cosmology data structure
+ * @param du_dt The time derivative of the internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_physical_internal_energy_dt(struct part *restrict p,
+                                      const struct cosmology *restrict cosmo,
+                                      float du_dt) {
+  p->entropy_dt =
+      gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, du_dt);
+}
+
 /**
  * @brief Computes the hydro time-step of a given particle
  *

src/runner.c

@@ -433,6 +433,7 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
   const struct cooling_function_data *cooling_func = e->cooling_func;
   const struct phys_const *constants = e->physical_constants;
   const struct unit_system *us = e->internal_units;
+  const struct hydro_props *hydro_props = e->hydro_properties;
   const double time_base = e->time_base;
   const integertime_t ti_current = e->ti_current;
   struct part *restrict parts = c->parts;
@@ -459,19 +460,25 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
 
       if (part_is_active(p, e)) {
 
-        double dt_cool;
+        double dt_cool, dt_therm;
+        ;
         if (with_cosmology) {
           const integertime_t ti_step = get_integer_timestep(p->time_bin);
           const integertime_t ti_begin =
               get_integer_time_begin(ti_current + 1, p->time_bin);
           dt_cool =
               cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
+          dt_therm = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
+                                                     ti_begin + ti_step);
+
         } else {
           dt_cool = get_timestep(p->time_bin, time_base);
+          dt_therm = get_timestep(p->time_bin, time_base);
         }
 
         /* Let's cool ! */
-        cooling_cool_part(constants, us, cosmo, cooling_func, p, xp, dt_cool);
+        cooling_cool_part(constants, us, cosmo, hydro_props, cooling_func, p,
+                          xp, dt_cool, dt_therm);
       }
     }
   }

src/timestep.h

@@ -126,8 +126,9 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
   /* Compute the next timestep (cooling condition) */
   float new_dt_cooling = FLT_MAX;
   if (e->policy & engine_policy_cooling)
-    new_dt_cooling = cooling_timestep(e->cooling_func, e->physical_constants,
-                                      e->cosmology, e->internal_units, p);
+    new_dt_cooling =
+        cooling_timestep(e->cooling_func, e->physical_constants, e->cosmology,
+                         e->internal_units, e->hydro_properties, p);
 
   /* Compute the next timestep (gravity condition) */
   float new_dt_grav = FLT_MAX, new_dt_self_grav = FLT_MAX,