Commit b934599d authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Back-port the changes to the COLIBRE cooling to better handle subgrid...

Back-port the changes to the COLIBRE cooling to better handle subgrid properties when the entropy floor is high
parent d02e9667
......@@ -513,14 +513,18 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
const double gas_P_phys =
gas_pressure_from_internal_energy(gas_rho_phys, gas_u_phys);
/* Assume primordial abundance and solar metallicity
/* Assume primordial abundance and solar metallicity pattern
* (Yes, that is inconsitent but makes no difference) */
const double logZZsol = 0.;
const double XH = 0.75;
float abundance_ratio[colibre_cooling_N_elementtypes];
for (int i = 0; i < colibre_cooling_N_elementtypes; ++i)
abundance_ratio[i] = 1.f;
/* Get the gas temperature */
const float gas_T = cooling_get_temperature_from_gas(
constants, cosmo, cooling, gas_rho_phys, logZZsol, XH, gas_u_phys);
constants, cosmo, cooling, gas_rho_phys, logZZsol, XH, gas_u_phys,
/*HII_region=*/0);
const float log10_gas_T = log10f(gas_T);
/* Get the temperature on the EOS at this physical density */
......@@ -533,9 +537,10 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
/* Compute the subgrid density assuming pressure
* equilibirum if on the entropy floor */
const double rho_sub = compute_subgrid_density(
const double rho_sub = compute_subgrid_property(
cooling, constants, floor_props, cosmo, gas_rho_phys, logZZsol, XH,
gas_P_phys, log10_gas_T, log10_T_EOS_max);
gas_P_phys, log10_gas_T, log10_T_EOS_max, /*HII_region=*/0,
abundance_ratio, 0.f, colibre_compute_subgrid_density);
/* Record what we used */
bp->rho_subgrid_gas = rho_sub;
......
......@@ -95,6 +95,10 @@ void cooling_update(const struct cosmology *cosmo,
* @brief Compute the internal energy of a #part based on the cooling function
* but for a given temperature.
*
* This is used e.g. for particles in HII regions that are set to a constant
* temperature, but their internal energies should reflect the particle
* composition .
*
* @param phys_const #phys_const data structure.
* @param hydro_props The properties of the hydro scheme.
* @param us The internal system of units.
......@@ -122,6 +126,8 @@ float cooling_get_internalenergy_for_temperature(
chemistry_get_metal_mass_fraction_for_cooling(p);
const float XH = metal_fraction[chemistry_element_H];
/* Normal case --> Interpolate the table */
/* Convert Hydrogen mass fraction into Hydrogen number density */
const float rho = hydro_get_physical_density(p, cosmo);
const double n_H = rho * XH / phys_const->const_proton_mass;
......@@ -169,20 +175,28 @@ float cooling_get_internalenergy_for_temperature(
* metallicity.
* @param XH The Hydrogen abundance of the gas.
* @param u_phys Internal energy of the gas in internal physical units.
* @param HII_region Is this patch of gas in an HII region?
*/
float cooling_get_temperature_from_gas(
const struct phys_const *phys_const, const struct cosmology *cosmo,
const struct cooling_function_data *cooling, const float rho_phys,
const float logZZsol, const float XH, const float u_phys) {
const float logZZsol, const float XH, const float u_phys,
const int HII_region) {
if (HII_region)
error("HII regions are not implemented in the EAGLE-XL flavour");
/* Convert to CGS */
const double u_cgs = u_phys * cooling->internal_energy_to_cgs;
/* Get density in Hydrogen number density */
const double n_H = rho_phys * XH / phys_const->const_proton_mass;
const double n_H_cgs = n_H * cooling->number_density_to_cgs;
/* Normal case --> Interpolate the table */
/* compute hydrogen number density, metallicity and redshift indices and
* offsets */
float d_red, d_met, d_n_H;
int red_index, met_index, n_H_index;
......@@ -205,7 +219,8 @@ float cooling_get_temperature_from_gas(
/**
* @brief Compute the temperature of a #part based on the cooling function.
*
* The temperature returned is consistent with the cooling rates.
* The temperature returned is consistent with the cooling rates or
* is the temperature of an HII region if the particle is flagged as such.
*
* @param phys_const #phys_const data structure.
* @param hydro_props The properties of the hydro scheme.
......@@ -229,17 +244,15 @@ float cooling_get_temperature(const struct phys_const *phys_const,
"--temperature runtime flag?");
#endif
/* Get physical internal energy */
/* 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);
/* Get the Hydrogen mass fraction */
float const *metal_fraction =
chemistry_get_metal_mass_fraction_for_cooling(p);
const float XH = metal_fraction[chemistry_element_H];
/* Convert Hydrogen mass fraction into Hydrogen number density */
const float rho_phys = hydro_get_physical_density(p, cosmo);
/* Get this particle's metallicity ratio to solar.
*
* Note that we do not need the individual element's ratios that
......@@ -247,8 +260,11 @@ float cooling_get_temperature(const struct phys_const *phys_const,
float dummy[colibre_cooling_N_elementtypes];
const float logZZsol = abundance_ratio_to_solar(p, cooling, dummy);
/* Are we in an HII region? */
const int HII_region = 0; /* No HII regions in the EAGLE-XL flavour */
return cooling_get_temperature_from_gas(phys_const, cosmo, cooling, rho_phys,
logZZsol, XH, u_phys);
logZZsol, XH, u_phys, HII_region);
}
/**
......@@ -721,8 +737,8 @@ float cooling_get_particle_subgrid_HI_fraction(
const float rho_phys = hydro_get_physical_density(p, cosmo);
/* Get the total metallicity in units of solar */
float dummy[colibre_cooling_N_elementtypes];
const float logZZsol = abundance_ratio_to_solar(p, cooling, dummy);
float abundance_ratio[colibre_cooling_N_elementtypes];
const float logZZsol = abundance_ratio_to_solar(p, cooling, abundance_ratio);
/* Get the Hydrogen abundance */
const float *const metal_fraction =
......@@ -735,14 +751,19 @@ float cooling_get_particle_subgrid_HI_fraction(
/* Get physical internal energy */
const float u_phys = hydro_get_physical_internal_energy(p, xp, cosmo);
/* Are we in an HII region? */
const int HII_region = 0; /* No HII regions in the EAGLE-XL flavour */
/* Get the particle's temperature */
const float T = cooling_get_temperature_from_gas(
phys_const, cosmo, cooling, rho_phys, logZZsol, XH, u_phys);
phys_const, cosmo, cooling, rho_phys, logZZsol, XH, u_phys, HII_region);
const float log10_T = log10f(T);
return compute_subgrid_HI_fraction(cooling, phys_const, floor_props, cosmo,
rho_phys, logZZsol, XH, P_phys, log10_T,
log10_T_EOS_max);
return compute_subgrid_property(
cooling, phys_const, floor_props, cosmo, rho_phys, logZZsol, XH, P_phys,
log10_T, log10_T_EOS_max, HII_region, abundance_ratio,
log10(u_phys * cooling->internal_energy_to_cgs),
colibre_compute_subgrid_HI_fraction);
}
/**
......@@ -777,8 +798,8 @@ float cooling_get_particle_subgrid_HII_fraction(
const float rho_phys = hydro_get_physical_density(p, cosmo);
/* Get the total metallicity in units of solar */
float dummy[colibre_cooling_N_elementtypes];
const float logZZsol = abundance_ratio_to_solar(p, cooling, dummy);
float abundance_ratio[colibre_cooling_N_elementtypes];
const float logZZsol = abundance_ratio_to_solar(p, cooling, abundance_ratio);
/* Get the Hydrogen abundance */
const float *const metal_fraction =
......@@ -791,14 +812,19 @@ float cooling_get_particle_subgrid_HII_fraction(
/* Get physical internal energy */
const float u_phys = hydro_get_physical_internal_energy(p, xp, cosmo);
/* Are we in an HII region? */
const int HII_region = 0; /* No HII regions in the EAGLE-XL flavour */
/* Get the particle's temperature */
const float T = cooling_get_temperature_from_gas(
phys_const, cosmo, cooling, rho_phys, logZZsol, XH, u_phys);
phys_const, cosmo, cooling, rho_phys, logZZsol, XH, u_phys, HII_region);
const float log10_T = log10f(T);
return compute_subgrid_HII_fraction(cooling, phys_const, floor_props, cosmo,
rho_phys, logZZsol, XH, P_phys, log10_T,
log10_T_EOS_max);
return compute_subgrid_property(
cooling, phys_const, floor_props, cosmo, rho_phys, logZZsol, XH, P_phys,
log10_T, log10_T_EOS_max, HII_region, abundance_ratio,
log10(u_phys * cooling->internal_energy_to_cgs),
colibre_compute_subgrid_HII_fraction);
}
/**
......@@ -833,8 +859,8 @@ float cooling_get_particle_subgrid_H2_fraction(
const float rho_phys = hydro_get_physical_density(p, cosmo);
/* Get the total metallicity in units of solar */
float dummy[colibre_cooling_N_elementtypes];
const float logZZsol = abundance_ratio_to_solar(p, cooling, dummy);
float abundance_ratio[colibre_cooling_N_elementtypes];
const float logZZsol = abundance_ratio_to_solar(p, cooling, abundance_ratio);
/* Get the Hydrogen abundance */
const float *const metal_fraction =
......@@ -847,14 +873,19 @@ float cooling_get_particle_subgrid_H2_fraction(
/* Get physical internal energy */
const float u_phys = hydro_get_physical_internal_energy(p, xp, cosmo);
/* Are we in an HII region? */
const int HII_region = 0; /* No HII regions in the EAGLE-XL flavour */
/* Get the particle's temperature */
const float T = cooling_get_temperature_from_gas(
phys_const, cosmo, cooling, rho_phys, logZZsol, XH, u_phys);
phys_const, cosmo, cooling, rho_phys, logZZsol, XH, u_phys, HII_region);
const float log10_T = log10f(T);
return compute_subgrid_H2_fraction(cooling, phys_const, floor_props, cosmo,
rho_phys, logZZsol, XH, P_phys, log10_T,
log10_T_EOS_max);
return compute_subgrid_property(
cooling, phys_const, floor_props, cosmo, rho_phys, logZZsol, XH, P_phys,
log10_T, log10_T_EOS_max, HII_region, abundance_ratio,
log10(u_phys * cooling->internal_energy_to_cgs),
colibre_compute_subgrid_H2_fraction);
}
/**
......@@ -888,8 +919,8 @@ float cooling_get_particle_subgrid_temperature(
const float rho_phys = hydro_get_physical_density(p, cosmo);
/* Get the total metallicity in units of solar */
float dummy[colibre_cooling_N_elementtypes];
const float logZZsol = abundance_ratio_to_solar(p, cooling, dummy);
float abundance_ratio[colibre_cooling_N_elementtypes];
const float logZZsol = abundance_ratio_to_solar(p, cooling, abundance_ratio);
/* Get the Hydrogen abundance */
const float *const metal_fraction =
......@@ -902,14 +933,19 @@ float cooling_get_particle_subgrid_temperature(
/* Get physical internal energy */
const float u_phys = hydro_get_physical_internal_energy(p, xp, cosmo);
/* Are we in an HII region? */
const int HII_region = 0; /* No HII regions in the EAGLE-XL flavour */
/* Get the particle's temperature */
const float T = cooling_get_temperature_from_gas(
phys_const, cosmo, cooling, rho_phys, logZZsol, XH, u_phys);
phys_const, cosmo, cooling, rho_phys, logZZsol, XH, u_phys, HII_region);
const float log10_T = log10f(T);
return compute_subgrid_temperature(cooling, phys_const, floor_props, cosmo,
rho_phys, logZZsol, XH, P_phys, log10_T,
log10_T_EOS_max);
return compute_subgrid_property(
cooling, phys_const, floor_props, cosmo, rho_phys, logZZsol, XH, P_phys,
log10_T, log10_T_EOS_max, HII_region, abundance_ratio,
log10(u_phys * cooling->internal_energy_to_cgs),
colibre_compute_subgrid_temperature);
}
/**
......@@ -945,8 +981,8 @@ float cooling_get_particle_subgrid_density(
const float rho_phys = hydro_get_physical_density(p, cosmo);
/* Get the total metallicity in units of solar */
float dummy[colibre_cooling_N_elementtypes];
const float logZZsol = abundance_ratio_to_solar(p, cooling, dummy);
float abundance_ratio[colibre_cooling_N_elementtypes];
const float logZZsol = abundance_ratio_to_solar(p, cooling, abundance_ratio);
/* Get the Hydrogen abundance */
const float *const metal_fraction =
......@@ -956,17 +992,21 @@ float cooling_get_particle_subgrid_density(
/* Get the particle pressure */
const float P_phys = hydro_get_physical_pressure(p, cosmo);
/* Get physical internal energy */
const float u_phys = hydro_get_physical_internal_energy(p, xp, cosmo);
/* Get the particle's temperature */
const float T = cooling_get_temperature_from_gas(
phys_const, cosmo, cooling, rho_phys, logZZsol, XH, u_phys);
const float T = cooling_get_temperature(phys_const, hydro_props, us, cosmo,
cooling, p, xp);
const float log10_T = log10f(T);
return compute_subgrid_density(cooling, phys_const, floor_props, cosmo,
rho_phys, logZZsol, XH, P_phys, log10_T,
log10_T_EOS_max);
/* Are we in an HII region? */
const int HII_region = 0; /* No HII regions in the EAGLE-XL flavour */
const float u_start = hydro_get_physical_internal_energy(p, xp, cosmo);
const double u_0_cgs = u_start * cooling->internal_energy_to_cgs;
return compute_subgrid_property(
cooling, phys_const, floor_props, cosmo, rho_phys, logZZsol, XH, P_phys,
log10_T, log10_T_EOS_max, HII_region, abundance_ratio, log10(u_0_cgs),
colibre_compute_subgrid_density);
}
/**
......@@ -997,8 +1037,8 @@ void cooling_set_particle_subgrid_properties(
const float rho_phys = hydro_get_physical_density(p, cosmo);
/* Get the total metallicity in units of solar */
float dummy[colibre_cooling_N_elementtypes];
const float logZZsol = abundance_ratio_to_solar(p, cooling, dummy);
float abundance_ratio[colibre_cooling_N_elementtypes];
const float logZZsol = abundance_ratio_to_solar(p, cooling, abundance_ratio);
/* Get the Hydrogen abundance */
const float *const metal_fraction =
......@@ -1010,18 +1050,25 @@ void cooling_set_particle_subgrid_properties(
/* Get physical internal energy */
const float u_phys = hydro_get_physical_internal_energy(p, xp, cosmo);
const double u_cgs = u_phys * cooling->internal_energy_to_cgs;
/* Are we in an HII region? */
const int HII_region = 0; /* No HII regions in the EAGLE-XL flavour */
/* Get the particle's temperature */
const float T = cooling_get_temperature_from_gas(
phys_const, cosmo, cooling, rho_phys, logZZsol, XH, u_phys);
phys_const, cosmo, cooling, rho_phys, logZZsol, XH, u_phys, HII_region);
const float log10_T = log10f(T);
p->cooling_data.subgrid_temp = compute_subgrid_temperature(
p->cooling_data.subgrid_temp = compute_subgrid_property(
cooling, phys_const, floor_props, cosmo, rho_phys, logZZsol, XH, P_phys,
log10_T, log10_T_EOS_max, HII_region, abundance_ratio, log10(u_cgs),
colibre_compute_subgrid_temperature);
p->cooling_data.subgrid_dens = compute_subgrid_property(
cooling, phys_const, floor_props, cosmo, rho_phys, logZZsol, XH, P_phys,
log10_T, log10_T_EOS_max);
p->cooling_data.subgrid_dens =
compute_subgrid_density(cooling, phys_const, floor_props, cosmo, rho_phys,
logZZsol, XH, P_phys, log10_T, log10_T_EOS_max);
log10_T, log10_T_EOS_max, HII_region, abundance_ratio, log10(u_cgs),
colibre_compute_subgrid_density);
}
/**
......@@ -1153,6 +1200,8 @@ void cooling_init_backend(struct swift_params *parameter_file,
cooling->number_density_to_cgs =
units_cgs_conversion_factor(us, UNIT_CONV_NUMBER_DENSITY);
cooling->number_density_from_cgs = 1. / cooling->number_density_to_cgs;
cooling->density_to_cgs = units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
cooling->density_from_cgs = 1. / cooling->density_to_cgs;
/* Store some constants in CGS units */
const float units_kB[5] = {1, 2, -2, 0, -1};
......@@ -1164,6 +1213,7 @@ void cooling_init_backend(struct swift_params *parameter_file,
cooling->log10_kB_cgs = log10(kB_cgs);
cooling->inv_proton_mass_cgs = 1. / proton_mass_cgs;
cooling->proton_mass_cgs = proton_mass_cgs;
cooling->T_CMB_0 = phys_const->const_T_CMB_0 *
units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
......
......@@ -26,6 +26,7 @@
/* Local includes. */
#include "cooling_struct.h"
#include "cooling_tables.h"
struct part;
struct xpart;
......@@ -64,7 +65,8 @@ void cooling_first_init_part(const struct phys_const *phys_const,
float cooling_get_temperature_from_gas(
const struct phys_const *phys_const, const struct cosmology *cosmo,
const struct cooling_function_data *cooling, const float rho_phys,
const float XH, const float logZZsol, const float u_phys);
const float XH, const float logZZsol, const float u_phys,
const int HII_region);
float cooling_get_temperature(const struct phys_const *phys_const,
const struct hydro_props *hydro_props,
......@@ -115,13 +117,22 @@ void cooling_set_particle_subgrid_properties(
const struct cooling_function_data *cooling, struct part *p,
struct xpart *xp);
double compute_subgrid_density(
double get_thermal_equilibrium_pressure(
const struct cooling_function_data *cooling,
const float abundance_ratio[colibre_cooling_N_elementtypes],
const double log_u_cgs, const float log10nH_local, const double rho_cgs,
const double redshift, const int ired, const int imet, const float dred,
const float dmet);
double compute_subgrid_property(
const struct cooling_function_data *cooling,
const struct phys_const *phys_const,
const struct entropy_floor_properties *floor_props,
const struct cosmology *cosmo, const float rho_phys, const float logZZsol,
const float XH, const float P_phys, const float log10_T,
const float log10_T_EOS_max);
const float log10_T_EOS_max, const int HII_region,
const float abundance_ratio[colibre_cooling_N_elementtypes],
const double log_u_cgs, const enum colibre_subgrid_properties isub);
float cooling_get_radiated_energy(const struct xpart *xp);
......@@ -141,4 +152,39 @@ void cooling_print_backend(const struct cooling_function_data *cooling);
void cooling_clean(struct cooling_function_data *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 does not do anything.
*
* @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.
*/
__attribute__((always_inline)) INLINE static void cooling_convert_quantities(
struct part *restrict p, struct xpart *restrict 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) {}
/**
* @brief Updates cooling properties of particle hit by feedback.
*
* For this cooling module, this routine does not do anything.
*
* @param xp Pointer to the extended particle data.
*/
__attribute__((always_inline)) INLINE static void
cooling_update_feedback_particle(struct xpart *restrict xp) {}
#endif /* SWIFT_COOLING_COLIBRE_H */
......@@ -165,9 +165,18 @@ struct cooling_function_data {
/*! Number density conversion from CGS to internal units (for quick access) */
double number_density_from_cgs;
/*! Density conversion from internal units to CGS (for quick access) */
double density_to_cgs;
/*! Density conversion from CGS to internal units (for quick access) */
double density_from_cgs;
/*! Inverse of proton mass in cgs (for quick access) */
double inv_proton_mass_cgs;
/*! Proton mass in cgs (for quick access) */
double proton_mass_cgs;
/*! Logarithm base 10 of the Boltzmann constant in CGS (for quick access) */
double log10_kB_cgs;
......
This diff is collapsed.
......@@ -55,6 +55,17 @@
/*! Number of different elements in the tables */
#define colibre_cooling_N_elementtypes 12
/**
* @brief Subgrid properties to calculate
*/
enum colibre_subgrid_properties {
colibre_compute_subgrid_density,
colibre_compute_subgrid_temperature,
colibre_compute_subgrid_HI_fraction,
colibre_compute_subgrid_HII_fraction,
colibre_compute_subgrid_H2_fraction
};
/**
* @brief Elements present in the tables
*/
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment