diff --git a/doc/RTD/source/SubgridModels/GEAR/index.rst b/doc/RTD/source/SubgridModels/GEAR/index.rst index c1ac89f70be3bc2a75a99e2fb02ef515d7d67d5d..a3ada8807f5b85bd86783fd4034040d0ae98e6cd 100644 --- a/doc/RTD/source/SubgridModels/GEAR/index.rst +++ b/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``. diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml index 6c07371a47d7b0b15113e80b602ba5e197661bec..0947c0da43928544de0f2a5f1f973748d8dec09e 100644 --- a/examples/parameter_example.yml +++ b/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 diff --git a/src/cooling/grackle/cooling.c b/src/cooling/grackle/cooling.c index 215ec4cd9db8570dc1603a9f5f46ca8bb3afc473..317490ef2f32b7f600f210c8146e6b9feaa1a6db 100644 --- a/src/cooling/grackle/cooling.c +++ b/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, diff --git a/src/cooling/grackle/cooling_io.h b/src/cooling/grackle/cooling_io.h index 46ecaff67bd9b31a715de92679ef106ca60f9254..1ecdad94589e6e3c97946f23fcd6b94a39329e14 100644 --- a/src/cooling/grackle/cooling_io.h +++ b/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"); diff --git a/src/cooling/grackle/cooling_properties.h b/src/cooling/grackle/cooling_properties.h index 93b602cffdb4b6130bb773ffa66b7c0c0a6dd394..f08083b15b3a14f6b49dfc6a6df15dd6084518d7 100644 --- a/src/cooling/grackle/cooling_properties.h +++ b/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;