Merge branch 'GrackleCooling4' into 'master'

Small changes in the Grackle cooling

See merge request !1811

doc/RTD/source/SubgridModels/GEAR/index.rst

@@ -91,18 +91,30 @@ The self shielding method is defined by ``GrackleCooling:self_shielding_method``
 .. code:: YAML
 
   GrackleCooling:
-    cloudy_table: CloudyData_UVB=HM2012.h5       # Name of the Cloudy Table (available on the grackle bitbucket repository)
-    with_UV_background: 1                        # Enable or not the UV background
-    redshift: 0                                  # Redshift to use (-1 means time based redshift)
-    with_metal_cooling: 1                        # Enable or not the metal cooling
-    provide_volumetric_heating_rates: 0          # (optional) User provide volumetric heating rates
-    provide_specific_heating_rates: 0            # (optional) User provide specific heating rates
-    max_steps: 10000                             # (optional) Max number of step when computing the initial composition
-    convergence_limit: 1e-2                      # (optional) Convergence threshold (relative) for initial composition
-    thermal_time_myr: 5                          # (optional) Time (in Myr) for adiabatic cooling after a feedback event.
-    self_shielding_method: -1                    # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR)
-    self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding
-
+  cloudy_table: CloudyData_UVB=HM2012.h5       # Name of the Cloudy Table (available on the grackle bitbucket repository)
+  with_UV_background: 1                        # Enable or not the UV background
+  redshift: 0                                  # Redshift to use (-1 means time based redshift)
+  with_metal_cooling: 1                        # Enable or not the metal cooling
+  provide_volumetric_heating_rates: 0          # (optional) User provide volumetric heating rates
+  provide_specific_heating_rates: 0            # (optional) User provide specific heating rates
+  max_steps: 10000                             # (optional) Max number of step when computing the initial composition
+  convergence_limit: 1e-2                      # (optional) Convergence threshold (relative) for initial composition
+  thermal_time_myr: 5                          # (optional) Time (in Myr) for adiabatic cooling after a feedback event.
+  self_shielding_method: -1                    # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR)
+  self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding
+  HydrogenFractionByMass : 1.                  # Hydrogen fraction by mass (default is 0.76)
+  use_radiative_transfer : 1                   # Arrays of ionization and heating rates are provided
+  RT_heating_rate_cgs    : 0                   # heating         rate in units of / nHI_cgs 
+  RT_HI_ionization_rate_cgs  : 0               # HI ionization   rate in cgs [1/s]
+  RT_HeI_ionization_rate_cgs : 0               # HeI ionization  rate in cgs [1/s]
+  RT_HeII_ionization_rate_cgs: 0               # HeII ionization rate in cgs [1/s]
+  RT_H2_dissociation_rate_cgs: 0               # H2 dissociation rate in cgs [1/s]
+  volumetric_heating_rates_cgs: 0              # Volumetric heating rate in cgs  [erg/s/cm3]
+  specific_heating_rates_cgs: 0                # Specific heating rate in cgs    [erg/s/g]
+  H2_three_body_rate : 1                       # Specific the H2 formation three body rate (0->5,see Grackle documentation)
+  H2_cie_cooling : 0                           # Enable/disable H2 collision-induced emission cooling from Ripamonti & Abel (2004)
+  cmb_temperature_floor : 1                    # Enable/disable an effective CMB temperature floor
+  
 .. note::
    A simple example running SWIFT with Grackle can be find in ``examples/Cooling/CoolingBox``. A more advanced example combining heating and cooling (with heating and ionization sources) is given in ``examples/Cooling/CoolingHeatingBox``.
 

examples/parameter_example.yml

@@ -545,8 +545,9 @@ GrackleCooling:
   RT_H2_dissociation_rate_cgs: 0               # H2 dissociation rate in cgs [1/s]
   volumetric_heating_rates_cgs: 0              # Volumetric heating rate in cgs  [erg/s/cm3]
   specific_heating_rates_cgs: 0                # Specific heating rate in cgs    [erg/s/g]
-  
-  
+  H2_three_body_rate : 1                       # Specific the H2 formation three body rate (0->5,see Grackle documentation)
+  H2_cie_cooling : 0                           # Enable/disable H2 collision-induced emission cooling from Ripamonti & Abel (2004)
+  cmb_temperature_floor : 1                    # Enable/disable an effective CMB temperature floor
   
 
 

src/cooling/grackle/cooling.c

@@ -179,45 +179,7 @@ void cooling_compute_equilibrium(const struct phys_const* phys_const,
                                  const struct cooling_function_data* cooling,
                                  const struct part* p, struct xpart* xp) {
 
-  /* TODO: this can fail spectacularly and needs to be replaced. */
-
-  /* get temporary data */
-  struct part p_tmp = *p;
-  struct cooling_function_data cooling_tmp = *cooling;
-  cooling_tmp.chemistry_data.with_radiative_cooling = 0;
-  /* need density for computation, therefore quick estimate */
-  p_tmp.rho = 0.2387 * p_tmp.mass / pow(p_tmp.h, 3);
-
-  /* compute time step */
-  const double alpha = 0.01;
-  double dt = fabs(cooling_time(phys_const, us, hydro_properties, cosmo,
-                                &cooling_tmp, &p_tmp, xp));
-  cooling_new_energy(phys_const, us, cosmo, hydro_properties, &cooling_tmp,
-                     &p_tmp, xp, dt, dt);
-  dt = alpha * fabs(cooling_time(phys_const, us, hydro_properties, cosmo,
-                                 &cooling_tmp, &p_tmp, xp));
-
-  /* init simple variables */
-  int step = 0;
-  const int max_step = cooling_tmp.max_step;
-  const float conv_limit = cooling_tmp.convergence_limit;
-  struct xpart old;
-
-  do {
-    /* update variables */
-    step += 1;
-    old = *xp;
-
-    /* update chemistry */
-    cooling_new_energy(phys_const, us, cosmo, hydro_properties, &cooling_tmp,
-                       &p_tmp, xp, dt, dt);
-  } while (step < max_step && !cooling_converged(xp, &old, conv_limit));
-
-  if (step == max_step)
-    error(
-        "A particle element fraction failed to converge."
-        "You can change 'GrackleCooling:MaxSteps' or "
-        "'GrackleCooling:ConvergenceLimit' to avoid this problem");
+  return;
 }
 
 /**
@@ -249,15 +211,13 @@ void cooling_first_init_part(const struct phys_const* phys_const,
    * a better determination will be done in cooling_post_init_part */
 
   /* primordial chemistry >= 1 */
+  /* assume neutral gas */
   xp->cooling_data.HI_frac = grackle_data->HydrogenFractionByMass;
   xp->cooling_data.HII_frac = zero;
-  xp->cooling_data.HeI_frac = zero;
+  xp->cooling_data.HeI_frac = 1. - grackle_data->HydrogenFractionByMass;
   xp->cooling_data.HeII_frac = zero;
-  xp->cooling_data.HeIII_frac = 1. - grackle_data->HydrogenFractionByMass;
-  xp->cooling_data.e_frac = xp->cooling_data.HII_frac +
-                            0.25 * xp->cooling_data.HeII_frac +
-                            0.5 * xp->cooling_data.HeIII_frac;
-
+  xp->cooling_data.HeIII_frac = zero;
+  xp->cooling_data.e_frac = zero;
 #endif  // MODE >= 1
 
 #if COOLING_GRACKLE_MODE >= 2
@@ -301,9 +261,9 @@ void cooling_post_init_part(const struct phys_const* phys_const,
   // message("rho = %g energy = %g",rho,energy);
 
 #if COOLING_GRACKLE_MODE > 0
-  /* TODO: this can fail spectacularly and needs to be replaced. */
-  // cooling_compute_equilibrium(phys_const, us, hydro_props, cosmo, cooling,
-  // p,xp);
+  /* The function below currently does nothing. Will have to be updated. */
+  cooling_compute_equilibrium(phys_const, us, hydro_props, cosmo, cooling, p,
+                              xp);
 #endif
 }
 
@@ -366,6 +326,12 @@ void cooling_print_backend(const struct cooling_function_data* cooling) {
           cooling->chemistry_data.with_radiative_cooling);
   message("grackle_chemistry_data.primordial_chemistry = %d",
           cooling->chemistry_data.primordial_chemistry);
+  message("grackle_chemistry_data.three_body_rate = %d",
+          cooling->chemistry_data.three_body_rate);
+  message("grackle_chemistry_data.cmb_temperature_floor = %d",
+          cooling->chemistry_data.cmb_temperature_floor);
+  message("grackle_chemistry_data.cie_cooling = %d",
+          cooling->chemistry_data.cie_cooling);
   message("grackle_chemistry_data.dust_chemistry = %d",
           cooling->chemistry_data.dust_chemistry);
   message("grackle_chemistry_data.metal_cooling = %d",
@@ -1149,9 +1115,19 @@ void cooling_init_grackle(struct cooling_function_data* cooling) {
   chemistry->primordial_chemistry = cooling->primordial_chemistry;
   chemistry->metal_cooling = cooling->with_metal_cooling;
   chemistry->UVbackground = cooling->with_uv_background;
+  chemistry->three_body_rate = cooling->H2_three_body_rate;
+  chemistry->cmb_temperature_floor = cooling->cmb_temperature_floor;
+  chemistry->cie_cooling = cooling->H2_cie_cooling;
   chemistry->grackle_data_file = cooling->cloudy_table;
 
   /* radiative transfer */
+#if COOLING_GRACKLE_MODE == 0
+  if (cooling->use_radiative_transfer)
+    error(
+        "The parameter use_radiative_transfer cannot be set to 1 in Grackle "
+        "mode 0 !");
+#endif
+
   chemistry->use_radiative_transfer = cooling->use_radiative_transfer;
 
   if (cooling->volumetric_heating_rates > 0)
@@ -1173,7 +1149,10 @@ void cooling_init_grackle(struct cooling_function_data* cooling) {
         "heating rates, not both");
 
   /* self shielding */
-  chemistry->self_shielding_method = cooling->self_shielding_method;
+  if (cooling->self_shielding_method <= 0)
+    chemistry->self_shielding_method = 0;
+  else
+    chemistry->self_shielding_method = cooling->self_shielding_method;
 
   if (local_initialize_chemistry_data(&cooling->chemistry_data,
                                       &cooling->chemistry_rates,

src/cooling/grackle/cooling_io.h

@@ -156,6 +156,15 @@ __attribute__((always_inline)) INLINE static void cooling_read_parameters(
     error("Cannot run primordial chemistry %i when compiled with %i",
           cooling->primordial_chemistry, COOLING_GRACKLE_MODE);
 
+  cooling->H2_three_body_rate = parser_get_opt_param_int(
+      parameter_file, "GrackleCooling:H2_three_body_rate", 0);
+
+  cooling->H2_cie_cooling = parser_get_opt_param_int(
+      parameter_file, "GrackleCooling:H2_cie_cooling", 0);
+
+  cooling->cmb_temperature_floor = parser_get_opt_param_int(
+      parameter_file, "GrackleCooling:cmb_temperature_floor", 1);
+
   cooling->with_uv_background =
       parser_get_param_int(parameter_file, "GrackleCooling:with_UV_background");
 

src/cooling/grackle/cooling_properties.h

@@ -44,6 +44,16 @@ struct cooling_function_data {
   /*! Chemistry network */
   int primordial_chemistry;
 
+  /*! Set the three-body reaction rate (see grackle documentation) */
+  int H2_three_body_rate;
+
+  /*! Enable/disable H2 collision-induced emission cooling from Ripamonti & Abel
+   * (2004) */
+  int H2_cie_cooling;
+
+  /*! Enable/disable CMB temperature floor */
+  int cmb_temperature_floor;
+
   /*! Redshift to use for the UV backgroud (-1 to use cosmological one) */
   double redshift;