/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see .
*
******************************************************************************/
/**
* @file src/cooling/CHIMES/cooling.c
* @brief CHIMES cooling functions
*/
/* Config parameters. */
#include
/* Some standard headers. */
#include
#include
#include
#include
/* Local includes. */
#include "adiabatic_index.h"
#include "chemistry.h"
#include "chimes/src/chimes_interpol.h"
#include "colibre_subgrid.h"
#include "cooling_properties.h"
#include "dust.h"
#include "entropy_floor.h"
#include "error.h"
#include "exp10.h"
#include "hydro.h"
#include "io_properties.h"
#include "parser.h"
#include "part.h"
#include "physical_constants.h"
#include "space.h"
#include "star_formation.h"
#include "units.h"
const char *CHIMES_tables_version_string = "14Oct2022";
/* Pre-declarations to simplify the code */
double chimes_mu(const struct cooling_function_data *cooling,
const struct hydro_props *hydro_props, const struct part *p,
const struct xpart *xp);
void chimes_update_element_abundances(
const struct phys_const *phys_const, const struct unit_system *us,
const struct cosmology *cosmo, const struct cooling_function_data *cooling,
const struct hydro_props *hydro_props, const struct part *p,
const struct xpart *xp, struct gasVariables *ChimesGasVars,
const double u_0, const int mode);
void chimes_update_gas_vars(
const double u_cgs, const struct phys_const *phys_const,
const struct unit_system *us, const struct cosmology *cosmo,
const struct hydro_props *hydro_properties,
const struct entropy_floor_properties *floor_props,
const struct cooling_function_data *cooling, const struct dustevo_props *dp,
const struct part *p, const struct xpart *xp,
struct gasVariables *ChimesGasVars, const float dt_cgs);
float cooling_get_temperature_on_entropy_floor(
const struct phys_const *phys_const, const struct hydro_props *hydro_props,
const struct unit_system *us, const struct cosmology *cosmo,
const struct cooling_function_data *cooling, const struct part *p,
const struct xpart *xp);
double calculate_colibre_N_ref(const struct phys_const *phys_const,
const struct unit_system *us,
const struct cosmology *cosmo,
const struct cooling_function_data *cooling,
const struct hydro_props *hydro_properties,
const struct part *p, const struct xpart *xp,
const double temperature);
void cooling_set_HIIregion_chimes_abundances(
struct gasVariables *ChimesGasVars,
const struct cooling_function_data *cooling);
void cooling_set_FB_particle_chimes_abundances(
struct gasVariables *ChimesGasVars,
const struct cooling_function_data *cooling);
/**
* @brief Returns whether the gas particle is in an HII region or not.
*/
int is_HII_region(const struct xpart *xp) {
return xp->tracers_data.HIIregion_timer_gas > 0.;
}
/**
* @brief Computes the mean molecular weight for an HII region at a given H
* abundance.
*
* @param cooling The properties of the cooling scheme.
* @param XH The Hydrogen abundance of the gas.
*/
double get_mu_HII_region(const struct cooling_function_data *cooling,
const double XH) {
return 4.0 / ((1.0 + cooling->HIIregion_fion) * (1.0 + (3.0 * XH)));
}
/**
* @brief Compute the temperature from mu and the internal energy in CGS.
*/
double get_T_from_u_cgs(const double u_cgs, const double mu,
const struct cooling_function_data *cooling) {
return u_cgs * hydro_gamma_minus_one * mu *
cooling->proton_mass_over_boltzmann_k_cgs;
}
/**
* @brief Compute the temperature from mu and the internal energy in internal
* units.
*/
double get_T_from_u_internal(const double u, const double mu,
const struct phys_const *phys_const) {
return u * hydro_gamma_minus_one * mu * phys_const->const_proton_mass /
phys_const->const_boltzmann_k;
}
/**
* @brief Compute the the internal energy in CGS from mu and the temperature.
*/
double get_u_cgs_from_T(const double T, const double mu,
const struct cooling_function_data *cooling) {
return T * hydro_one_over_gamma_minus_one /
(cooling->proton_mass_over_boltzmann_k_cgs * mu);
}
/**
* @brief Compute the the internal energy in internal units from mu and the
* temperature.
*/
double get_u_internal_from_T(const double T, const double mu,
const struct phys_const *phys_const) {
return T * hydro_one_over_gamma_minus_one * phys_const->const_boltzmann_k /
(phys_const->const_proton_mass * mu);
}
/**
* @brief Return the hydrogen mass fraction of a particle.
*
* Deals with the case where we don't have a chemistry model.
*/
double get_hydrogen_mass_fraction(const struct part *p,
const struct hydro_props *hydro_props) {
#if defined(CHEMISTRY_COLIBRE) || defined(CHEMISTRY_EAGLE)
const float *const metal_fraction =
chemistry_get_metal_mass_fraction_for_cooling(p);
return metal_fraction[chemistry_element_H];
#else
/* Without COLIBRE or EAGLE chemistry, the metal abundances are unavailable.
* --> Set to primordial abundances. */
return hydro_props->hydrogen_mass_fraction;
#endif /* CHEMISTRY_COLIBRE || CHEMISTRY_EAGLE */
}
/**
* @brief Get the dust ratio for a given particle
* and for the different dust modes
*/
double get_dust_ratio(const struct part *p,
const struct cooling_function_data *cooling,
const struct dustevo_props *dp, const double Z) {
#if defined(DUST_NONE)
return Z;
#else
if (dp->pair_to_cooling) {
float dfrac = 0.;
for (int grain = 0; grain < dust_grain_species_count; grain++) {
double grainfrac = p->dust_data.grain_mass_fraction[grain];
#if defined(DUST_T20MS)
if (grain > (dust_grain_species_count / 2)) {
/* the second half of the grain array has a smaller radius.
* Grain processes scale with cross section, so we apply
* a boost factor based on the size ratios. */
grainfrac *= dust_get_large_to_small_ratio();
}
#endif
dfrac += grainfrac;
}
return dfrac / dust_get_local_dust_to_gas_ratio();
} else {
return Z;
}
#endif
}
double get_dust_boost_factor(const double nH_cgs,
const struct cooling_function_data *cooling) {
if (nH_cgs <= cooling->dust_boost_nH_min_cgs) {
/* Below the min density for boost --> no boost */
return 1.0f;
} else if (nH_cgs >= cooling->dust_boost_nH_max_cgs) {
/* Above the max density for boost --> max boost */
return cooling->max_dust_boost_factor;
} else {
/* In the transition region --> Compute the boost factor */
const float log_nH = log10f(nH_cgs);
const float log_nH_min = log10f(cooling->dust_boost_nH_min_cgs);
const float log_nH_max = log10f(cooling->dust_boost_nH_max_cgs);
const float log_max_boost = log10f(cooling->max_dust_boost_factor);
const float delta_log_nH = chimes_max(log_nH_max - log_nH_min, FLT_MIN);
/* Linear ramp (in log) between min and max */
const float log_boost =
((log_nH - log_nH_min) / delta_log_nH) * log_max_boost;
return exp10f(log_boost);
}
}
double get_CR_rate(const double N_ref_cgs, const double nH_cgs,
const struct cosmology *cosmo,
const struct phys_const *phys_const,
const struct cooling_function_data *cooling) {
/* Set the cosmic ray rate */
if (cooling->UV_field_choice != UV_field_COLIBRE) {
/* Constant cosmic ray rate */
return cooling->cosmic_ray_rate;
} else {
/* The cosmic ray rate scales with the reference column density N_ref.
* Based on Ploeckinger & Schaye (2020) */
const double cr_coldensratio =
N_ref_cgs / cooling->cr_transition_column_density_cgs;
double cosmic_ray_rate;
if (cr_coldensratio < 1.) {
/* N_ref < N_trans: power law index 1 */
cosmic_ray_rate = pow(cr_coldensratio, cooling->cr_powerlaw_slope_1);
} else {
/* N_ref >=N_trans: power law index 2 */
cosmic_ray_rate = pow(cr_coldensratio, cooling->cr_powerlaw_slope_2);
}
cosmic_ray_rate *= cooling->cr_rate_at_transition_cgs;
const double log_cosmic_ray_rate_cgs = log10(cosmic_ray_rate);
/* Now, apply some saturation */
const double log_cr_rate_min =
log10(cooling->cr_rate_at_transition_cgs *
pow(cooling->nref_column_density_min_cgs /
cooling->cr_transition_column_density_cgs,
cooling->cr_powerlaw_slope_1));
/* low density cutoff between cr_nmax and cr_nmin
* with k regulating the steepness */
const double cr_nmax_over_dens =
cosmo->critical_density * cooling->cr_cutoff_over_density *
cooling->nref_XH_const / phys_const->const_proton_mass;
const double cr_nmax_over_dens_cgs =
cr_nmax_over_dens * cooling->number_density_to_cgs;
/* Max density */
const double cr_nmax =
fmax(cr_nmax_over_dens_cgs, cooling->cr_cutoff_phys_density_cgs);
/* Min density: (1. / 10^(cutoff_width_dex)) * cr_max */
const double cr_nmin = exp10(-cooling->cr_cutoff_width_dex) * cr_nmax;
const double numerator = (log_cosmic_ray_rate_cgs - log_cr_rate_min);
const double argument = sqrt(cr_nmin / cr_nmax) * cr_nmax / nH_cgs;
const double denominator =
1. + pow(argument, -cooling->cr_cutoff_steepness_k);
const double log_cosmic_ray_rate =
log_cosmic_ray_rate_cgs - numerator / denominator;
return exp10(log_cosmic_ray_rate);
}
}
/**
* @brief Find the index of the current redshift along the redshift dimension
* of the cooling tables.
*
* Since the redshift table is not evenly spaced, compare z with each
* table value in decreasing order starting with the previous redshift index
*
* The returned difference is expressed in units of the table separation. This
* means dx = (x - table[i]) / (table[i+1] - table[i]). It is always between
* 0 and 1.
*
* @param z Redshift we are searching for.
* @param z_index (return) Index of the redshift in the table.
* @param dz (return) Difference in redshift between z and table[z_index].
* @param cooling #cooling_function_data structure containing redshift table.
*/
__attribute__((always_inline)) INLINE void get_redshift_index(
const float z, int *z_index, float *dz,
struct colibre_cooling_tables *table) {
/* If z is larger than the maximum table redshift, use the last value */
if (z >= table->Redshifts[table->number_of_redshifts - 1]) {
*z_index = table->number_of_redshifts - 2;
*dz = 1.0;
}
/* At the end, just use the last value */
else if (z <= table->Redshifts[0]) {
*z_index = 0;
*dz = 0.0;
}
/* Normal case: search... */
else {
/* start at the previous index and search */
for (int i = table->previous_z_index; i >= 0; i--) {
if (z > table->Redshifts[i]) {
*z_index = i;
table->previous_z_index = i;
*dz = (z - table->Redshifts[i]) /
(table->Redshifts[i + 1] - table->Redshifts[i]);
break;
}
}
}
}
/**
* @brief Checks CHIMES eqm tables headers.
*
* If we are using the redshift-dependent
* CHIMES eqm abundances with the COLIBRE
* radiation field and shielding models,
* we need to check that the parameters
* used in the tables match the
* hard-coded values.
*
* @param fname_chimes filename of the eqm table.
* @param cooling #cooling_function_data struct.
*/
void chimes_check_colibre_eqm_tables_header(
const char *fname_chimes, const struct cooling_function_data *cooling) {
const hid_t file_id = H5Fopen(fname_chimes, H5F_ACC_RDONLY, H5P_DEFAULT);
if (file_id < 0)
error("Unable to open equilibrium abundance data file '%s'!", fname_chimes);
/* Check table version */
const hid_t attr = H5Aopen_by_name(file_id, "/Header", "version_date",
H5P_DEFAULT, H5P_DEFAULT);
if (attr < 0)
error(
"Error reading the version string. You are likely using old tables "
"that do not match this code version.");
const hid_t type = H5Aget_type(attr);
char *version_string = NULL;
const hid_t status = H5Aread(attr, type, &version_string);
if (status < 0)
error("Error reading version_date in CHIMES eqm table header!");
if (strcmp(version_string, CHIMES_tables_version_string) != 0)
error(
"The table and code version do not match, please use table version "
"'%s' (you are using version '%s').",
CHIMES_tables_version_string, version_string);
H5free_memory(version_string);
H5Tclose(type);
H5Aclose(attr);
H5Fclose(file_id);
}
/**
* @brief Initialises properties stored in the cooling_function_data struct
*
* @param parameter_file The parsed parameter file
* @param us Internal system of units data structure
* @param phys_const #phys_const data structure
* @param hydro_props The properties of the hydro scheme.
* @param cooling #cooling_function_data struct to initialize
*/
void cooling_init_backend(struct swift_params *parameter_file,
const struct unit_system *us,
const struct phys_const *phys_const,
const struct hydro_props *hydro_props,
struct cooling_function_data *cooling,
struct dustevo_props *dp) {
char chimes_data_dir[256];
char string_buffer[196];
/* The following CHIMES variables control the options involved with coupling
* the thermo-chemistry solver in CHIMES to a Radiative Transfer (RT) solver.
* However, we have not yet implemented the coupling between CHIMES and the
* RT modules in Swift. For now, we therefore simply switch off the RT
* coupling options in CHIMES. */
cooling->ChimesGlobalVars.rt_update_flux = 0;
cooling->ChimesGlobalVars.rt_use_on_the_spot_approx = 0;
cooling->ChimesGlobalVars.rt_density_absoluteTolerance = 1.0e-10f;
cooling->ChimesGlobalVars.rt_flux_absoluteTolerance = 1.0e-10f;
/* read the parameters */
/* Set paths to CHIMES data files. */
parser_get_param_string(parameter_file, "CHIMESCooling:data_path",
chimes_data_dir);
sprintf(cooling->ChimesGlobalVars.MainDataTablePath,
"%s/chimes_main_data.hdf5", chimes_data_dir);
parser_get_param_string(parameter_file, "CHIMESCooling:EqmAbundanceTable",
string_buffer);
if ((strcmp("None", string_buffer) == 0) ||
(strcmp("none", string_buffer) == 0))
sprintf(cooling->ChimesGlobalVars.EqAbundanceTablePath, "%s",
string_buffer);
else
sprintf(cooling->ChimesGlobalVars.EqAbundanceTablePath,
"%s/EqAbundancesTables/%s", chimes_data_dir, string_buffer);
parser_get_param_string(parameter_file, "CHIMESCooling:UV_field",
string_buffer);
if (strcmp(string_buffer, "none") == 0) {
cooling->UV_field_choice = UV_field_no_UV;
} else if (strcmp(string_buffer, "user-defined") == 0) {
cooling->UV_field_choice = UV_field_user_defined;
} else if (strcmp(string_buffer, "COLIBRE") == 0) {
cooling->UV_field_choice = UV_field_COLIBRE;
} else {
error(
"Invalid choice of 'CHIMESCooling:UV_field'. Must be 'none', "
"'user-defined', or 'COLIBRE'");
}
parser_get_param_string(parameter_file, "CHIMESCooling:UVB_z_dependence",
string_buffer);
if (strcmp(string_buffer, "constant") == 0) {
cooling->UVB_dependence_choice = UVB_constant;
} else if (strcmp(string_buffer, "cross sections") == 0) {
cooling->UVB_dependence_choice = UVB_zdep_cross_secs;
} else if (strcmp(string_buffer, "COLIBRE") == 0) {
cooling->UVB_dependence_choice = UVB_zdep_cross_secs_and_eqm;
} else {
error(
"Invalid choice of 'CHIMESCooling:UVB_z_dependence'. Must be "
"'constant', 'cross sections', or 'COLIBRE'");
}
parser_get_param_string(parameter_file, "CHIMESCooling:Shielding_model",
string_buffer);
if (strcmp(string_buffer, "none") == 0) {
cooling->Shielding_choice = Shielding_none;
} else if (strcmp(string_buffer, "Jeans length") == 0) {
cooling->Shielding_choice = Shielding_Jeans_length;
} else if (strcmp(string_buffer, "COLIBRE") == 0) {
cooling->Shielding_choice = Shielding_COLIBRE;
} else {
error(
"Invalid choice of 'CHIMESCooling:Shielding_model'. Must be 'none', "
"'Jeans length', or 'COLIBRE'");
}
cooling->rad_field_norm_factor = 0.0;
if (cooling->UVB_dependence_choice == UVB_constant) {
cooling->ChimesGlobalVars.redshift_dependent_UVB_index = -1;
cooling->ChimesGlobalVars.use_redshift_dependent_eqm_tables = 0;
} else {
if (cooling->UV_field_choice == UV_field_no_UV)
error(
"CHIMES ERROR: redshift-dependent UVB has been switched on, but "
"UV_field_choice is 'no_UV'. These options are incompatible.");
cooling->ChimesGlobalVars.redshift_dependent_UVB_index = 0;
if (cooling->UVB_dependence_choice == UVB_zdep_cross_secs)
cooling->ChimesGlobalVars.use_redshift_dependent_eqm_tables = 0;
else if (cooling->UVB_dependence_choice == UVB_zdep_cross_secs_and_eqm)
cooling->ChimesGlobalVars.use_redshift_dependent_eqm_tables = 1;
else
error("Invalid choice of 'CHIMESCooling:UVB_z_dependence");
}
if (cooling->UV_field_choice == UV_field_no_UV) {
cooling->ChimesGlobalVars.N_spectra = 0;
} else if (cooling->UV_field_choice == UV_field_user_defined) {
cooling->rad_field_norm_factor = (ChimesFloat)parser_get_opt_param_double(
parameter_file, "CHIMESCooling:rad_field_norm_factor", 1.0);
cooling->ChimesGlobalVars.N_spectra = 1;
parser_get_param_string(parameter_file, "CHIMESCooling:PhotoIonTable",
string_buffer);
sprintf(cooling->ChimesGlobalVars.PhotoIonTablePath[0], "%s/%s",
chimes_data_dir, string_buffer);
} else if (cooling->UV_field_choice == UV_field_COLIBRE) {
cooling->ChimesGlobalVars.N_spectra = 2;
parser_get_param_string(parameter_file, "CHIMESCooling:PhotoIonTable_UVB",
string_buffer);
sprintf(cooling->ChimesGlobalVars.PhotoIonTablePath[0], "%s/%s",
chimes_data_dir, string_buffer);
parser_get_param_string(parameter_file, "CHIMESCooling:PhotoIonTable_ISRF",
string_buffer);
sprintf(cooling->ChimesGlobalVars.PhotoIonTablePath[1], "%s/%s",
chimes_data_dir, string_buffer);
} else {
error("Invalid 'UV_field_choice'");
}
if (cooling->Shielding_choice == Shielding_none)
cooling->ChimesGlobalVars.cellSelfShieldingOn = 0;
else if (cooling->Shielding_choice == Shielding_Jeans_length)
cooling->ChimesGlobalVars.cellSelfShieldingOn = 1;
else if (cooling->Shielding_choice == Shielding_COLIBRE)
cooling->ChimesGlobalVars.cellSelfShieldingOn = 1;
if (cooling->Shielding_choice == Shielding_Jeans_length) {
/* Maximum shielding length (in kpc).
*If negative, do not impose a maximum. */
cooling->max_shielding_length_kpc = parser_get_opt_param_double(
parameter_file, "CHIMESCooling:max_shielding_length_kpc", -1.0);
/* Convert the maximum shielding length to cm. */
cooling->max_shielding_length_cm =
cooling->max_shielding_length_kpc * 1.0e3 *
units_cgs_conversion_factor(us, UNIT_CONV_LENGTH) *
phys_const->const_parsec;
}
if (cooling->Shielding_choice != Shielding_none) {
/* This parameter is needed if shielding is switched on */
/* Factor to re-scale shielding length. */
cooling->shielding_length_factor = parser_get_opt_param_double(
parameter_file, "CHIMESCooling:shielding_length_factor", 1.0);
}
if (cooling->Shielding_choice == Shielding_COLIBRE) {
/* The maximum shielding length can depend on redshift */
/* in the COLIBRE shielding model */
cooling->lshmax_lowz_kpc = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:LSHmax_lowz_kpc", 50.0);
cooling->lshmax_highz_kpc = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:LSHmax_highz_kpc", 10.0);
cooling->lshmax_transition_redshift = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:LSHmax_transition_redshift", 7.0);
cooling->lshmax_transition_width_dz = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:LSHmax_transition_width_dz", 0.2);
cooling->lshmax_transition_steepness_k = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:LSHmax_transition__steepness_k", 100.);
}
if ((cooling->Shielding_choice == Shielding_COLIBRE) ||
(cooling->UV_field_choice == UV_field_COLIBRE)) {
/* Read in the parameters needed for calculating Nref;
* The reference column density is used for the shielding length
* if Shielding_flag == COLIBRE
* and for the ISRF if UV_field_flag == COLIBRE */
/* Constant hydrogen mass fraction to calculate Nref */
cooling->nref_XH_const = parser_get_opt_param_double(
parameter_file, "CHIMESCooling:NREF_XH_const", 0.75);
/* Constant mean particle mass to calculate Nref */
cooling->nref_mu_const = parser_get_opt_param_double(
parameter_file, "CHIMESCooling:NREF_mu_const", 1.24);
/* Minimum column density for Nref */
cooling->nref_column_density_min_cgs = parser_get_opt_param_double(
parameter_file, "CHIMESCooling:NREF_column_density_min_cgs", 3.086e15);
/* Maximum column density for Nref */
cooling->nref_column_density_max_cgs = parser_get_opt_param_double(
parameter_file, "CHIMESCooling:NREF_column_density_max_cgs", 1.0e24);
/* Minimum temperature for transition to nref_column_density_min_cgs */
cooling->nref_trans_temperature_min_K = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:NREF_trans_temperature_min_K", 1.e4);
/* Maximum temperature for transition to nref_column_density_min_cgs */
cooling->nref_trans_temperature_max_K = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:NREF_trans_temperature_max_K", 1.e5);
/* Steepness of the transition to nref_column_density_min_cgs */
cooling->nref_trans_steepness_k = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:NREF_trans_steepness_k", 2.);
/* Maximum length of reference column density */
cooling->nref_max_length_kpc = parser_get_opt_param_double(
parameter_file, "CHIMESCooling:NREF_column_length_max_kpc", 50.);
/* Convert the maximum length of reference column density to cm */
cooling->nref_max_length_cgs =
cooling->nref_max_length_kpc * 1.0e3 *
units_cgs_conversion_factor(us, UNIT_CONV_LENGTH) *
phys_const->const_parsec;
}
if (cooling->UV_field_choice == UV_field_COLIBRE) {
/* cosmic ray rate in s-1 at the transition column density */
cooling->cr_rate_at_transition_cgs = parser_get_opt_param_double(
parameter_file, "CHIMESCooling:CR_rate_at_transition_cgs", 2.e-16);
/* cosmic ray transition column density from power law 1 to 2 in cm-2 */
cooling->cr_transition_column_density_cgs = parser_get_opt_param_double(
parameter_file, "CHIMESCooling:CR_transition_column_density_cgs",
1.e21);
/* cosmic rays: slope of power law 1 */
cooling->cr_powerlaw_slope_1 = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:CR_powerlaw_slope_1", 1.4);
/* cosmic rays: slope of power law 2 */
cooling->cr_powerlaw_slope_2 = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:CR_powerlaw_slope_2", 0.0);
/* cosmic rays: overdensity (rel to rho_crit) cutoff */
cooling->cr_cutoff_over_density = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:CR_cutoff_over_density", 100.);
/* cosmic rays: physical density cutoff in cm-3 */
cooling->cr_cutoff_phys_density_cgs = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:CR_cutoff_phys_density_cgs", 1.e-2);
/* cosmic rays: width of cutoff in dex (gas density) */
cooling->cr_cutoff_width_dex = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:CR_cutoff_width_dex", 1.);
/* cosmic rays: steepness of the cutoff */
cooling->cr_cutoff_steepness_k = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:CR_cutoff_steepness_k", 2.);
/* Ionising photon density relative to that in the cross sections file
* at the transition column density */
cooling->isrf_norm_at_transition = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:ISRF_norm_at_transition", 19.);
/* ISRF transition column density from power law 1 to 2 in cm-2 */
cooling->isrf_transition_column_density_cgs = parser_get_opt_param_double(
parameter_file, "CHIMESCooling:ISRF_transition_column_density_cgs",
1.e22);
/* ISRF: slope of power law 1 */
cooling->isrf_powerlaw_slope_1 = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:ISRF_powerlaw_slope_1", 1.4);
/* ISRF: slope of power law 2 */
cooling->isrf_powerlaw_slope_2 = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:ISRF_powerlaw_slope_2", 0.0);
/* ISRF: overdensity (rel to rho_crit) cutoff */
cooling->isrf_cutoff_over_density = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:ISRF_cutoff_over_density", 100.);
/* ISRF: physical cutoff in cm-3 */
cooling->isrf_cutoff_phys_density_cgs = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:ISRF_cutoff_phys_density_cgs", 1.e-2);
/* ISRF: width of cutoff in dex (gas density) */
cooling->isrf_cutoff_width_dex = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:ISRF_cutoff_width_dex", 1.);
/* ISRF: steepness of the cutoff */
cooling->isrf_cutoff_steepness_k = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:ISRF_cutoff_steepness_k", 2.);
} else {
/* Cosmic ray ionisation rate of HI. */
cooling->cosmic_ray_rate = (ChimesFloat)parser_get_param_double(
parameter_file, "CHIMESCooling:cosmic_ray_rate_cgs");
}
/* In CHIMES, the UVB is switched off above
* the redshift of H-reionisation. */
cooling->ChimesGlobalVars.reionisation_redshift =
(ChimesFloat)parser_get_param_float(parameter_file,
"CHIMESCooling:UVB_cutoff_z");
/* Parameters used for the COLIBRE ISRF and shielding length. These have
* just been hard-coded for now
* The values are taken from Ploeckinger & Schaye (2020). */
cooling->N_H0 = 3.63e20; /* cgs */
/* Parameter to control use of turbulent Jeans length.
* Only used with the Colibre UV field and shielding options. */
cooling->colibre_use_turbulent_jeans_length = parser_get_param_int(
parameter_file, "CHIMESCooling:colibre_use_turbulent_jeans_length");
/* Turbulent 1D velocity dispersion, in km/s. */
cooling->sigma_turb_km_p_s = (ChimesFloat)parser_get_param_float(
parameter_file, "CHIMESCooling:turbulent_velocity_dispersion_km_p_s");
parser_get_param_string(parameter_file, "CHIMESCooling:init_abundance_mode",
string_buffer);
if (strcmp(string_buffer, "single") == 0) {
cooling->init_abundance_choice = init_abundance_single;
} else if (strcmp(string_buffer, "read") == 0) {
cooling->init_abundance_choice = init_abundance_read;
} else if (strcmp(string_buffer, "compute") == 0) {
cooling->init_abundance_choice = init_abundance_compute;
} else {
error(
"Invalid choice of 'CHIMESCooling:init_abundance_mode'. Must be "
"'single', 'read', or 'compute'");
}
if (cooling->init_abundance_choice == init_abundance_single)
cooling->InitIonState =
parser_get_param_int(parameter_file, "CHIMESCooling:InitIonState");
else
cooling->InitIonState = 0;
/* CHIMES tolerance parameters */
cooling->ChimesGlobalVars.relativeTolerance =
(ChimesFloat)parser_get_param_double(parameter_file,
"CHIMESCooling:relativeTolerance");
cooling->ChimesGlobalVars.absoluteTolerance =
(ChimesFloat)parser_get_param_double(parameter_file,
"CHIMESCooling:absoluteTolerance");
cooling->ChimesGlobalVars.explicitTolerance =
(ChimesFloat)parser_get_param_double(parameter_file,
"CHIMESCooling:explicitTolerance");
cooling->ChimesGlobalVars.scale_metal_tolerances =
(ChimesFloat)parser_get_param_double(
parameter_file, "CHIMESCooling:scale_metal_tolerances");
/* Maximum temperature for the molecular network */
cooling->ChimesGlobalVars.Tmol_K = (ChimesFloat)parser_get_param_double(
parameter_file, "CHIMESCooling:Tmol_K");
/* For now, the CMB temperature in CHIMES is just set
* to the present day value. For cosmological runs, this
* will need to be updated for the current redshift. */
/* CMB temperature at redshift zero. */
cooling->T_CMB_0 = phys_const->const_T_CMB_0 *
units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
/* Initially set the CMB temperature in CHIMES to the redshift zero value.
* This will be updated to the current redshift in cooling_update(). */
cooling->ChimesGlobalVars.cmb_temperature = (ChimesFloat)cooling->T_CMB_0;
/* Equilibrium mode:
* 0 --> Evolve in non-equilibrium.
* 1 --> Evolve with equilibrium abundances.
*/
cooling->ChemistryEqmMode =
parser_get_param_int(parameter_file, "CHIMESCooling:ChemistryEqmMode");
/* Flag switch thermal evolution on/off:
* 0 --> Fixed temperature (but chemical abundances can still be evolved).
* 1 --> Enable thermal evolution.
*/
cooling->ThermEvolOn =
parser_get_param_int(parameter_file, "CHIMESCooling:ThermEvolOn");
/* Flag to switch on extra debug print
* statements in CHIMES.
* 0 --> No extra debug prints.
* 1 --> If CVode returns a non-zero flag (i.e. it produces a
* CVODE error/warning), then we print out everything in
* the ChimesGasVars struct.
*/
cooling->ChimesGlobalVars.chimes_debug =
parser_get_param_int(parameter_file, "CHIMESCooling:chimes_debug");
/* The following CHIMES parameters do not need to be modified
* by the user. These are just hard-coded for now. */
cooling->ChimesGlobalVars.grain_temperature = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:dust_grain_temperature_K", 10.f);
/* Physical velocity divergence isn't easily accessible, so just
* run with static molecular cooling for now. */
cooling->ChimesGlobalVars.StaticMolCooling = 1;
/* Optional parameters to define S and Ca
* relative to Si. */
cooling->S_over_Si_ratio_in_solar = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:S_over_Si_in_solar", 1.f);
cooling->Ca_over_Si_ratio_in_solar = parser_get_opt_param_float(
parameter_file, "CHIMESCooling:Ca_over_Si_in_solar", 1.f);
/* Solar mass fractions of Si, S and Ca are hard-coded here,
* because we might not have access to the COLIBRE cooling tables to
* get these. */
cooling->Si_solar_mass_fraction = 6.64948e-4;
cooling->S_solar_mass_fraction = 3.09171e-4;
cooling->Ca_solar_mass_fraction = 6.41451e-5;
/* CHIMES uses a solar metallicity of 0.0129. */
cooling->Zsol = 0.0129;
/* Dust depletion factors in the solar neighbourhood.
* Taken from Ploeckinger & Schaye (2020) */
cooling->f_dust0_C = 0.34385;
cooling->f_dust0_O = 0.31766;
cooling->f_dust0_Mg = 0.94338;
cooling->f_dust0_Si = 0.94492;
cooling->f_dust0_Ca = 0.9999;
cooling->f_dust0_Fe = 0.99363;
/* Parameter that controls whether to reduce gas-phase metal abundances
* due to dust depletion. */
cooling->colibre_metal_depletion = parser_get_param_int(
parameter_file, "CHIMESCooling:colibre_metal_depletion");
/* Distance from EOS to use thermal equilibrium temperature for subgrid
* props, and to evolve the chemistry in equilibrium. */
cooling->dlogT_EOS = parser_get_param_float(
parameter_file, "CHIMESCooling:delta_logTEOS_subgrid_properties");
/* Flag to use the subgrid density and
* temperature from the COLIBRE cooling tables. */
cooling->use_colibre_subgrid_EOS = parser_get_param_int(
parameter_file, "CHIMESCooling:use_colibre_subgrid_EOS");
/* Flag to set particles to chemical equilibrium
* if they have been heated by feedback. */
cooling->set_FB_particles_to_eqm = parser_get_param_int(
parameter_file, "CHIMESCooling:set_FB_particles_to_eqm");
/* Threshold in dt / t_cool above which we are in the rapid cooling regime.
* If negative, we never use this scheme (i.e. always drift the internal
* energies). */
cooling->rapid_cooling_threshold = parser_get_param_double(
parameter_file, "CHIMESCooling:rapid_cooling_threshold");
/* Properties of the HII region model */
cooling->HIIregion_fion = parser_get_opt_param_float(
parameter_file, "COLIBREFeedback:HIIregion_ionization_fraction", 1.0);
cooling->HIIregion_temp = parser_get_opt_param_float(
parameter_file, "COLIBREFeedback:HIIregion_temperature_K", 1.0e4);
/* Maximum boost factor for reactions on the surfaces
* of dust grains. */
cooling->max_dust_boost_factor = parser_get_param_float(
parameter_file, "CHIMESCooling:max_dust_boost_factor");
/* Density below which the dust boost factor is
* equal to unity, in units of cm^-3. */
cooling->dust_boost_nH_min_cgs = parser_get_param_float(
parameter_file, "CHIMESCooling:dust_boost_nH_min_cgs");
/* Density above which the dust boost factor is
* equal to the maximum, in units of cm^-3. */
cooling->dust_boost_nH_max_cgs = parser_get_param_float(
parameter_file, "CHIMESCooling:dust_boost_nH_max_cgs");
/* Flag to indicate whether or not we redistribute FB heated
* particle's dust in the cooling step */
cooling->destroy_FB_heated_dust = parser_get_param_int(
parameter_file, "CHIMESCooling:destroy_FB_heated_dust");
/* Check that the density thresholds for the dust boost factor are sensible.
* They must both be positive (as we will be taking logs), and nH_min must be
* less than nH_max. */
if ((cooling->dust_boost_nH_min_cgs <= 0.0f) ||
(cooling->dust_boost_nH_max_cgs <= 0.0f) ||
(cooling->dust_boost_nH_min_cgs >= cooling->dust_boost_nH_max_cgs))
error(
"CHIMES Error: dust_boost_nH_min_cgs = %.2e, dust_boost_nH_max_cgs = "
"%.2e. Both must be positive, and dust_boost_nH_min_cgs must be less "
"than dust_boost_nH_max_cgs.",
cooling->dust_boost_nH_min_cgs, cooling->dust_boost_nH_max_cgs);
/* Switch for Hybrid cooling */
cooling->ChimesGlobalVars.hybrid_cooling_mode =
parser_get_param_int(parameter_file, "CHIMESCooling:use_hybrid_cooling");
/* Determine which metal elements to include in the CHIMES network.
* Note that H and He are always included. */
cooling->ChimesGlobalVars.element_included[0] =
parser_get_param_int(parameter_file, "CHIMESCooling:IncludeCarbon");
cooling->ChimesGlobalVars.element_included[1] =
parser_get_param_int(parameter_file, "CHIMESCooling:IncludeNitrogen");
cooling->ChimesGlobalVars.element_included[2] =
parser_get_param_int(parameter_file, "CHIMESCooling:IncludeOxygen");
cooling->ChimesGlobalVars.element_included[3] =
parser_get_param_int(parameter_file, "CHIMESCooling:IncludeNeon");
cooling->ChimesGlobalVars.element_included[4] =
parser_get_param_int(parameter_file, "CHIMESCooling:IncludeMagnesium");
cooling->ChimesGlobalVars.element_included[5] =
parser_get_param_int(parameter_file, "CHIMESCooling:IncludeSilicon");
cooling->ChimesGlobalVars.element_included[6] =
parser_get_param_int(parameter_file, "CHIMESCooling:IncludeSulphur");
cooling->ChimesGlobalVars.element_included[7] =
parser_get_param_int(parameter_file, "CHIMESCooling:IncludeCalcium");
cooling->ChimesGlobalVars.element_included[8] =
parser_get_param_int(parameter_file, "CHIMESCooling:IncludeIron");
/* The COLIBRE cooling tables are needed if we are using hybrid cooling or
* if we are using the subgrid density and temperature on the EOS. */
if ((cooling->ChimesGlobalVars.hybrid_cooling_mode == 1) ||
(cooling->use_colibre_subgrid_EOS == 1)) {
/* Path to colibre cooling table */
parser_get_param_string(parameter_file, "CHIMESCooling:colibre_table_path",
cooling->colibre_table.cooling_table_path);
/* Filebase for split cooling tables */
parser_get_param_string(parameter_file,
"CHIMESCooling:colibre_table_filebase",
cooling->colibre_table.cooling_table_filebase);
/* Set the S/Si and Ca/Si ratios to
* the values already read in for CHIMES. */
cooling->colibre_table.S_over_Si_ratio_in_solar =
cooling->S_over_Si_ratio_in_solar;
cooling->colibre_table.Ca_over_Si_ratio_in_solar =
cooling->Ca_over_Si_ratio_in_solar;
/* Physical constants that we will need, in cgs units */
const float dimension_k[5] = {1, 2, -2, 0, -1};
const float dimension_G[5] = {-1, 3, -2, 0, 0};
cooling->boltzmann_k_cgs =
phys_const->const_boltzmann_k *
units_general_cgs_conversion_factor(us, dimension_k);
cooling->newton_G_cgs =
phys_const->const_newton_G *
units_general_cgs_conversion_factor(us, dimension_G);
cooling->proton_mass_cgs = phys_const->const_proton_mass *
units_cgs_conversion_factor(us, UNIT_CONV_MASS);
cooling->proton_mass_over_boltzmann_k_cgs =
cooling->proton_mass_cgs / cooling->boltzmann_k_cgs;
cooling->colibre_table.log10_kB_cgs = log10(cooling->boltzmann_k_cgs);
cooling->colibre_table.inv_proton_mass_cgs = 1. / cooling->proton_mass_cgs;
cooling->colibre_table.proton_mass_cgs = cooling->proton_mass_cgs;
/* Compute conversion factors */
cooling->internal_energy_to_cgs =
units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
cooling->density_to_cgs =
units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
cooling->number_density_to_cgs =
units_cgs_conversion_factor(us, UNIT_CONV_NUMBER_DENSITY);
cooling->time_to_cgs = units_cgs_conversion_factor(us, UNIT_CONV_TIME);
cooling->internal_energy_over_time_to_cgs = units_cgs_conversion_factor(
us, UNIT_CONV_ENERGY_PER_UNIT_MASS_PER_UNIT_TIME);
cooling->internal_energy_from_cgs = 1. / cooling->internal_energy_to_cgs;
cooling->internal_energy_over_time_from_cgs =
1. / cooling->internal_energy_over_time_to_cgs;
/* Same for the COLIBRE tables */
cooling->colibre_table.pressure_to_cgs =
units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
cooling->colibre_table.internal_energy_to_cgs =
units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
cooling->colibre_table.internal_energy_from_cgs =
1. / cooling->colibre_table.internal_energy_to_cgs;
cooling->colibre_table.number_density_to_cgs =
units_cgs_conversion_factor(us, UNIT_CONV_NUMBER_DENSITY);
cooling->colibre_table.number_density_from_cgs =
1. / cooling->colibre_table.number_density_to_cgs;
cooling->colibre_table.density_to_cgs =
units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
cooling->colibre_table.density_from_cgs =
1. / cooling->colibre_table.density_to_cgs;
/* Pre-factor for the Compton y terms */
cooling->y_compton_factor =
phys_const->const_thomson_cross_section *
phys_const->const_boltzmann_k /
(phys_const->const_speed_light_c * phys_const->const_speed_light_c *
phys_const->const_electron_mass);
/* Recover the minimal energy allowed (in internal units) */
const double u_min = hydro_props->minimal_internal_energy;
/* Convert to CGS units */
cooling->colibre_table.umin_cgs =
u_min * cooling->colibre_table.internal_energy_to_cgs;
/* If running with dust, read if pairing to cooling? */
#if defined(DUST_NONE)
dp->pair_to_cooling = 0;
#else
dp->pair_to_cooling = parser_get_opt_param_int(
parameter_file, "DustEvolution:pair_to_cooling", 1);
if (cooling->colibre_metal_depletion && dp->pair_to_cooling) {
error(
"May not run with DustEvolution:pair_to_cooling:1 and "
"CHIMESCooling:colibre_metal_depletion:1 simultaneously, as "
"this double-counts dust depletion effects. Please try switching "
"one of these flags off in the configuration file.");
}
#endif
/* Read the Colibre table. */
message("Reading Colibre cooling table.");
get_cooling_redshifts(&(cooling->colibre_table));
/* Read in cooling table header */
char fname[colibre_table_path_name_length + colibre_table_filebase_length +
12];
sprintf(fname, "%s/%s_z0.00.hdf5",
cooling->colibre_table.cooling_table_path,
cooling->colibre_table.cooling_table_filebase);
read_cooling_header(fname, &(cooling->colibre_table));
/* Allocate memory for the tables */
allocate_cooling_tables(&(cooling->colibre_table));
/* Set the redshift indices to invalid values */
cooling->colibre_table.z_index = -10;
/* set previous_z_index and to last value of redshift table*/
cooling->colibre_table.previous_z_index =
cooling->colibre_table.number_of_redshifts - 2;
}
/* Set redshift to a very high value, just
* while we initialise the CHIMES module.
* It will be set to the correct redshift in
* the cooling_update() routine. */
cooling->ChimesGlobalVars.redshift = 1000.0;
if ((cooling->UVB_dependence_choice == UVB_zdep_cross_secs_and_eqm) &&
(cooling->UV_field_choice == UV_field_COLIBRE) &&
(cooling->Shielding_choice == Shielding_COLIBRE)) {
/* If we are using the redshift-dependent CHIMES eqm abundances with the
* COLIBRE radiation field and shielding models, we need to check that the
* parameters used in the tables match the hard-coded values. */
char fname_chimes[520];
sprintf(fname_chimes, "%s/z%.3f_eqm.hdf5",
cooling->ChimesGlobalVars.EqAbundanceTablePath, 0.0);
/* Implement better version control check but not while we are
* still testing */
/* chimes_check_colibre_eqm_tables_header(fname_chimes, cooling); */
}
/* Initialise the CHIMES module. */
message("Initialising CHIMES cooling module.");
init_chimes(&cooling->ChimesGlobalVars);
/* Check that the size of the CHIMES network
* determined by the various Include
* parameters matches the size specified
* through the configure script. */
if (cooling->ChimesGlobalVars.totalNumberOfSpecies != CHIMES_NETWORK_SIZE)
error(
"The size of the CHIMES network based on the Include "
"parameters is %d. However, Swift was configured and compiled using "
"--with-cooling=CHIMES_%d. If you are sure that you have "
"correctly set which metals to include in the network in the parameter "
"file, then you will need to re-compile Swift using ./configure "
"--with-cooling=CHIMES_%d.",
cooling->ChimesGlobalVars.totalNumberOfSpecies, CHIMES_NETWORK_SIZE,
cooling->ChimesGlobalVars.totalNumberOfSpecies);
/* Construct the array of species names
* for the reduced network. */
for (int idx = 0; idx < CHIMES_TOTSIZE; idx++) {
if (cooling->ChimesGlobalVars.speciesIndices[idx] >= 0)
strcpy(cooling->chimes_species_names_reduced[cooling->ChimesGlobalVars
.speciesIndices[idx]],
chimes_species_names[idx]);
}
if (cooling->ChimesGlobalVars.hybrid_cooling_mode == 1) {
/* Create data structure for hybrid cooling,
* and store pointer to the Colibre table. */
cooling->ChimesGlobalVars.hybrid_data =
(void *)malloc(sizeof(struct global_hybrid_data_struct));
struct global_hybrid_data_struct *myData;
myData = (struct global_hybrid_data_struct *)
cooling->ChimesGlobalVars.hybrid_data;
myData->table = &(cooling->colibre_table);
myData->Zsol = cooling->Zsol;
/* Set the hybrid cooling function pointers. */
cooling->ChimesGlobalVars.hybrid_cooling_fn =
&colibre_metal_cooling_rate_temperature;
cooling->ChimesGlobalVars.allocate_gas_hybrid_data_fn =
&chimes_allocate_gas_hybrid_data;
cooling->ChimesGlobalVars.free_gas_hybrid_data_fn =
&chimes_free_gas_hybrid_data;
}
}
/**
* @brief Prints the properties of the cooling model to stdout.
*
* @param cooling #cooling_function_data struct.
*/
void cooling_print_backend(const struct cooling_function_data *cooling) {
message("Cooling function is 'SPLITCHIMES'.");
}
/**
* @brief Common operations performed on the cooling function at a
* given time-step or redshift.
*
* @param cosmo The current cosmological model.
* @param pressure_floor Properties of the pressure floor.
* @param cooling The #cooling_function_data used in the run.
* @param s The #space containing all the particles.
* @param dustevo The dust evolution properties
* @param time The current system time
*/
void cooling_update(const struct cosmology *cosmo,
const struct pressure_floor_props *pressure_floor,
struct cooling_function_data *cooling, struct space *s,
struct dustevo_props *dustevo, const double time) {
/* Update redshift stored in ChimesGlobalVars. */
cooling->ChimesGlobalVars.redshift = cosmo->z;
/* Update T_CMB in CHIMES to current redshift. */
cooling->ChimesGlobalVars.cmb_temperature =
(ChimesFloat)cooling->T_CMB_0 * (1.0 + cosmo->z);
/* Update redshift-dependent UVB. */
if (cooling->UVB_dependence_choice == UVB_zdep_cross_secs ||
cooling->UVB_dependence_choice == UVB_zdep_cross_secs_and_eqm)
interpolate_redshift_dependent_UVB(&cooling->ChimesGlobalVars);
/* The COLIBRE cooling tables are needed
* if we are using hybrid cooling or
* if we are using the subgrid density
* and temperature on the EOS. */
if ((cooling->ChimesGlobalVars.hybrid_cooling_mode == 1) ||
(cooling->use_colibre_subgrid_EOS == 1)) {
/* Current redshift */
const float redshift = cosmo->z;
/* What is the current table index along the redshift axis? */
int z_index = -1;
float dz = 0.f;
get_redshift_index(redshift, &z_index, &dz, &(cooling->colibre_table));
cooling->colibre_table.dz = dz;
/* Do we already have the correct tables loaded? */
if (cooling->colibre_table.z_index == z_index) return;
/* Normal case: two tables bracketing the current z */
const int low_z_index = z_index;
const int high_z_index = z_index + 1;
get_cooling_table(&(cooling->colibre_table), low_z_index, high_z_index,
dustevo->logfD, dustevo->pair_to_cooling);
/* Store the currently loaded index */
cooling->colibre_table.z_index = z_index;
}
}
/**
* @brief Update ChimesGasVars structure.
*
* Updates the ChimesGasVars structure with the various
* thermodynamics quantities of the gas particle that
* will be needed for the CHIMES chemistry solver.
*
* @param u_cgs The internal energy, in cgs units.
* @param phys_const The physical constants in internal units.
* @param us The internal system of units.
* @param cosmo The current cosmological model.
* @param hydro_properties the hydro_props struct
* @param floor_props Properties of the entropy floor.
* @param cooling The #cooling_function_data used in the run.
* @param p Pointer to the particle data.
* @param xp Pointer to the extended particle data.
* @param ChimesGasVars CHIMES gasVariables structure.
* @param dt_cgs The cooling time-step of this particle.
*/
void chimes_update_gas_vars(
const double u_cgs, const struct phys_const *phys_const,
const struct unit_system *us, const struct cosmology *cosmo,
const struct hydro_props *hydro_properties,
const struct entropy_floor_properties *floor_props,
const struct cooling_function_data *cooling, const struct dustevo_props *dp,
const struct part *p, const struct xpart *xp,
struct gasVariables *ChimesGasVars, const float dt_cgs) {
/* Get the mean molecular weight */
const double mu = chimes_mu(cooling, hydro_properties, p, xp);
/* Get the Hydrogen mass fraction */
const ChimesFloat XH = get_hydrogen_mass_fraction(p, hydro_properties);
/* Get the metallicity (metal mass fraction / solar metal mass frac.) */
#if defined(CHEMISTRY_COLIBRE) || defined(CHEMISTRY_EAGLE)
ChimesGasVars->metallicity =
(ChimesFloat)chemistry_get_total_metal_mass_fraction_for_cooling(p);
ChimesGasVars->metallicity /= cooling->Zsol;
ChimesGasVars->metallicity = max(ChimesGasVars->metallicity, FLT_MIN);
#else
/* Without COLIBRE or EAGLE chemistry, the metal abundances are unavailable.
* --> Set to primordial abundances. */
ChimesGasVars->metallicity = 0.0;
#endif /* CHEMISTRY_COLIBRE || CHEMISTRY_EAGLE */
dust_check_for_NaNs(p)
/* Get the dust ratio */
ChimesGasVars->dust_ratio =
get_dust_ratio(p, cooling, dp, ChimesGasVars->metallicity);
/* 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);
const double u_floor_cgs = u_floor * cooling->internal_energy_to_cgs;
/* Set the temperature floor */
if (xp->tracers_data.HIIregion_timer_gas > 0.) {
/* Particles is an HII region. */
/* Temperature on the entropy floor, but using the HII region mean
* molecular weight. */
const double mu_HII = get_mu_HII_region(cooling, XH);
const double T_floor = get_T_from_u_cgs(u_floor_cgs, mu_HII, cooling);
/* Set the temperature floor in CHIMES to be either T_floor or
* HIIregion_temp, whichever is higher. */
ChimesGasVars->TempFloor =
(ChimesFloat)chimes_max(T_floor, cooling->HIIregion_temp);
/* Check TempFloor doesn't go below minimal_temperature. */
ChimesGasVars->TempFloor =
chimes_max(ChimesGasVars->TempFloor,
(ChimesFloat)hydro_properties->minimal_temperature);
} else {
/* Particle is *not* an HII region. */
/* Temperature on the entropy floor using current mean molecular weight. */
const double T_floor = get_T_from_u_cgs(u_floor_cgs, mu, cooling);
/* If using the entropy floor, passing it as a negative number means
* CHIMES will interpret it as a minimum thermal energy, in erg/cm^3,
* rather than a minimum temperature. */
if (T_floor < hydro_properties->minimal_temperature) {
ChimesGasVars->TempFloor =
(ChimesFloat)hydro_properties->minimal_temperature;
} else {
ChimesGasVars->TempFloor =
(ChimesFloat)(-u_floor_cgs * rho * cooling->density_to_cgs);
}
}
ChimesGasVars->temp_floor_mode = 1;
/* Select energy to use based on where we sit w.r.t. the entropy floor */
double u_actual_cgs;
if (u_cgs < u_floor_cgs) {
/* Particle is below the entropy floor.
*
* Set internal energy to the floor and chemistry will be evolved in
* equilibrium. */
u_actual_cgs = u_floor_cgs;
ChimesGasVars->ForceEqOn = 1;
} else if (u_cgs < exp10(cooling->dlogT_EOS) * u_floor_cgs) {
/* Particle is above the entropy floor, but close enough that we will need
* to evolve the chemistry in equilibrium. */
u_actual_cgs = u_cgs;
ChimesGasVars->ForceEqOn = 1;
} else {
/* Particle is well above the entropy floor.
*
* Evolve chemistry as usual, according to the user-provided parameter. */
u_actual_cgs = u_cgs;
ChimesGasVars->ForceEqOn = cooling->ChemistryEqmMode;
}
/* Now convert this to a temperature */
ChimesGasVars->temperature =
(ChimesFloat)get_T_from_u_cgs(u_actual_cgs, mu, cooling);
/* Recover the Hydrogen number density */
const double nH =
hydro_get_physical_density(p, cosmo) * XH / phys_const->const_proton_mass;
const double nH_cgs = (nH * cooling->number_density_to_cgs);
ChimesGasVars->nH_tot = (ChimesFloat)nH_cgs;
/* Recover the dust boost factor to use based on the H number density */
ChimesGasVars->dust_boost_factor = get_dust_boost_factor(nH_cgs, cooling);
ChimesGasVars->hydro_timestep = (ChimesFloat)dt_cgs;
ChimesGasVars->constant_heating_rate = 0.0;
ChimesGasVars->ThermEvolOn = cooling->ThermEvolOn;
ChimesGasVars->divVel = 0.0;
/* N_ref is used by both the COLIBRE ISRF and the COLIBRE shielding length,
* and to scale metal depletion. Based on Ploeckinger & Schaye (2020). */
double N_ref_cgs = 0.0;
if ((cooling->UV_field_choice == UV_field_COLIBRE) ||
(cooling->Shielding_choice == Shielding_COLIBRE) ||
(cooling->colibre_metal_depletion)) {
N_ref_cgs = calculate_colibre_N_ref(phys_const, us, cosmo, cooling,
hydro_properties, p, xp,
ChimesGasVars->temperature);
}
/* Get the CR rate we use */
ChimesGasVars->cr_rate =
get_CR_rate(N_ref_cgs, nH_cgs, cosmo, phys_const, cooling);
/* Set the radiation field by copying over the sectrum parameters from the
* chimes tables. */
if (cooling->UV_field_choice == UV_field_user_defined) {
/* Single, constant radiation field. */
ChimesGasVars->G0_parameter[0] = chimes_table_spectra.G0_parameter[0];
ChimesGasVars->H2_dissocJ[0] = chimes_table_spectra.H2_dissocJ[0];
ChimesGasVars->isotropic_photon_density[0] =
chimes_table_spectra.isotropic_photon_density[0];
ChimesGasVars->isotropic_photon_density[0] *=
cooling->rad_field_norm_factor;
} else if (cooling->UV_field_choice == UV_field_COLIBRE) {
/* Extra-galactic UVB */
ChimesGasVars->G0_parameter[0] = chimes_table_spectra.G0_parameter[0];
ChimesGasVars->H2_dissocJ[0] = chimes_table_spectra.H2_dissocJ[0];
ChimesGasVars->isotropic_photon_density[0] =
chimes_table_spectra.isotropic_photon_density[0];
/* ISRF, scales with Nref */
ChimesGasVars->G0_parameter[1] = chimes_table_spectra.G0_parameter[1];
ChimesGasVars->H2_dissocJ[1] = chimes_table_spectra.H2_dissocJ[1];
const double log_isrf_min =
log10(chimes_table_spectra.isotropic_photon_density[1] *
cooling->isrf_norm_at_transition *
pow(cooling->nref_column_density_min_cgs /
cooling->isrf_transition_column_density_cgs,
cooling->isrf_powerlaw_slope_1));
const double isrf_coldensratio =
N_ref_cgs / cooling->isrf_transition_column_density_cgs;
double isophotdens;
if (isrf_coldensratio < 1.) {
/* N_ref < N_trans: power law index 1 */
isophotdens = chimes_table_spectra.isotropic_photon_density[1] *
cooling->isrf_norm_at_transition *
pow(isrf_coldensratio, cooling->isrf_powerlaw_slope_1);
} else {
/* N_ref >=N_trans: power law index 2 */
isophotdens = chimes_table_spectra.isotropic_photon_density[1] *
cooling->isrf_norm_at_transition *
pow(isrf_coldensratio, cooling->isrf_powerlaw_slope_2);
}
const double log_isophotdens_cgs = log10(isophotdens);
/* Now, apply some saturation */
/* low density cutoff between isrf_nmax and isrf_nmin
* with k regulating the steepness */
const double isrf_nmax_over_dens =
cosmo->critical_density * cooling->isrf_cutoff_over_density *
cooling->nref_XH_const / phys_const->const_proton_mass;
const double isrf_nmax_over_dens_cgs =
isrf_nmax_over_dens * cooling->number_density_to_cgs;
/* Max density */
const double isrf_nmax =
max(isrf_nmax_over_dens_cgs, cooling->isrf_cutoff_phys_density_cgs);
/* Min density: (1. / 10^(cutoff_width_dex)) * isrf_max */
const double isrf_nmin = exp10(-cooling->isrf_cutoff_width_dex) * isrf_nmax;
const double numerator = log_isophotdens_cgs - log_isrf_min;
const double argument =
sqrt(isrf_nmin / isrf_nmax) * isrf_nmax / ChimesGasVars->nH_tot;
const double denominator =
1. + pow(argument, -cooling->isrf_cutoff_steepness_k);
const double log_isophotdens =
log_isophotdens_cgs - numerator / denominator;
ChimesGasVars->isotropic_photon_density[1] = exp10(log_isophotdens);
}
double N_scale_rad_cgs = N_ref_cgs;
if (cooling->colibre_metal_depletion) {
/* Scale dust_ratio by N_scale_rad, but only if N_scale_rad < N_H0 */
if (N_scale_rad_cgs < cooling->N_H0)
ChimesGasVars->dust_ratio *=
(ChimesFloat)pow(N_scale_rad_cgs / cooling->N_H0, 1.4);
}
/* Set shielding length */
if (cooling->Shielding_choice == Shielding_none) {
ChimesGasVars->cell_size = 0;
} else if (cooling->Shielding_choice == Shielding_Jeans_length) {
/* Jeans length */
ChimesGasVars->cell_size = (ChimesFloat)sqrt(
M_PI * hydro_gamma * cooling->boltzmann_k_cgs *
((double)ChimesGasVars->temperature) /
(mu * cooling->newton_G_cgs *
(cooling->proton_mass_cgs * ((double)ChimesGasVars->nH_tot) / XH) *
cooling->proton_mass_cgs));
ChimesGasVars->cell_size *= cooling->shielding_length_factor;
if ((cooling->max_shielding_length_cm > 0.0) &&
(ChimesGasVars->cell_size > cooling->max_shielding_length_cm))
ChimesGasVars->cell_size = (ChimesFloat)cooling->max_shielding_length_cm;
} else if (cooling->Shielding_choice == Shielding_COLIBRE) {
double shielding_length_cgs = cooling->shielding_length_factor * N_ref_cgs /
((double)ChimesGasVars->nH_tot);
/* limit shielding length by redshift dependent maximum value */
const double numerator =
(cooling->lshmax_lowz_kpc - cooling->lshmax_highz_kpc);
const double argument = sqrt((cooling->lshmax_transition_redshift -
cooling->lshmax_transition_width_dz) *
(cooling->lshmax_transition_redshift)) /
fmax(cosmo->z, 1.e-10);
const double denominator =
1. + pow(fmin(argument, 10.), cooling->lshmax_transition_steepness_k);
const double shielding_length_max_kpc =
cooling->lshmax_lowz_kpc - numerator / denominator;
double shielding_length_max_cgs =
shielding_length_max_kpc * 1.0e3 *
units_cgs_conversion_factor(us, UNIT_CONV_LENGTH) *
phys_const->const_parsec;
ChimesGasVars->cell_size =
(ChimesFloat)(fmin(shielding_length_cgs, shielding_length_max_cgs));
}
/* Doppler broadening parameter, for H2 self-shielding.
* Value is set in the parameter file. */
ChimesGasVars->doppler_broad = M_SQRT2 * cooling->sigma_turb_km_p_s;
ChimesGasVars->InitIonState = cooling->InitIonState;
/* If using hybrid cooling, we need to set the abundance_ratio array using
* the corresponding routine from COLIBRE. */
if (cooling->ChimesGlobalVars.hybrid_cooling_mode == 1) {
struct gas_hybrid_data_struct *myData;
myData = (struct gas_hybrid_data_struct *)ChimesGasVars->hybrid_data;
abundance_ratio_to_solar(p, &(cooling->colibre_table),
myData->abundance_ratio);
}
}
/**
* @brief Update CHIMES element abundances.
*
* Updates the element abundances based on the particle's current mass
* fractions. If the element abundances have changed since the last time this
* routine was called, for example due to enrichment or metal diffusion, then
* the ion and molecule abundances will now be inconsistent with their
* corresponding total element abundances. Therefore, if mode == 1, we call the
* check_constraint_equations() routine from the CHIMES module, which adds all
* species of each element and compares to the total abundance of that element.
* If they are inconsistent, the species abundances are adjusted, preserving
* their relative fractions. If a metal is now present that was not present
* before (e.g. if a particle was primordial but has recently been enriched),
* then that metal is introduced as fully neutral.
*
* We do not call check_constraint_equations() if mode == 0 (for example, when
* we initialise the abundance arrays for the first time).
*
* @param phys_const #phys_const data structure.
* @param us The internal system of units.
* @param cosmo #cosmology data structure.
* @param cooling #cooling_function_data struct.
* @param hydro_properties The properties of the hydro scheme
* @param p #part data.
* @param xp Pointer to the #xpart data.
* @param ChimesGasVars CHIMES gasVariables structure.
* @param u_0 internal energy (in code units).
* @param mode Set to zero if particle not fully initialised.
*/
void chimes_update_element_abundances(
const struct phys_const *phys_const, const struct unit_system *us,
const struct cosmology *cosmo, const struct cooling_function_data *cooling,
const struct hydro_props *hydro_props, const struct part *p,
const struct xpart *xp, struct gasVariables *ChimesGasVars,
const double u_0, const int mode) {
const struct globalVariables *ChimesGlobalVars = &cooling->ChimesGlobalVars;
/* Special case when there is no chemistry */
#if !defined(CHEMISTRY_COLIBRE) && !defined(CHEMISTRY_EAGLE)
/* Without COLIBRE or EAGLE chemistry,the metal abundances are unavailable.
* --> Set to primordial abundances. */
ChimesGasVars->element_abundances[0] = 0.0833; /* He */
for (int i = 1; i < 10; i++) ChimesGasVars->element_abundances[i] = 0.0;
if (mode == 1) check_constraint_equations(ChimesGasVars, ChimesGlobalVars);
return;
#endif /* !CHEMISTRY_COLIBRE && !CHEMISTRY_EAGLE */
/* Normal case */
/* Get element mass fractions */
const float *const metal_fraction =
chemistry_get_metal_mass_fraction_for_cooling(p);
const float XH = metal_fraction[chemistry_element_H];
ChimesGasVars->element_abundances[0] =
(ChimesFloat)metal_fraction[chemistry_element_He] / (4.0 * XH);
ChimesGasVars->element_abundances[1] =
(ChimesFloat)metal_fraction[chemistry_element_C] / (12.0 * XH);
ChimesGasVars->element_abundances[2] =
(ChimesFloat)metal_fraction[chemistry_element_N] / (14.0 * XH);
ChimesGasVars->element_abundances[3] =
(ChimesFloat)metal_fraction[chemistry_element_O] / (16.0 * XH);
ChimesGasVars->element_abundances[4] =
(ChimesFloat)metal_fraction[chemistry_element_Ne] / (20.0 * XH);
ChimesGasVars->element_abundances[5] =
(ChimesFloat)metal_fraction[chemistry_element_Mg] / (24.0 * XH);
ChimesGasVars->element_abundances[6] =
(ChimesFloat)metal_fraction[chemistry_element_Si] / (28.0 * XH);
ChimesGasVars->element_abundances[9] =
(ChimesFloat)metal_fraction[chemistry_element_Fe] / (56.0 * XH);
ChimesGasVars->element_abundances[7] =
(ChimesFloat)metal_fraction[chemistry_element_Si] *
cooling->S_over_Si_ratio_in_solar *
(cooling->S_solar_mass_fraction / cooling->Si_solar_mass_fraction) /
(32.0 * XH);
ChimesGasVars->element_abundances[8] =
(ChimesFloat)metal_fraction[chemistry_element_Si] *
cooling->Ca_over_Si_ratio_in_solar *
(cooling->Ca_solar_mass_fraction / cooling->Si_solar_mass_fraction) /
(40.0 * XH);
/* When chimes_update_element_abundances() is called on a particle that is
* being initialised for the first time, the chimes_abundance array has not
* yet been initialised. We therefore set mu = 1. */
const double mu = (mode == 0) ? 1.0 : chimes_mu(cooling, hydro_props, p, xp);
if (cooling->colibre_metal_depletion) {
const double u_cgs = u_0 * cooling->internal_energy_to_cgs;
const double T = get_T_from_u_cgs(u_cgs, mu, cooling);
/* Reference column density */
const double N_ref = calculate_colibre_N_ref(phys_const, us, cosmo, cooling,
hydro_props, p, xp, T);
/* Dust depletion factor */
const double factor = fmin(pow(N_ref / cooling->N_H0, 1.4), 1.0);
/* Reduce gas-phase metal abundances due to dust depletion. */
ChimesGasVars->element_abundances[1] *=
(ChimesFloat)(1.0 - (cooling->f_dust0_C * factor));
ChimesGasVars->element_abundances[3] *=
(ChimesFloat)(1.0 - (cooling->f_dust0_O * factor));
ChimesGasVars->element_abundances[5] *=
(ChimesFloat)(1.0 - (cooling->f_dust0_Mg * factor));
ChimesGasVars->element_abundances[6] *=
(ChimesFloat)(1.0 - (cooling->f_dust0_Si * factor));
ChimesGasVars->element_abundances[8] *=
(ChimesFloat)(1.0 - (cooling->f_dust0_Ca * factor));
ChimesGasVars->element_abundances[9] *=
(ChimesFloat)(1.0 - (cooling->f_dust0_Fe * factor));
}
/* Zero the abundances of any elements that not included in the network. */
for (int i = 0; i < 9; i++)
ChimesGasVars->element_abundances[i + 1] *=
(ChimesFloat)ChimesGlobalVars->element_included[i];
if (mode == 1) check_constraint_equations(ChimesGasVars, ChimesGlobalVars);
}
/**
* @brief Apply the cooling function to a particle.
*
* We use the CHIMES module to integrate the cooling rate
* and chemical abundances over the time-step.
*
* @param phys_const The physical constants in internal units.
* @param us The internal system of units.
* @param cosmo The current cosmological model.
* @param hydro_properties the hydro_props struct
* @param floor_props Properties of the entropy floor.
* @param pressure_floor Properties of the pressure floor.
* @param cooling The #cooling_function_data used in the run.
* @param p Pointer to the particle data.
* @param xp Pointer to the extended particle data.
* @param dt The cooling time-step of this particle.
* @param dt_therm The hydro time-step of this particle.
* @param time Time since Big Bang
*/
void cooling_cool_part(const struct phys_const *phys_const,
const struct unit_system *us,
const struct cosmology *cosmo,
const struct hydro_props *hydro_properties,
const struct entropy_floor_properties *floor_props,
const struct pressure_floor_props *pressure_floor,
const struct cooling_function_data *cooling,
const struct dustevo_props *dp, struct part *p,
struct xpart *xp, const float dt, const float dt_therm,
const double time) {
/* No cooling happens over zero time */
if (dt == 0.) {
/* But we still set the subgrid properties to a valid state */
cooling_set_subgrid_properties(phys_const, us, cosmo, hydro_properties,
floor_props, cooling, dp, p, xp);
return;
}
ticks tic_chimes = getticks();
/* Verify whether we are still in an HII region or not
* (i.e. has the counter expired?) */
if ((time > xp->tracers_data.HIIregion_timer_gas) &&
(xp->tracers_data.HIIregion_timer_gas > 0.)) {
/* Not anymore an HII region */
xp->tracers_data.HIIregion_timer_gas = -1.;
xp->tracers_data.HIIregion_starid = -1;
}
/* CHIMES structures. */
const struct globalVariables *ChimesGlobalVars = &cooling->ChimesGlobalVars;
/* Allocate memory to arrays within ChimesGasVars. */
struct gasVariables ChimesGasVars;
allocate_gas_abundances_memory(&ChimesGasVars, ChimesGlobalVars);
/* Copy abundances over from xp to ChimesGasVars. */
for (int i = 0; i < ChimesGlobalVars->totalNumberOfSpecies; i++)
ChimesGasVars.abundances[i] = xp->cooling_data.chimes_abundances[i];
/* Get internal energy at the last kick step */
const float u_start = hydro_get_physical_internal_energy(p, xp, cosmo);
/* 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
* increase) */
double u_0 = (u_start + hydro_du_dt * dt_therm);
/* Prior to computing updated element abundances, do we destroy dust? */
if (xp->cooling_data.heated_by_FB && cooling->destroy_FB_heated_dust) {
dust_redistribute_masses_cooling_step(p, dp);
dust_check_for_NaNs(p);
}
/* Update element abundances from metal mass fractions.
*
* We need to do this here, and not later on in chimes_update_gas_vars(),
* because the element abundances need to be set in ChimesGasVars before we
* can calculate the mean molecular weight. */
chimes_update_element_abundances(phys_const, us, cosmo, cooling,
hydro_properties, p, xp, &ChimesGasVars, u_0,
/*mode=*/1);
/* Check for minimal energy limit at the start of the step.
*
* Note that, to maintain consistency with the temperature
* floor that is imposed within CHIMES, we compute the minimal energy from the
* minimal temperature based on the particle's actual mean molecular weight
* mu, rather than an assumed constant mu.
*/
const double mu_start = chimes_mu(cooling, hydro_properties, p, xp);
const double minimal_internal_energy_start = get_u_internal_from_T(
hydro_properties->minimal_temperature, mu_start, phys_const);
u_0 = max(u_0, minimal_internal_energy_start);
/* Convert to CGS units */
const double u_0_cgs = u_0 * cooling->internal_energy_to_cgs;
const double dt_cgs = dt * cooling->time_to_cgs;
/* Update the ChimesGasVars structure with the
* particle's thermodynamic variables. */
chimes_update_gas_vars(u_0_cgs, phys_const, us, cosmo, hydro_properties,
floor_props, cooling, dp, p, xp, &ChimesGasVars,
dt_cgs);
/* Store the updated dust boost factor */
xp->cooling_data.dust_boost_factor = ChimesGasVars.dust_boost_factor;
/* Special treatment for particles recently heated by feedback */
if (xp->cooling_data.heated_by_FB && cooling->set_FB_particles_to_eqm) {
/* Set the particle abundances to equilibrium values */
cooling_set_FB_particle_chimes_abundances(&ChimesGasVars, cooling);
}
if (xp->cooling_data.heated_by_FB &&
(cooling->set_FB_particles_to_eqm || cooling->destroy_FB_heated_dust)) {
/* Reset the heated by FB flag, Particles that were heated are
treated as regular particles from now on.*/
xp->cooling_data.heated_by_FB = 0;
}
/* If particle is an HII region, immediately heat it up to 1e4 K (if it was
* colder), and force it to be evolved with equilibrium cooling. */
if (is_HII_region(xp)) {
ChimesGasVars.temperature = chimes_max(
ChimesGasVars.temperature, (ChimesFloat)cooling->HIIregion_temp);
ChimesGasVars.ForceEqOn = 1;
}
/* Check that the temperature passed to CHIMES is not unreasonably high. */
if (ChimesGasVars.temperature > 1.0e12) {
chimes_print_gas_vars(stderr, &ChimesGasVars, ChimesGlobalVars);
error(
"CHIMES ERROR: the temperature that is being passed to CHIMES is %.6e. "
"This is too high!",
ChimesGasVars.temperature);
}
/* Make a backup of the initial gas vars. */
struct gasVariables ChimesGasVars_save;
ChimesGasVars_save = ChimesGasVars;
allocate_gas_abundances_memory(&ChimesGasVars_save, ChimesGlobalVars);
memcpy(ChimesGasVars_save.isotropic_photon_density,
ChimesGasVars.isotropic_photon_density,
sizeof(ChimesFloat) * ChimesGlobalVars->N_spectra);
memcpy(ChimesGasVars_save.G0_parameter, ChimesGasVars.G0_parameter,
sizeof(ChimesFloat) * ChimesGlobalVars->N_spectra);
memcpy(ChimesGasVars_save.H2_dissocJ, ChimesGasVars.H2_dissocJ,
sizeof(ChimesFloat) * ChimesGlobalVars->N_spectra);
memcpy(ChimesGasVars_save.abundances, ChimesGasVars.abundances,
sizeof(ChimesFloat) * ChimesGlobalVars->totalNumberOfSpecies);
/* Call CHIMES to integrate the chemistry and cooling over the time-step. */
const int chimes_return_code =
chimes_network(&ChimesGasVars, ChimesGlobalVars);
/* Check whether the abundances became unphysical before enforcing the
* constraint equations within CHIMES.
* If they did, repeat CHIMES with equilibrium cooling. */
if (chimes_return_code == 1) {
if ((strcmp(ChimesGlobalVars->EqAbundanceTablePath, "None") == 0) ||
(strcmp(ChimesGlobalVars->EqAbundanceTablePath, "none") == 0)) {
/* Special case if no tables were provided... */
error(
"CHIMES ERROR: chimes_network() returned an error code. This "
"suggests we might need to use stricter tolerances in CHIMES.");
}
/* Be verbose about the failure! */
chimes_print_gas_vars(stderr, &ChimesGasVars_save, ChimesGlobalVars);
chimes_print_gas_vars(stderr, &ChimesGasVars, ChimesGlobalVars);
fprintf(stderr,
"CHIMES Warning: chimes_network() returned an error code. We will "
"repeat the CHIMES integration with equilibrium cooling.");
/* Copy over original CHIMES abundances and temperature. */
memcpy(ChimesGasVars.abundances, ChimesGasVars_save.abundances,
sizeof(ChimesFloat) * ChimesGlobalVars->totalNumberOfSpecies);
ChimesGasVars.temperature = ChimesGasVars_save.temperature;
/* Use equilibrium abundances. */
ChimesGasVars.ForceEqOn = 1;
/* Integrate CHIMES in equilibrium.
* (Note we assume this does *not* fail) */
chimes_network(&ChimesGasVars, ChimesGlobalVars);
}
/* Check that the temperature returned is not unreasonably high. */
if (ChimesGasVars.temperature > 1.0e12) {
chimes_print_gas_vars(stderr, &ChimesGasVars_save, ChimesGlobalVars);
chimes_print_gas_vars(stderr, &ChimesGasVars, ChimesGlobalVars);
error(
"CHIMES ERROR: the temperature that is being returned from CHIMES is "
"%.6e. This is too high!",
ChimesGasVars.temperature);
}
/* Check that radiative heating alone has not increased the temperature to
* an unphysical level.
*
* Radiative heating alone should not increase the temperature
* above 1e5 K, so above this we should have net cooling from CHIMES.
* Although changes in mean molecular weight could
* make the temperature increase by a factor of a few. */
const float T_rad_heat_max = 1.0e5f;
const float T_relative_change_limit = 10.0f;
if ((ChimesGasVars.temperature > T_rad_heat_max) &&
(ChimesGasVars.temperature / ChimesGasVars_save.temperature >
T_relative_change_limit)) {
chimes_print_gas_vars(stderr, &ChimesGasVars_save, ChimesGlobalVars);
chimes_print_gas_vars(stderr, &ChimesGasVars, ChimesGlobalVars);
error(
"CHIMES ERROR: After the CHIMES integration the temperature has "
"increased from %.6e to %.6e. Radiative heating alone should not "
"increase the temperature above 1e5 K, so this suggests we might need "
"to use stricter tolerances in CHIMES.",
ChimesGasVars_save.temperature, ChimesGasVars.temperature);
}
/* If particle is in an HII region, manually set it to be ionised,
* according to the HIIregion_fion parameter. */
if (is_HII_region(xp)) {
cooling_set_HIIregion_chimes_abundances(&ChimesGasVars, cooling);
}
/* Calculate net cooling rate at the end of the time-step.
*
* If the particle is in an HII region, set the cooling rate to zero.
* If the particle is on the EOS, the cooling rate will be
* overwritten when we compute the subgrid properties below. */
const float rho_phys_cgs =
hydro_get_physical_density(p, cosmo) * cooling->density_to_cgs;
if (is_HII_region(xp)) {
xp->cooling_data.net_cooling_rate = 0.0f;
} else {
/* cooling rate in erg/cm^3/s. */
float net_cooling_rate_cgs = (float)chimes_cooling_rate_external(
&ChimesGasVars, ChimesGlobalVars, /*mode=*/0);
/* Convert to erg/g/s */
net_cooling_rate_cgs /= rho_phys_cgs;
/* Convert to code units. */
xp->cooling_data.net_cooling_rate =
net_cooling_rate_cgs * cooling->internal_energy_over_time_from_cgs;
}
/* Recover the Hydrogen mass fraction */
const double XH = get_hydrogen_mass_fraction(p, hydro_properties);
/* Copy abundances from ChimesGasVars back to the particle. */
for (int i = 0; i < ChimesGlobalVars->totalNumberOfSpecies; i++)
xp->cooling_data.chimes_abundances[i] = ChimesGasVars.abundances[i];
/* Calculate mean molecular weight at the end of the time step. */
const float mu_end = chimes_mu(cooling, hydro_properties, p, xp);
/* Compute the internal energy at the end of the
* step using the final temperature from CHIMES. */
const double u_final_cgs =
get_u_cgs_from_T((double)ChimesGasVars.temperature, mu_end, cooling);
/* Convert back to internal units */
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 */
/* Recompute new minimal internal energy
* from the new mean molecular weight. */
const double minimal_internal_energy_end = get_u_internal_from_T(
hydro_properties->minimal_temperature, mu_end, phys_const);
u_final = max(u_final, minimal_internal_energy_end);
/* 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);
u_final = max(u_final, u_floor);
/* 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);
/* Determine if we are in the slow- or rapid-cooling regime,
* by comparing dt / t_cool to the rapid_cooling_threshold.
*
* Note that dt / t_cool = fabs(delta_u) / u_start. */
const double dt_over_t_cool = fabs(delta_u) / max(u_start, u_floor);
/* If rapid_cooling_threshold < 0, always use the slow-cooling regime. */
if ((cooling->rapid_cooling_threshold >= 0.0) &&
(dt_over_t_cool >= cooling->rapid_cooling_threshold)) {
/* Rapid-cooling regime. */
/* Update the particle's u and du/dt */
hydro_set_physical_internal_energy(p, xp, cosmo, u_final);
hydro_set_drifted_physical_internal_energy(p, cosmo, pressure_floor,
u_final);
hydro_set_physical_internal_energy_dt(p, cosmo, 0.);
} else {
/* Slow-cooling regime. */
/* Update du/dt so that we can subsequently drift internal energy. */
const float cooling_du_dt = delta_u / dt_therm;
/* Update the internal energy time derivative */
hydro_set_physical_internal_energy_dt(p, cosmo, cooling_du_dt);
}
if (is_HII_region(xp)) {
/* Check whether final internal energy has reached the HII region
* temperature. */
const double mu_HII = get_mu_HII_region(cooling, XH);
const double u_HII_cgs =
get_u_cgs_from_T(cooling->HIIregion_temp, mu_HII, cooling);
const double u_HII = u_HII_cgs * cooling->internal_energy_from_cgs;
/* The energy has reached, or is below, the HII target energy */
if (u_final <= u_HII) {
/* Immediately set the particle to u_HII */
hydro_set_physical_internal_energy(p, xp, cosmo, u_HII);
hydro_set_drifted_physical_internal_energy(p, cosmo, pressure_floor,
u_HII);
/* Internal energy should stay constant for the timestep */
hydro_set_physical_internal_energy_dt(p, cosmo, 0.f);
}
}
ticks toc_chimes = getticks();
if (clocks_from_ticks(toc_chimes - tic_chimes) > 1000.f) {
warning("Particle %lld took more than 1000ms in CHIMES!", p->id);
chimes_print_gas_vars(stderr, &ChimesGasVars_save, ChimesGlobalVars);
}
/* Store the radiated energy */
xp->cooling_data.radiated_energy -= hydro_get_mass(p) * (u_final - u_0);
/* Free CHIMES memory. */
free_gas_abundances_memory(&ChimesGasVars, ChimesGlobalVars);
free_gas_abundances_memory(&ChimesGasVars_save, ChimesGlobalVars);
/* Set subgrid properties. If the particle is below logT_EOS_max, the CHIMES
* abundances will be re-computed from the subgrid temperature and density. */
cooling_set_subgrid_properties(phys_const, us, cosmo, hydro_properties,
floor_props, cooling, dp, p, xp);
}
/**
* @brief Returns the total radiated energy by this particle.
*
* @param xp The extended particle data
*/
float cooling_get_radiated_energy(const struct xpart *xp) {
return xp->cooling_data.radiated_energy;
}
/**
* @brief Split the cooling content of a particle into n pieces
*
* @param p The #part.
* @param xp The #xpart.
* @param n The number of pieces to split into.
*/
void cooling_split_part(struct part *p, struct xpart *xp, double n) {
xp->cooling_data.radiated_energy /= n;
}
/**
*
* @brief Calculate mean molecular weight from CHIMES abundances and
* also deals with special cases (HII regions, feedback).
*
* @param cooling The #cooling_function_data used in the run.
* @param hydro_props The properties of the hydro scheme.
* @param p Pointer to the particle data.
* @param xp Pointer to the extended particle data.
*/
double chimes_mu(const struct cooling_function_data *cooling,
const struct hydro_props *hydro_props, const struct part *p,
const struct xpart *xp) {
/* Special cases first */
/* Particles heated by FB and we run with heated particles at equilibrium */
if (xp->cooling_data.heated_by_FB && cooling->set_FB_particles_to_eqm) {
return hydro_props->mu_ionised;
}
/* Particles in an HII region */
if (is_HII_region(xp)) {
/* Recover the Hydrogen mass fraction */
const double XH = get_hydrogen_mass_fraction(p, hydro_props);
return get_mu_HII_region(cooling, XH);
}
/* Now, normal case: Compute mu from the particle's metal content and CHIMES
* abundances */
const struct globalVariables *ChimesGlobalVars = &cooling->ChimesGlobalVars;
double numerator = 1.0;
/* Get element mass fractions */
#if defined(CHEMISTRY_COLIBRE) || defined(CHEMISTRY_EAGLE)
const float *const metal_fraction =
chemistry_get_metal_mass_fraction_for_cooling(p);
const float XH = metal_fraction[chemistry_element_H];
numerator += (double)metal_fraction[chemistry_element_He] / XH;
if (ChimesGlobalVars->element_included[0] == 1)
numerator += (double)metal_fraction[chemistry_element_C] / XH;
if (ChimesGlobalVars->element_included[1] == 1)
numerator += (double)metal_fraction[chemistry_element_N] / XH;
if (ChimesGlobalVars->element_included[2] == 1)
numerator += (double)metal_fraction[chemistry_element_O] / XH;
if (ChimesGlobalVars->element_included[3] == 1)
numerator += (double)metal_fraction[chemistry_element_Ne] / XH;
if (ChimesGlobalVars->element_included[4] == 1)
numerator += (double)metal_fraction[chemistry_element_Mg] / XH;
if (ChimesGlobalVars->element_included[5] == 1)
numerator += (double)metal_fraction[chemistry_element_Si] / XH;
if (ChimesGlobalVars->element_included[8] == 1)
numerator += (double)metal_fraction[chemistry_element_Fe] / XH;
if (ChimesGlobalVars->element_included[6] == 1)
numerator +=
(double)metal_fraction[chemistry_element_Si] *
cooling->S_over_Si_ratio_in_solar *
(cooling->S_solar_mass_fraction / cooling->Si_solar_mass_fraction) / XH;
if (ChimesGlobalVars->element_included[7] == 1)
numerator +=
(double)metal_fraction[chemistry_element_Si] *
cooling->Ca_over_Si_ratio_in_solar *
(cooling->Ca_solar_mass_fraction / cooling->Si_solar_mass_fraction) /
XH;
#else
/* Without COLIBRE or EAGLE chemistry, the metal abundances are unavailable.
* -->Set to primordial abundances. */
numerator += 0.0833 * 4.0;
#endif /* CHEMISTRY_COLIBRE || CHEMISTRY_EAGLE */
double denominator = 0.0;
for (int i = 0; i < ChimesGlobalVars->totalNumberOfSpecies; i++)
denominator += (double)xp->cooling_data.chimes_abundances[i];
return numerator / denominator;
}
/**
* @brief Compute the temperature of a #part based on the cooling function.
*
* This uses the mean molecular weight computed from the
* CHIMES abundance array.
*
* @param phys_const #phys_const data structure.
* @param hydro_props The properties of the hydro scheme.
* @param us The internal system of units.
* @param cosmo #cosmology data structure.
* @param cooling #cooling_function_data struct.
* @param p #part data.
* @param xp Pointer to the #xpart data.
*/
float cooling_get_temperature(const struct phys_const *phys_const,
const struct hydro_props *hydro_props,
const struct unit_system *us,
const struct cosmology *cosmo,
const struct cooling_function_data *cooling,
const struct part *p, const struct xpart *xp) {
/* Mean molecular weight */
const double mu = chimes_mu(cooling, hydro_props, p, xp);
/* Internal energy, in code units. */
const double u = hydro_get_physical_internal_energy(p, xp, cosmo);
/* Convert to cgs. */
const double u_cgs = u * cooling->internal_energy_to_cgs;
/* Return particle temperature */
return (float)get_T_from_u_cgs(u_cgs, mu, cooling);
}
/**
* @brief Compute the temperature of a #part on the entropy floor.
*
* This uses the mean molecular weight computed from the
* CHIMES equilibrium abundance tables, not the particle's
* actual abundance array. This is useful for particles
* on the entropy floor.
*
* @param phys_const #phys_const data structure.
* @param hydro_props The properties of the hydro scheme.
* @param us The internal system of units.
* @param cosmo #cosmology data structure.
* @param cooling #cooling_function_data struct.
* @param p #part data.
* @param xp Pointer to the #xpart data.
*/
float cooling_get_temperature_on_entropy_floor(
const struct phys_const *phys_const, const struct hydro_props *hydro_props,
const struct unit_system *us, const struct cosmology *cosmo,
const struct cooling_function_data *cooling, const struct part *p,
const struct xpart *xp) {
/* Internal energy, in code units. */
const double u = hydro_get_physical_internal_energy(p, xp, cosmo);
/* Convert to cgs. */
const double u_cgs = u * cooling->internal_energy_to_cgs;
/* Recover the Hydrogen mass fraction */
const double XH = get_hydrogen_mass_fraction(p, hydro_props);
/* Helium abundance. */
#if defined(CHEMISTRY_COLIBRE) || defined(CHEMISTRY_EAGLE)
const float *const metal_fraction =
chemistry_get_metal_mass_fraction_for_cooling(p);
float metallicity = chemistry_get_total_metal_mass_fraction_for_cooling(p);
if (metallicity <= 0.0) metallicity = FLT_MIN;
const ChimesFloat x_He =
(ChimesFloat)metal_fraction[chemistry_element_He] / (4.0 * XH);
#else
/* Without COLIBRE or EAGLE chemistry, the metal abundances are unavailable.
* --> Set to primordial abundances. */
const float metallicity = FLT_MIN;
const ChimesFloat x_He = 0.0833;
#endif
/* Interpolate equilibrium table to get H2 and electron fractions. */
const ChimesFloat nH = (ChimesFloat)((hydro_get_physical_density(p, cosmo) *
XH / phys_const->const_proton_mass) *
cooling->number_density_to_cgs);
const ChimesFloat Z_over_Zsol = (ChimesFloat)metallicity / cooling->Zsol;
/* Initial guess for T, assuming mu = 1 */
const ChimesFloat T = get_T_from_u_cgs(u_cgs, /*mu=*/1.0, cooling);
/* Take log's, but protect against log(0) */
const ChimesFloat log_nH = log10(chimes_max(nH, FLT_MIN));
const ChimesFloat log_Z_over_Zsol = log10(chimes_max(Z_over_Zsol, FLT_MIN));
const ChimesFloat log_T = log10(chimes_max(T, FLT_MIN));
const ChimesFloat log_u = log10(chimes_max(u_cgs, FLT_MIN));
int T_index, nH_index, Z_index;
ChimesFloat dT, dnH, dZ;
const int N_T = chimes_table_eqm_abundances.N_Temperatures;
const int N_nH = chimes_table_eqm_abundances.N_Densities;
const int N_Z = chimes_table_eqm_abundances.N_Metallicities;
chimes_get_table_index(chimes_table_eqm_abundances.Temperatures, N_T, log_T,
&T_index, &dT);
chimes_get_table_index(chimes_table_eqm_abundances.Densities, N_nH, log_nH,
&nH_index, &dnH);
chimes_get_table_index(chimes_table_eqm_abundances.Metallicities, N_Z,
log_Z_over_Zsol, &Z_index, &dZ);
/* Position of electrons and H2 in the CHIMES abundance arrays. */
const int elec_ind = cooling->ChimesGlobalVars.speciesIndices[sp_elec];
const int H2_ind = cooling->ChimesGlobalVars.speciesIndices[sp_H2];
/* Find the temperature bins in the eqm tables that bracket u_cgs */
dT = 0.0;
ChimesFloat mu, x_H2, x_elec, u_low, u_hi, T_low, T_hi;
x_elec = exp10(chimes_interpol_4d_fix_x(
chimes_table_eqm_abundances.Abundances, elec_ind, T_index, nH_index,
Z_index, dT, dnH, dZ, N_T, N_nH, N_Z));
x_H2 = exp10(chimes_interpol_4d_fix_x(chimes_table_eqm_abundances.Abundances,
H2_ind, T_index, nH_index, Z_index, dT,
dnH, dZ, N_T, N_nH, N_Z));
mu = (1.0 + (4.0 * x_He)) / (1.0 + x_He + x_elec - x_H2);
T_low = exp10(chimes_table_eqm_abundances.Temperatures[T_index]);
u_low = get_u_cgs_from_T(T_low, mu, cooling);
x_elec = exp10(chimes_interpol_4d_fix_x(
chimes_table_eqm_abundances.Abundances, elec_ind, T_index + 1, nH_index,
Z_index, dT, dnH, dZ, N_T, N_nH, N_Z));
x_H2 = exp10(chimes_interpol_4d_fix_x(chimes_table_eqm_abundances.Abundances,
H2_ind, T_index + 1, nH_index, Z_index,
dT, dnH, dZ, N_T, N_nH, N_Z));
mu = (1.0 + (4.0 * x_He)) / (1.0 + x_He + x_elec - x_H2);
T_hi = exp10(chimes_table_eqm_abundances.Temperatures[T_index + 1]);
u_hi = get_u_cgs_from_T(T_hi, mu, cooling);
if (u_cgs < u_low) {
while (u_cgs < u_low) {
T_index -= 1;
if (T_index < 0)
error(
"CHIMES failed to convert internal energy to temperature for "
"snapshot, as T_index < 0.");
x_elec = exp10(chimes_interpol_4d_fix_x(
chimes_table_eqm_abundances.Abundances, elec_ind, T_index, nH_index,
Z_index, dT, dnH, dZ, N_T, N_nH, N_Z));
x_H2 = exp10(chimes_interpol_4d_fix_x(
chimes_table_eqm_abundances.Abundances, H2_ind, T_index, nH_index,
Z_index, dT, dnH, dZ, N_T, N_nH, N_Z));
mu = (1.0 + (4.0 * x_He)) / (1.0 + x_He + x_elec - x_H2);
T_low = exp10(chimes_table_eqm_abundances.Temperatures[T_index]);
u_low = get_u_cgs_from_T(T_low, mu, cooling);
x_elec = exp10(chimes_interpol_4d_fix_x(
chimes_table_eqm_abundances.Abundances, elec_ind, T_index + 1,
nH_index, Z_index, dT, dnH, dZ, N_T, N_nH, N_Z));
x_H2 = exp10(chimes_interpol_4d_fix_x(
chimes_table_eqm_abundances.Abundances, H2_ind, T_index + 1, nH_index,
Z_index, dT, dnH, dZ, N_T, N_nH, N_Z));
mu = (1.0 + (4.0 * x_He)) / (1.0 + x_He + x_elec - x_H2);
T_hi = exp10(chimes_table_eqm_abundances.Temperatures[T_index + 1]);
u_hi = get_u_cgs_from_T(T_hi, mu, cooling);
}
} else if (u_cgs > u_hi) {
while (u_cgs > u_hi) {
T_index += 1;
if (T_index > N_T - 2)
error(
"CHIMES failed to convert internal energy to temperature for "
"snapshot, as T_index > N_T - 2.");
x_elec = exp10(chimes_interpol_4d_fix_x(
chimes_table_eqm_abundances.Abundances, elec_ind, T_index, nH_index,
Z_index, dT, dnH, dZ, N_T, N_nH, N_Z));
x_H2 = exp10(chimes_interpol_4d_fix_x(
chimes_table_eqm_abundances.Abundances, H2_ind, T_index, nH_index,
Z_index, dT, dnH, dZ, N_T, N_nH, N_Z));
mu = (1.0 + (4.0 * x_He)) / (1.0 + x_He + x_elec - x_H2);
T_low = exp10(chimes_table_eqm_abundances.Temperatures[T_index]);
u_low = get_u_cgs_from_T(T_low, mu, cooling);
x_elec = exp10(chimes_interpol_4d_fix_x(
chimes_table_eqm_abundances.Abundances, elec_ind, T_index + 1,
nH_index, Z_index, dT, dnH, dZ, N_T, N_nH, N_Z));
x_H2 = exp10(chimes_interpol_4d_fix_x(
chimes_table_eqm_abundances.Abundances, H2_ind, T_index + 1, nH_index,
Z_index, dT, dnH, dZ, N_T, N_nH, N_Z));
mu = (1.0 + (4.0 * x_He)) / (1.0 + x_He + x_elec - x_H2);
T_hi = exp10(chimes_table_eqm_abundances.Temperatures[T_index + 1]);
u_hi = get_u_cgs_from_T(T_hi, mu, cooling);
}
}
/* Finally, interpolate mu between
* u_low and u_hi. */
if (!((u_cgs >= u_low) && (u_cgs <= u_hi)))
error(
"CHIMES failed to convert internal energy to temperature for snapshot, "
"as u_low and u_hi do not bracket u_cgs.");
dT = (log_u - log10(u_low)) / (log10(u_hi) - log10(u_low));
x_elec = exp10(chimes_interpol_4d_fix_x(
chimes_table_eqm_abundances.Abundances, elec_ind, T_index, nH_index,
Z_index, dT, dnH, dZ, N_T, N_nH, N_Z));
x_H2 = exp10(chimes_interpol_4d_fix_x(chimes_table_eqm_abundances.Abundances,
H2_ind, T_index, nH_index, Z_index, dT,
dnH, dZ, N_T, N_nH, N_Z));
mu = (1.0 + (4.0 * x_He)) / (1.0 + x_He + x_elec - x_H2);
/* Return particle temperature */
return (float)get_T_from_u_cgs(u_cgs, mu, cooling);
}
/**
* @brief Compute the temperature of a #part for the snapshot.
*
* This function returns the temperature that will be written
* out to the snapshot. If it is above the entropy floor,
* it uses the mean molecular weight calculated from the
* particle's abundance array. On the entropy floor, it
* uses the equilibrium abundance tables.
*
* @param phys_const #phys_const data structure.
* @param hydro_props The properties of the hydro scheme.
* @param us The internal system of units.
* @param cosmo #cosmology data structure.
* @param cooling #cooling_function_data struct.
* @param p #part data.
* @param xp Pointer to the #xpart data.
*/
float cooling_get_snapshot_temperature(
const struct phys_const *phys_const, const struct hydro_props *hydro_props,
const struct entropy_floor_properties *floor_props,
const struct unit_system *us, const struct cosmology *cosmo,
const struct cooling_function_data *cooling, const struct part *p,
const struct xpart *xp) {
/* Limit imposed by the entropy floor */
const double A_EOS = entropy_floor(p, cosmo, floor_props);
const double rho = hydro_get_physical_density(p, cosmo);
const double u_EOS = gas_internal_energy_from_entropy(rho, A_EOS);
const double u_EOS_max = u_EOS * exp10(cooling->dlogT_EOS);
/* Particle's internal energy */
const double u = hydro_get_physical_internal_energy(p, xp, cosmo);
if ((u < u_EOS_max) && !is_HII_region(xp)) {
/* If the particle is within dlogT_EOS of the entropy floor, and is not an
* HII region, its CHIMES abundance array would have been set according to
* the subgrid properties.
* We therefore need to calculate the temperature based on the equilibrium
* abundance tables at the particle's actual temperature, not subgrid.
*/
return cooling_get_temperature_on_entropy_floor(phys_const, hydro_props, us,
cosmo, cooling, p, xp);
} else {
/* Above u_EOS_max, we use the particle's actual abundance array; */
return cooling_get_temperature(phys_const, hydro_props, us, cosmo, cooling,
p, xp);
}
}
/**
* @brief Compute the electron number density of a #part based on the cooling
* function.
*
* The electron density returned is in physical internal units.
*
* @param phys_const #phys_const data structure.
* @param hydro_props The properties of the hydro scheme.
* @param floor_props The properties of the entropy floor.
* @param us The internal system of units.
* @param cosmo #cosmology data structure.
* @param cooling #cooling_function_data struct.
* @param dust_props The properties of the dust model.
* @param p #part data.
* @param xp Pointer to the #xpart data.
*/
double cooling_get_electron_density(
const struct phys_const *phys_const, const struct hydro_props *hydro_props,
const struct entropy_floor_properties *floor_props,
const struct unit_system *us, const struct cosmology *cosmo,
const struct cooling_function_data *cooling,
const struct dustevo_props *dust_props, const struct part *p,
const struct xpart *xp) {
const struct colibre_cooling_tables *table = &(cooling->colibre_table);
const struct globalVariables *ChimesGlobalVars = &(cooling->ChimesGlobalVars);
/* Create ChimesGasVars struct */
struct gasVariables ChimesGasVars;
allocate_gas_abundances_memory(&ChimesGasVars, ChimesGlobalVars);
/* Copy abundances over from xp to ChimesGasVars */
for (int i = 0; i < ChimesGlobalVars->totalNumberOfSpecies; i++)
ChimesGasVars.abundances[i] = xp->cooling_data.chimes_abundances[i];
/* Particle's internal energy */
const double u = hydro_get_physical_internal_energy(p, xp, cosmo);
const double u_cgs = u * cooling->internal_energy_to_cgs;
/* Update element abundances */
chimes_update_element_abundances(phys_const, us, cosmo, cooling, hydro_props,
p, xp, &ChimesGasVars, u,
/*mode=*/1);
/* Special treatment for particles recently heated by feedback */
if (xp->cooling_data.heated_by_FB && cooling->set_FB_particles_to_eqm) {
/* Update the ChimesGasVars with the particle's thermo variables. */
chimes_update_gas_vars(u_cgs, phys_const, us, cosmo, hydro_props,
floor_props, cooling, dust_props, p, xp,
&ChimesGasVars, 1.0);
/* Set the particle abundances to equilibrium values */
cooling_set_FB_particle_chimes_abundances(&ChimesGasVars, cooling);
}
/* Special treatment for particles that are part of HII regions */
else if (is_HII_region(xp)) {
/* Set HII region abundance array */
cooling_set_HIIregion_chimes_abundances(&ChimesGasVars, cooling);
}
/* Recover the Hydrogen mass fraction */
const double XH = get_hydrogen_mass_fraction(p, hydro_props);
/* Calculate hydrogen number density nH
* in internal units. */
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);
const double T_cgs = cooling_get_temperature(phys_const, hydro_props, us,
cosmo, cooling, p, xp) *
units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
/* Calculate electron density in internal units
* using CHIMES abundance array. */
const int elec_ind = ChimesGlobalVars->speciesIndices[sp_elec];
const double n_e_non_eq = n_H * ChimesGasVars.abundances[elec_ind];
/* Free CHIMES memory. */
free_gas_abundances_memory(&ChimesGasVars, ChimesGlobalVars);
double n_e_metal;
if (ChimesGlobalVars->hybrid_cooling_mode == 1) {
/* When we use hybrid cooling, we need to add
* the contribution of metals that are not included
* in the non-eq network, which we read in from
* the equilibrium cooling tables. */
/* Get the total metallicity in units of solar */
float abundance_ratio[colibre_cooling_N_elementtypes];
const float logZZsol =
abundance_ratio_to_solar(p, &(cooling->colibre_table), abundance_ratio);
// Set weights for electron densities.
float weights_electron[colibre_cooling_N_electrontypes];
for (int i = 0; i < colibre_cooling_N_electrontypes; i++) {
if (i < colibre_cooling_N_elementtypes - 1)
weights_electron[i] = abundance_ratio[i];
else
weights_electron[i] = 1.f;
}
/* Get indices of T, nH, and metallicity*/
int T_index, n_H_index, met_index, red_index;
float d_T, d_n_H, d_met;
cooling_get_index_1d(table->Temp, colibre_cooling_N_temperature,
log10(T_cgs), &T_index, &d_T);
cooling_get_index_1d(table->Metallicity, colibre_cooling_N_metallicity,
logZZsol, &met_index, &d_met);
cooling_get_index_1d(table->nH, colibre_cooling_N_density, log10(n_H_cgs),
&n_H_index, &d_n_H);
if (table->dz > 0.5)
red_index = 1;
else
red_index = 0;
if (met_index == 0) {
/* Primordial abundances; metal
* contribution is zero. */
n_e_metal = 0.0;
} else {
/* From metals, with actual metal ratios.
* We also need to exclude any metals
* that are already included in CHIMES. */
for (int i = element_C; i < element_OA; i++) {
if (ChimesGlobalVars->element_included[i - 2] == 1)
weights_electron[i] = 0.f;
}
n_e_metal =
n_H * interpolation3d_fix_x_plus_summation(
table->Telectron_fraction, weights_electron, element_C,
colibre_cooling_N_electrontypes - 4, red_index, T_index,
met_index, n_H_index, d_T, d_met, d_n_H,
colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_electrontypes);
}
} else {
/* If we are not using hybrid cooling, we
* only include the electrons that are in
* the non-eq network. */
n_e_metal = 0.0;
}
return n_e_non_eq + n_e_metal;
}
/**
* @brief Compute the electron pressure of a #part based on the cooling
* function.
*
* Returns the total electron pressure of the particle, P_e * V, in code units.
*
* @param phys_const #phys_const data structure.
* @param hydro_props The properties of the hydro scheme.
* @param floor_props The properties of the entropy floor.
* @param us The internal system of units.
* @param cosmo #cosmology data structure.
* @param cooling #cooling_function_data struct.
* @param dust_props The properties of the dust model.
* @param p #part data.
* @param xp Pointer to the #xpart data.
*/
double cooling_get_electron_pressure(
const struct phys_const *phys_const, const struct hydro_props *hydro_props,
const struct entropy_floor_properties *floor_props,
const struct unit_system *us, const struct cosmology *cosmo,
const struct cooling_function_data *cooling,
const struct dustevo_props *dust_props, const struct part *p,
const struct xpart *xp) {
/* Get internal energy in physical frame */
const float u_phys = hydro_get_physical_internal_energy(p, xp, cosmo);
/* Calculate temperature in internal units. */
const double mu = chimes_mu(cooling, hydro_props, p, xp);
const double temperature = get_T_from_u_internal(u_phys, mu, phys_const);
/* Calculate total electron density in internal units using
* CHIMES abundance array and elements evolved in equilibrium. */
const double n_e =
cooling_get_electron_density(phys_const, hydro_props, floor_props, us,
cosmo, cooling, dust_props, p, xp);
return n_e * phys_const->const_boltzmann_k * temperature;
}
/**
* @brief Compute the y-Compton contribution of a #part based on the cooling
* function.
*
* This is eq. (2) of McCarthy et al. (2017).
*
* @param phys_const #phys_const data structure.
* @param hydro_props The properties of the hydro scheme.
* @param floor_props The properties of the entropy floor.
* @param us The internal system of units.
* @param cosmo #cosmology data structure.
* @param cooling #cooling_function_data struct.
* @param dust_props The properties of the dust model.
* @param p #part data.
* @param xp Pointer to the #xpart data.
*/
double cooling_get_ycompton(const struct phys_const *phys_const,
const struct hydro_props *hydro_props,
const struct entropy_floor_properties *floor_props,
const struct unit_system *us,
const struct cosmology *cosmo,
const struct cooling_function_data *cooling,
const struct dustevo_props *dust_props,
const struct part *p, const struct xpart *xp) {
/* Get quantities in physical frame */
const float u_phys = hydro_get_physical_internal_energy(p, xp, cosmo);
const float rho_phys = hydro_get_physical_density(p, cosmo);
const double m = hydro_get_mass(p);
/* Calculate temperature in internal units. */
const double mu = chimes_mu(cooling, hydro_props, p, xp);
const double temperature = get_T_from_u_internal(u_phys, mu, phys_const);
/* Calculate total electron density in internal units using
* CHIMES abundance array and elements evolved in equilibrium. */
const double n_e =
cooling_get_electron_density(phys_const, hydro_props, floor_props, us,
cosmo, cooling, dust_props, p, xp);
return cooling->y_compton_factor * temperature * m * n_e / rho_phys;
}
/**
* @brief Calculate the N_ref column density.
*
* This routine returns the column density N_ref
* as defined in Ploeckinger et al. (in prep),
* which is used to scale the ISRF, cosmic rays,
* dust depletion and shielding column density
* in COLIBRE.
*
* @param phys_const #phys_const data structure.
* @param us The internal system of units.
* @param cosmo #cosmology data structure.
* @param cooling #cooling_function_data struct.
* @param p #part data.
* @param xp Pointer to the #xpart data.
* @param temperature The gas temperature.
* @param mu Mean molecular weight.
*/
double calculate_colibre_N_ref(const struct phys_const *phys_const,
const struct unit_system *us,
const struct cosmology *cosmo,
const struct cooling_function_data *cooling,
const struct hydro_props *hydro_props,
const struct part *p, const struct xpart *xp,
const double temperature) {
/* Based on Ploeckinger & Schaye (2020). */
const float N_min = cooling->nref_column_density_min_cgs;
const float N_max = cooling->nref_column_density_max_cgs;
const float T_min = cooling->nref_trans_temperature_min_K;
const float T_max = cooling->nref_trans_temperature_max_K;
const double l_max_cgs = cooling->nref_max_length_cgs;
const float k = cooling->nref_trans_steepness_k;
const double XH_const = cooling->nref_XH_const;
const double mu_const = cooling->nref_mu_const;
/* Recover the Hydrogen mass fraction */
const double XH = get_hydrogen_mass_fraction(p, hydro_props);
/* Density */
const double nH =
hydro_get_physical_density(p, cosmo) * XH / phys_const->const_proton_mass;
const double nH_cgs = nH * cooling->number_density_to_cgs;
/* Jeans column density */
const double N_J_thermal =
nH_cgs * sqrt(hydro_gamma * cooling->boltzmann_k_cgs * temperature /
(mu_const * cooling->newton_G_cgs *
(cooling->proton_mass_cgs * nH_cgs / XH_const) *
cooling->proton_mass_cgs));
double N_J = N_J_thermal;
/* Special case where we use the turbulent Jeans scale */
if (cooling->colibre_use_turbulent_jeans_length) {
const double N_J_turb =
sqrt(hydro_gamma * XH_const * nH_cgs /
(cooling->proton_mass_cgs * cooling->newton_G_cgs)) *
(cooling->sigma_turb_km_p_s * 1.0e5);
N_J = chimes_max(N_J_thermal, N_J_turb);
}
double N_ref_prime = chimes_min(N_J, N_max);
if (l_max_cgs > 0.0)
N_ref_prime = chimes_min(N_ref_prime, l_max_cgs * nH_cgs);
/* Apply the sigmoid */
double N_ref = exp10(log10(N_ref_prime) -
(log10(N_ref_prime) - log10(N_min)) /
(1.0 + pow(sqrt(T_min * T_max) / temperature, k)));
return N_ref;
}
/**
* @brief Write a cooling struct to the given FILE as a stream of bytes.
*
* @param cooling the struct
* @param stream the file stream
*/
void cooling_struct_dump(const struct cooling_function_data *cooling,
struct dustevo_props *dp, FILE *stream) {
/* Zero all pointers in the colibre_table within
* the cooling_function_data struct. */
struct cooling_function_data cooling_copy = *cooling;
cooling_copy.colibre_table.Tcooling = NULL;
cooling_copy.colibre_table.Tcooling_reduced = NULL;
cooling_copy.colibre_table.Theating = NULL;
cooling_copy.colibre_table.Theating_reduced = NULL;
cooling_copy.colibre_table.Telectron_fraction = NULL;
cooling_copy.colibre_table.Redshifts = NULL;
cooling_copy.colibre_table.nH = NULL;
cooling_copy.colibre_table.Temp = NULL;
cooling_copy.colibre_table.Metallicity = NULL;
cooling_copy.colibre_table.LogAbundances = NULL;
cooling_copy.colibre_table.Abundances = NULL;
cooling_copy.colibre_table.Abundances_inv = NULL;
cooling_copy.colibre_table.atomicmass = NULL;
cooling_copy.colibre_table.atomicmass_inv = NULL;
cooling_copy.colibre_table.LogMassFractions = NULL;
cooling_copy.colibre_table.MassFractions = NULL;
cooling_copy.colibre_table.Zsol = NULL;
cooling_copy.colibre_table.Zsol_inv = NULL;
cooling_copy.ChimesGlobalVars.hybrid_data = NULL;
restart_write_blocks((void *)&cooling_copy,
sizeof(struct cooling_function_data), 1, stream,
"cooling", "cooling function");
}
/**
* @brief Restore a hydro_props struct from the given FILE as a stream of
* bytes.
*
* @param cooling the struct
* @param stream the file stream
* @param cosmo #cosmology structure
*/
void cooling_struct_restore(struct cooling_function_data *cooling,
struct dustevo_props *dp, FILE *stream,
const struct cosmology *cosmo) {
restart_read_blocks((void *)cooling, sizeof(struct cooling_function_data), 1,
stream, NULL, "cooling function");
/* Initialise the CHIMES module. */
if (engine_rank == 0) message("Initialising CHIMES cooling module.");
init_chimes(&cooling->ChimesGlobalVars);
if ((cooling->ChimesGlobalVars.hybrid_cooling_mode == 1) ||
(cooling->use_colibre_subgrid_EOS == 1)) {
/* Read the Colibre table. */
if (engine_rank == 0) message("Reading Colibre cooling table.");
get_cooling_redshifts(&(cooling->colibre_table));
char fname[colibre_table_path_name_length + colibre_table_filebase_length +
12];
sprintf(fname, "%s/%s_z0.00.hdf5",
cooling->colibre_table.cooling_table_path,
cooling->colibre_table.cooling_table_filebase);
read_cooling_header(fname, &(cooling->colibre_table));
/* Allocate memory for the tables */
allocate_cooling_tables(&(cooling->colibre_table));
/* Force a re-read of the cooling tables */
cooling->colibre_table.z_index = -10;
cooling->colibre_table.previous_z_index =
cooling->colibre_table.number_of_redshifts - 2;
cooling_update(cosmo, /*pfloor=*/NULL, cooling, /*space=*/NULL, dp,
/*time=*/0);
}
/* Check that the size of the CHIMES network
* determined by the various Include
* parameters matches the size specified
* through the configure script. */
if (cooling->ChimesGlobalVars.totalNumberOfSpecies != CHIMES_NETWORK_SIZE)
error(
"The size of the CHIMES network based on the Include "
"parameters is %d. However, Swift was configured and compiled using "
"--with-chimes-network-size=%d. If you are sure that you have "
"correctly set which metals to include in the network in the parameter "
"file, then you will need to re-compile Swift using ./configure "
"--with-cooling=CHIMES --with-chimes-network-size=%d.",
cooling->ChimesGlobalVars.totalNumberOfSpecies, CHIMES_NETWORK_SIZE,
cooling->ChimesGlobalVars.totalNumberOfSpecies);
if (cooling->ChimesGlobalVars.hybrid_cooling_mode == 1) {
/* Create data structure for hybrid cooling,
* and store pointer to the Colibre table. */
cooling->ChimesGlobalVars.hybrid_data =
(void *)malloc(sizeof(struct global_hybrid_data_struct));
struct global_hybrid_data_struct *myData;
myData = (struct global_hybrid_data_struct *)
cooling->ChimesGlobalVars.hybrid_data;
myData->table = &(cooling->colibre_table);
myData->Zsol = cooling->Zsol;
/* Set the hybrid cooling function pointers. */
cooling->ChimesGlobalVars.hybrid_cooling_fn =
&colibre_metal_cooling_rate_temperature;
cooling->ChimesGlobalVars.allocate_gas_hybrid_data_fn =
&chimes_allocate_gas_hybrid_data;
cooling->ChimesGlobalVars.free_gas_hybrid_data_fn =
&chimes_free_gas_hybrid_data;
}
}
/**
* @brief Converts cooling quantities of a particle at the start of a run
*
* This function is called once at the end of the engine_init_particle()
* routine (at the start of a calculation) after the densities of
* particles have been computed.
*
* For this cooling module, this routine is used to set the cooling
* properties of the (x-)particles to a valid start state, in particular
* the CHIMES abundance array.
*
* This is controlled by the cooling->init_abundance_mode as follows:
* 0 -- Set each element to one ionisation state, determined by the
* ChimesGlobalVars.InitIonState parameter.
* 1 -- Read abundances from eqm abundance tables.
* 2 -- Compute initial equilibrium abundances.
*
* @param p The particle to act upon
* @param xp The extended particle to act upon
* @param cosmo The cosmological model.
* @param hydro_props The constants used in the scheme.
* @param phys_const #phys_const data structure.
* @param us Internal system of units data structure.
* @param floor_props Properties of the entropy floor.
* @param cooling #cooling_function_data data structure.
*/
void cooling_convert_quantities(
struct part *p, struct xpart *xp, const struct cosmology *cosmo,
const struct hydro_props *hydro_props, const struct phys_const *phys_const,
const struct unit_system *us,
const struct entropy_floor_properties *floor_props,
const struct cooling_function_data *cooling,
const struct dustevo_props *dp) {
const struct globalVariables *ChimesGlobalVars = &cooling->ChimesGlobalVars;
/* Allocate memory to arrays within ChimesGasVars. */
struct gasVariables ChimesGasVars;
allocate_gas_abundances_memory(&ChimesGasVars, ChimesGlobalVars);
/* Set element abundances from
* metal mass fractions. */
const double u_0 = hydro_get_physical_internal_energy(p, xp, cosmo);
chimes_update_element_abundances(phys_const, us, cosmo, cooling, hydro_props,
p, xp, &ChimesGasVars, u_0, /*mode=*/0);
/* Zero particle's radiated energy. */
xp->cooling_data.radiated_energy = 0.f;
/* Set initial values for CHIMES abundance array. */
ChimesGasVars.InitIonState = cooling->InitIonState;
initialise_gas_abundances(&ChimesGasVars, ChimesGlobalVars);
if (cooling->init_abundance_choice == init_abundance_single) {
/* Copy abundances over to xp. */
for (int i = 0; i < ChimesGlobalVars->totalNumberOfSpecies; i++)
xp->cooling_data.chimes_abundances[i] = ChimesGasVars.abundances[i];
} else if ((cooling->init_abundance_choice == init_abundance_read) ||
(cooling->init_abundance_choice == init_abundance_compute)) {
/* Copy initial abundances over to xp. */
for (int i = 0; i < ChimesGlobalVars->totalNumberOfSpecies; i++)
xp->cooling_data.chimes_abundances[i] = ChimesGasVars.abundances[i];
/* Convert internal energy to cgs */
const double u_0_cgs = u_0 * cooling->internal_energy_to_cgs;
/* If computing the eqm (init_abundance_compute), we will integrate the
* chemistry ten times for 1 Gyr per iteration. Multiple iterations are
* required so that the shielding column densities can be updated between
* each iteration. If reading from tables, we only need 1 iteration. */
int n_iterations;
if (cooling->init_abundance_choice == init_abundance_read)
n_iterations = 1;
else /* cooling->init_abundance_choice == init_abundance_compute */
n_iterations = 10;
/* 1 Gyr in cgs */
const double dt_cgs = 1e9 * 365.25 * 24. * 3600.;
for (int i = 0; i < n_iterations; i++) {
/* Update element abundances. This
* accounts for the dust depletion
* factors. */
chimes_update_element_abundances(phys_const, us, cosmo, cooling,
hydro_props, p, xp, &ChimesGasVars, u_0,
/*mode=*/1);
/* Update ChimesGasVars with the particle's
* thermodynamic variables. */
chimes_update_gas_vars(u_0_cgs, phys_const, us, cosmo, hydro_props,
floor_props, cooling, dp, p, xp, &ChimesGasVars,
dt_cgs);
/* Set temperature evolution off, so that we
* compute equilibrium at fixed temperature. */
ChimesGasVars.ThermEvolOn = 0;
/* Determine whether to use eqm tables or compute the eqm state. */
if (cooling->init_abundance_choice == init_abundance_read)
ChimesGasVars.ForceEqOn = 1;
else
ChimesGasVars.ForceEqOn = 0;
/* Integrate to chemical equilibrium. */
const int chimes_return_code =
chimes_network(&ChimesGasVars, ChimesGlobalVars);
/* Check whether the abundances became
* unphysical before enforcing the
* constraint equations within CHIMES. */
if (chimes_return_code == 1) {
chimes_print_gas_vars(stderr, &ChimesGasVars, ChimesGlobalVars);
error(
"CHIMES ERROR: When computing the initial equilibrium abundances, "
"chimes_network() returned an error code. This suggests we might "
"need to use stricter tolerances in CHIMES.");
}
}
/* Calculate net cooling rate. */
const float rho_phys_cgs =
hydro_get_physical_density(p, cosmo) * cooling->density_to_cgs;
/* cooling rate in erg/cm^3/s. */
float net_cooling_rate_cgs = (float)chimes_cooling_rate_external(
&ChimesGasVars, ChimesGlobalVars, 0);
/* Convert to erg/g/s */
net_cooling_rate_cgs /= rho_phys_cgs;
/* Convert to code units. */
xp->cooling_data.net_cooling_rate =
net_cooling_rate_cgs * cooling->internal_energy_over_time_from_cgs;
/* Copy final abundances over to xp. */
for (int i = 0; i < ChimesGlobalVars->totalNumberOfSpecies; i++)
xp->cooling_data.chimes_abundances[i] = ChimesGasVars.abundances[i];
}
/* Free CHIMES memory. */
free_gas_abundances_memory(&ChimesGasVars, ChimesGlobalVars);
}
/**
* @brief Set the subgrid properties of the gas 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 hydro_props struct
* @param floor_props Properties of the entropy floor.
* @param cooling The #cooling_function_data used in the run.
* @param p Pointer to the particle data.
* @param xp Pointer to the extended particle data.
*/
void cooling_set_subgrid_properties(
const struct phys_const *phys_const, const struct unit_system *us,
const struct cosmology *cosmo, const struct hydro_props *hydro_props,
const struct entropy_floor_properties *floor_props,
const struct cooling_function_data *cooling, const struct dustevo_props *dp,
struct part *p, struct xpart *xp) {
/* Limit imposed by the entropy floor */
const double A_EOS = entropy_floor(p, cosmo, floor_props);
const double rho = hydro_get_physical_density(p, cosmo);
const double u_EOS = gas_internal_energy_from_entropy(rho, A_EOS);
const double u_EOS_max = u_EOS * exp10(cooling->dlogT_EOS);
const float log10_u_EOS_max_cgs = log10f(
u_EOS_max * cooling->colibre_table.internal_energy_to_cgs + FLT_MIN);
/* Particle's internal energy */
const double u = hydro_get_physical_internal_energy(p, xp, cosmo);
/* Particle's temperature */
const float T = cooling_get_snapshot_temperature(
phys_const, hydro_props, floor_props, us, cosmo, cooling, p, xp);
if (cooling->use_colibre_subgrid_EOS == 1) {
const float log10_T = log10f(T);
/* Recover the Hydrogen mass fraction */
const float XH = get_hydrogen_mass_fraction(p, hydro_props);
/* Get the particle pressure */
const float P_phys = hydro_get_physical_pressure(p, cosmo);
/* Are we in an HII region? */
const int HII_region = xp->tracers_data.HIIregion_timer_gas > 0.;
/* Get the total metallicity in units of solar */
float abundance_ratio[colibre_cooling_N_elementtypes];
const float logZZsol =
abundance_ratio_to_solar(p, &(cooling->colibre_table), abundance_ratio);
/* Compute the subgrid properties */
p->cooling_data.subgrid_temp = compute_subgrid_property(
cooling, phys_const, floor_props, cosmo, rho, logZZsol, XH, P_phys,
log10_T, log10_u_EOS_max_cgs, HII_region, abundance_ratio,
log10(u * cooling->colibre_table.internal_energy_to_cgs),
colibre_compute_subgrid_temperature);
p->cooling_data.subgrid_dens = compute_subgrid_property(
cooling, phys_const, floor_props, cosmo, rho, logZZsol, XH, P_phys,
log10_T, log10_u_EOS_max_cgs, HII_region, abundance_ratio,
log10(u * cooling->colibre_table.internal_energy_to_cgs),
colibre_compute_subgrid_density);
} else {
p->cooling_data.subgrid_temp = T;
p->cooling_data.subgrid_dens = rho;
}
if (xp->tracers_data.HIIregion_timer_gas > 0.) {
/* If the particle is flagged as an HII region,
* we need to set its abundance array accordingly. */
const struct globalVariables *ChimesGlobalVars = &cooling->ChimesGlobalVars;
/* Create ChimesGasVars struct */
struct gasVariables ChimesGasVars;
allocate_gas_abundances_memory(&ChimesGasVars, ChimesGlobalVars);
/* Copy abundances over from xp to ChimesGasVars */
for (int i = 0; i < ChimesGlobalVars->totalNumberOfSpecies; i++)
ChimesGasVars.abundances[i] = xp->cooling_data.chimes_abundances[i];
/* Update element abundances */
chimes_update_element_abundances(phys_const, us, cosmo, cooling,
hydro_props, p, xp, &ChimesGasVars, u,
/*mode=*/1);
/* Set HII region abundance array */
cooling_set_HIIregion_chimes_abundances(&ChimesGasVars, cooling);
/* Copy abundances from ChimesGasVars to xp. */
for (int i = 0; i < ChimesGlobalVars->totalNumberOfSpecies; i++)
xp->cooling_data.chimes_abundances[i] = ChimesGasVars.abundances[i];
/* Free CHIMES memory. */
free_gas_abundances_memory(&ChimesGasVars, ChimesGlobalVars);
} else if (u < u_EOS_max) {
/* If the particle is within dlogT_EOS of the entropy floor, we need to
* re-set its abundance array to be in equilibriu using the subgrid density
* and temperature.
*
* Note that, if use_colibe_subgrid_EOS == 0, the subgrid density and
* temperature have
* simply been set to the particle density and temperature. */
const struct globalVariables *ChimesGlobalVars = &cooling->ChimesGlobalVars;
struct gasVariables ChimesGasVars;
allocate_gas_abundances_memory(&ChimesGasVars, ChimesGlobalVars);
for (int i = 0; i < ChimesGlobalVars->totalNumberOfSpecies; i++)
ChimesGasVars.abundances[i] = xp->cooling_data.chimes_abundances[i];
chimes_update_element_abundances(phys_const, us, cosmo, cooling,
hydro_props, p, xp, &ChimesGasVars, u,
/*mode=*/1);
const double mu = chimes_mu(cooling, hydro_props, p, xp);
const double T_subgrid = p->cooling_data.subgrid_temp;
const double u_subgrid = get_u_internal_from_T(T_subgrid, mu, phys_const);
const double u_subgrid_cgs = u_subgrid * cooling->internal_energy_to_cgs;
chimes_update_gas_vars(u_subgrid_cgs, phys_const, us, cosmo, hydro_props,
floor_props, cooling, dp, p, xp, &ChimesGasVars,
1.0);
/* Recover the Hydrogen mass fraction */
const ChimesFloat XH =
(ChimesFloat)get_hydrogen_mass_fraction(p, hydro_props);
ChimesGasVars.nH_tot = XH * p->cooling_data.subgrid_dens *
cooling->number_density_to_cgs /
phys_const->const_proton_mass;
ChimesGasVars.temperature = p->cooling_data.subgrid_temp;
ChimesGasVars.ForceEqOn = 1;
ChimesGasVars.ThermEvolOn = 0;
chimes_network(&ChimesGasVars, ChimesGlobalVars);
/* Copy abundances from ChimesGasVars back to xp. */
for (int i = 0; i < ChimesGlobalVars->totalNumberOfSpecies; i++)
xp->cooling_data.chimes_abundances[i] = ChimesGasVars.abundances[i];
/* Set the net cooling rate to zero */
xp->cooling_data.net_cooling_rate = 0.0f;
/* Free CHIMES memory. */
free_gas_abundances_memory(&ChimesGasVars, ChimesGlobalVars);
} else if (xp->cooling_data.heated_by_FB &&
(cooling->set_FB_particles_to_eqm ||
cooling->destroy_FB_heated_dust)) {
/* Prior to computing updated element abundances, destroy dust. */
if (cooling->destroy_FB_heated_dust) {
dust_redistribute_masses_cooling_step(p, dp);
}
if (cooling->set_FB_particles_to_eqm) {
/* If the particle is not on, or near, the EOS but has been heated by
* feedback since it last went through the cooling routines, then we need
* to re-set its abundance array to be in equilibrium. */
const struct globalVariables *ChimesGlobalVars =
&cooling->ChimesGlobalVars;
/* Create ChimesGasVars struct */
struct gasVariables ChimesGasVars;
allocate_gas_abundances_memory(&ChimesGasVars, ChimesGlobalVars);
/* Copy abundances over from xp to ChimesGasVars */
for (int i = 0; i < ChimesGlobalVars->totalNumberOfSpecies; i++)
ChimesGasVars.abundances[i] = xp->cooling_data.chimes_abundances[i];
/* Update element abundances */
chimes_update_element_abundances(phys_const, us, cosmo, cooling,
hydro_props, p, xp, &ChimesGasVars, u,
/*mode=*/1);
const double u_cgs = u * cooling->internal_energy_to_cgs;
/* Update ChimesGasVars */
chimes_update_gas_vars(u_cgs, phys_const, us, cosmo, hydro_props,
floor_props, cooling, dp, p, xp, &ChimesGasVars,
1.0);
/* Set abundances to equilibrium */
cooling_set_FB_particle_chimes_abundances(&ChimesGasVars, cooling);
/* Copy abundances from ChimesGasVars to xp. */
for (int i = 0; i < ChimesGlobalVars->totalNumberOfSpecies; i++)
xp->cooling_data.chimes_abundances[i] = ChimesGasVars.abundances[i];
/* Re-calculate the net cooling rate. */
const float rho_phys_cgs =
hydro_get_physical_density(p, cosmo) * cooling->density_to_cgs;
/* cooling rate in erg/cm^3/s. */
float net_cooling_rate_cgs = (float)chimes_cooling_rate_external(
&ChimesGasVars, ChimesGlobalVars, /*mode=*/0);
/* Convert to erg/g/s */
net_cooling_rate_cgs /= rho_phys_cgs;
/* Convert to code units. */
xp->cooling_data.net_cooling_rate =
net_cooling_rate_cgs * cooling->internal_energy_over_time_from_cgs;
/* Free CHIMES memory. */
free_gas_abundances_memory(&ChimesGasVars, ChimesGlobalVars);
}
/* Re-set heated_by_FB flag. */
xp->cooling_data.heated_by_FB = 0;
}
}
float cooling_get_particle_subgrid_HI_fraction(
const struct unit_system *us, const struct phys_const *phys_const,
const struct cosmology *cosmo, const struct hydro_props *hydro_props,
const struct entropy_floor_properties *floor_props,
const struct cooling_function_data *cooling, const struct part *p,
const struct xpart *xp) {
const struct globalVariables *ChimesGlobalVars = &cooling->ChimesGlobalVars;
return xp->cooling_data
.chimes_abundances[ChimesGlobalVars->speciesIndices[sp_HI]];
}
float cooling_get_particle_subgrid_HII_fraction(
const struct unit_system *us, const struct phys_const *phys_const,
const struct cosmology *cosmo, const struct hydro_props *hydro_props,
const struct entropy_floor_properties *floor_props,
const struct cooling_function_data *cooling, const struct part *p,
const struct xpart *xp) {
const struct globalVariables *ChimesGlobalVars = &cooling->ChimesGlobalVars;
return xp->cooling_data
.chimes_abundances[ChimesGlobalVars->speciesIndices[sp_HII]];
}
float cooling_get_particle_subgrid_H2_fraction(
const struct unit_system *us, const struct phys_const *phys_const,
const struct cosmology *cosmo, const struct hydro_props *hydro_props,
const struct entropy_floor_properties *floor_props,
const struct cooling_function_data *cooling, const struct part *p,
const struct xpart *xp) {
const struct globalVariables *ChimesGlobalVars = &cooling->ChimesGlobalVars;
return xp->cooling_data
.chimes_abundances[ChimesGlobalVars->speciesIndices[sp_H2]];
}
/**
* @brief Compute the temperature of the gas.
*
* For particles on the entropy floor, we use pressure equilibrium to
* infer the properties of the particle.
*
* @param us The internal system of units.
* @param phys_const The physical constants.
* @param hydro_props The properties of the hydro scheme.
* @param cosmo The cosmological model.
* @param floor_props The properties of the entropy floor.
* @param cooling The properties of the cooling scheme.
* @param p The #part.
* @param xp The #xpart.
*/
float cooling_get_subgrid_temperature(
const struct unit_system *us, const struct phys_const *phys_const,
const struct cosmology *cosmo, const struct hydro_props *hydro_props,
const struct entropy_floor_properties *floor_props,
const struct cooling_function_data *cooling, const struct part *p,
const struct xpart *xp) {
/* Get the particle's temperature */
const float T = cooling_get_temperature(phys_const, hydro_props, us, cosmo,
cooling, p, xp);
if (cooling->use_colibre_subgrid_EOS == 1) {
/* Particle's internal energy */
const double u = hydro_get_physical_internal_energy(p, xp, cosmo);
const float log10_T = log10f(T);
/* Physical density of this particle */
const float rho_phys = hydro_get_physical_density(p, cosmo);
/* Recover the Hydrogen mass fraction */
const float XH = get_hydrogen_mass_fraction(p, hydro_props);
/* Get the particle pressure */
const float P_phys = hydro_get_physical_pressure(p, cosmo);
/* Are we in an HII region? */
const int HII_region = xp->tracers_data.HIIregion_timer_gas > 0.;
/* Get the total metallicity in units of solar */
float abundance_ratio[colibre_cooling_N_elementtypes];
const float logZZsol =
abundance_ratio_to_solar(p, &(cooling->colibre_table), abundance_ratio);
/* Limit imposed by the entropy floor */
const double A_EOS = entropy_floor(p, cosmo, floor_props);
const double u_EOS = gas_internal_energy_from_entropy(rho_phys, A_EOS);
const double u_EOS_max = u_EOS * exp10(cooling->dlogT_EOS);
const float log10_u_EOS_max_cgs = log10f(
u_EOS_max * cooling->colibre_table.internal_energy_to_cgs + FLT_MIN);
return compute_subgrid_property(
cooling, phys_const, floor_props, cosmo, rho_phys, logZZsol, XH, P_phys,
log10_T, log10_u_EOS_max_cgs, HII_region, abundance_ratio,
log10(u * cooling->colibre_table.internal_energy_to_cgs),
colibre_compute_subgrid_temperature);
} else {
return T;
}
}
/**
* @brief Compute the physical density of the gas.
*
* For particles on the entropy floor, we use pressure equilibrium to
* infer the properties of the particle.
*
* Note that we return the density in physical coordinates.
*
* @param us The internal system of units.
* @param phys_const The physical constants.
* @param hydro_props The properties of the hydro scheme.
* @param cosmo The cosmological model.
* @param floor_props The properties of the entropy floor.
* @param cooling The properties of the cooling scheme.
* @param p The #part.
* @param xp The #xpart.
*/
float cooling_get_subgrid_density(
const struct unit_system *us, const struct phys_const *phys_const,
const struct cosmology *cosmo, const struct hydro_props *hydro_props,
const struct entropy_floor_properties *floor_props,
const struct cooling_function_data *cooling, const struct part *p,
const struct xpart *xp) {
if (cooling->use_colibre_subgrid_EOS == 1) {
/* Particle's internal energy */
const double u = hydro_get_physical_internal_energy(p, xp, cosmo);
/* Particle's temperature */
const float T = cooling_get_temperature(phys_const, hydro_props, us, cosmo,
cooling, p, xp);
const float log10_T = log10f(T);
/* Physical density of this particle */
const float rho_phys = hydro_get_physical_density(p, cosmo);
/* Recover the Hydrogen mass fraction */
const float XH = get_hydrogen_mass_fraction(p, hydro_props);
/* Get the particle pressure */
const float P_phys = hydro_get_physical_pressure(p, cosmo);
/* Are we in an HII region? */
const int HII_region = xp->tracers_data.HIIregion_timer_gas > 0.;
/* Get the total metallicity in units of solar */
float abundance_ratio[colibre_cooling_N_elementtypes];
const float logZZsol =
abundance_ratio_to_solar(p, &(cooling->colibre_table), abundance_ratio);
/* Limit imposed by the entropy floor */
const double A_EOS = entropy_floor(p, cosmo, floor_props);
const double u_EOS = gas_internal_energy_from_entropy(rho_phys, A_EOS);
const double u_EOS_max = u_EOS * exp10(cooling->dlogT_EOS);
const float log10_u_EOS_max_cgs = log10f(
u_EOS_max * cooling->colibre_table.internal_energy_to_cgs + FLT_MIN);
return compute_subgrid_property(
cooling, phys_const, floor_props, cosmo, rho_phys, logZZsol, XH, P_phys,
log10_T, log10_u_EOS_max_cgs, HII_region, abundance_ratio,
log10(u * cooling->colibre_table.internal_energy_to_cgs),
colibre_compute_subgrid_density);
} else {
return hydro_get_physical_density(p, cosmo);
}
}
/**
* @brief Set CHIMES abundances for HII regions.
*
* Sets the CHIMES abundance array according to the
* HIIregion_fion parameter. This parameter gives the
* fraction of each element that is singly ionised,
* with the remainder set to neutral.
*
* @param ChimesGasVars CHIMES gasVariables structure.
* @param cooling The #cooling_function_data used in the run.
*/
void cooling_set_HIIregion_chimes_abundances(
struct gasVariables *ChimesGasVars,
const struct cooling_function_data *cooling) {
const struct globalVariables *ChimesGlobalVars = &cooling->ChimesGlobalVars;
/* First zero all abundances */
for (int i = 0; i < ChimesGlobalVars->totalNumberOfSpecies; i++)
ChimesGasVars->abundances[i] = 0.0;
/* Set H and He ionisation fractions. */
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_HII]] =
(ChimesFloat)cooling->HIIregion_fion;
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_HI]] =
(ChimesFloat)(1.0 -
ChimesGasVars
->abundances[ChimesGlobalVars->speciesIndices[sp_HII]]);
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_elec]] +=
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_HII]];
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_HeII]] =
(ChimesFloat)ChimesGasVars->element_abundances[0] *
cooling->HIIregion_fion;
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_HeI]] =
(ChimesFloat)ChimesGasVars->element_abundances[0] *
(1.0 - cooling->HIIregion_fion);
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_elec]] +=
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_HeII]];
/* Set metal ionisation fractions. */
if ((ChimesGlobalVars->element_included[0] == 1) &&
(ChimesGasVars->element_abundances[1] > 0.0)) {
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_CII]] =
(ChimesFloat)ChimesGasVars->element_abundances[1] *
cooling->HIIregion_fion;
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_CI]] =
(ChimesFloat)ChimesGasVars->element_abundances[1] *
(1.0 - cooling->HIIregion_fion);
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_elec]] +=
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_CII]];
}
if ((ChimesGlobalVars->element_included[1] == 1) &&
(ChimesGasVars->element_abundances[2] > 0.0)) {
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_NII]] =
(ChimesFloat)ChimesGasVars->element_abundances[2] *
cooling->HIIregion_fion;
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_NI]] =
(ChimesFloat)ChimesGasVars->element_abundances[2] *
(1.0 - cooling->HIIregion_fion);
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_elec]] +=
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_NII]];
}
if ((ChimesGlobalVars->element_included[2] == 1) &&
(ChimesGasVars->element_abundances[3] > 0.0)) {
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_OII]] =
(ChimesFloat)ChimesGasVars->element_abundances[3] *
cooling->HIIregion_fion;
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_OI]] =
(ChimesFloat)ChimesGasVars->element_abundances[3] *
(1.0 - cooling->HIIregion_fion);
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_elec]] +=
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_OII]];
}
if ((ChimesGlobalVars->element_included[3] == 1) &&
(ChimesGasVars->element_abundances[4] > 0.0)) {
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_NeII]] =
(ChimesFloat)ChimesGasVars->element_abundances[4] *
cooling->HIIregion_fion;
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_NeI]] =
(ChimesFloat)ChimesGasVars->element_abundances[4] *
(1.0 - cooling->HIIregion_fion);
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_elec]] +=
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_NeII]];
}
if ((ChimesGlobalVars->element_included[4] == 1) &&
(ChimesGasVars->element_abundances[5] > 0.0)) {
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_MgII]] =
(ChimesFloat)ChimesGasVars->element_abundances[5] *
cooling->HIIregion_fion;
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_MgI]] =
(ChimesFloat)ChimesGasVars->element_abundances[5] *
(1.0 - cooling->HIIregion_fion);
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_elec]] +=
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_MgII]];
}
if ((ChimesGlobalVars->element_included[5] == 1) &&
(ChimesGasVars->element_abundances[6] > 0.0)) {
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_SiII]] =
(ChimesFloat)ChimesGasVars->element_abundances[6] *
cooling->HIIregion_fion;
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_SiI]] =
(ChimesFloat)ChimesGasVars->element_abundances[6] *
(1.0 - cooling->HIIregion_fion);
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_elec]] +=
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_SiII]];
}
if ((ChimesGlobalVars->element_included[6] == 1) &&
(ChimesGasVars->element_abundances[7] > 0.0)) {
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_SII]] =
(ChimesFloat)ChimesGasVars->element_abundances[7] *
cooling->HIIregion_fion;
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_SI]] =
(ChimesFloat)ChimesGasVars->element_abundances[7] *
(1.0 - cooling->HIIregion_fion);
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_elec]] +=
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_SII]];
}
if ((ChimesGlobalVars->element_included[7] == 1) &&
(ChimesGasVars->element_abundances[8] > 0.0)) {
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_CaII]] =
(ChimesFloat)ChimesGasVars->element_abundances[8] *
cooling->HIIregion_fion;
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_CaI]] =
(ChimesFloat)ChimesGasVars->element_abundances[8] *
(1.0 - cooling->HIIregion_fion);
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_elec]] +=
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_CaII]];
}
if ((ChimesGlobalVars->element_included[8] == 1) &&
(ChimesGasVars->element_abundances[9] > 0.0)) {
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_FeII]] =
(ChimesFloat)ChimesGasVars->element_abundances[9] *
cooling->HIIregion_fion;
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_FeI]] =
(ChimesFloat)ChimesGasVars->element_abundances[9] *
(1.0 - cooling->HIIregion_fion);
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_elec]] +=
ChimesGasVars->abundances[ChimesGlobalVars->speciesIndices[sp_FeII]];
}
}
/**
* @brief Set CHIMES abundances for particles heated by feedback.
*
* Sets gas particle that has just been heated by feedback
* to be in chemical equilibrium at its new temperature.
*
* @param ChimesGasVars CHIMES gasVariables structure.
* @param cooling The #cooling_function_data used in the run.
*/
void cooling_set_FB_particle_chimes_abundances(
struct gasVariables *ChimesGasVars,
const struct cooling_function_data *cooling) {
const struct globalVariables *ChimesGlobalVars = &cooling->ChimesGlobalVars;
/* Save equilibrium and ThermEvol flags. */
int ForceEqOn_save = ChimesGasVars->ForceEqOn;
int ThermEvolOn_save = ChimesGasVars->ThermEvolOn;
/* Use equilibrium abundances. */
ChimesGasVars->ForceEqOn = 1;
/* Disable temperature evolution. */
ChimesGasVars->ThermEvolOn = 0;
/* Set abundance array to equilibrium. */
chimes_network(ChimesGasVars, ChimesGlobalVars);
/* Revert flags to saved values. */
ChimesGasVars->ForceEqOn = ForceEqOn_save;
ChimesGasVars->ThermEvolOn = ThermEvolOn_save;
}