Commit 0ffe8095 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'yet_another_agn_model' into 'master'

Yet another agn model

See merge request !1114
parents 22b1ff19 a30a6622
......@@ -509,6 +509,7 @@ EAGLEAGN:
use_subgrid_mass_from_ics: 1 # (Optional) Use subgrid masses specified in ICs [1, default], or initialise them to particle masses [0]?
with_subgrid_mass_check: 1 # (Optional) Verify that initial black hole subgrid masses are positive [1, default]. Only used if use_subgrid_mass_from_ics is 1.
multi_phase_bondi: 0 # Compute Bondi rates per neighbour particle?
subgrid_bondi: 0 # Compute Bondi rates using the subgrid extrapolation of the gas properties around the BH?
with_angmom_limiter: 1 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term?
viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term
radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
......
......@@ -185,6 +185,8 @@ __attribute__((always_inline)) INLINE static void black_holes_swallow_bpart(
* @param props The properties of the black hole scheme.
* @param constants The physical constants (in internal units).
* @param cosmo The cosmological model.
* @param cooling Properties of the cooling model.
* @param floor_props Properties of the entropy fllor.
* @param time Time since the start of the simulation (non-cosmo mode).
* @param with_cosmology Are we running with cosmology?
* @param dt The time-step size (in physical internal units).
......@@ -192,7 +194,9 @@ __attribute__((always_inline)) INLINE static void black_holes_swallow_bpart(
__attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
struct bpart* restrict bp, const struct black_holes_props* props,
const struct phys_const* constants, const struct cosmology* cosmo,
const double time, const int with_cosmology, const double dt) {}
const struct cooling_function_data* cooling,
const struct entropy_floor_properties* floor_props, const double time,
const int with_cosmology, const double dt) {}
/**
* @brief Finish the calculation of the new BH position.
......@@ -203,6 +207,7 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
* @param props The properties of the black hole scheme.
* @param constants The physical constants (in internal units).
* @param cosmo The cosmological model.
* @param dt The time-step size (in physical internal units).
*/
__attribute__((always_inline)) INLINE static void black_holes_end_reposition(
struct bpart* restrict bp, const struct black_holes_props* props,
......
......@@ -22,6 +22,7 @@
/* Local includes */
#include "black_holes_properties.h"
#include "black_holes_struct.h"
#include "cooling.h"
#include "cosmology.h"
#include "dimension.h"
#include "gravity.h"
......@@ -108,6 +109,9 @@ __attribute__((always_inline)) INLINE static void black_holes_init_bpart(
bp->density.wcount_dh = 0.f;
bp->rho_gas = 0.f;
bp->sound_speed_gas = 0.f;
bp->internal_energy_gas = 0.f;
bp->rho_subgrid_gas = -1.f;
bp->sound_speed_subgrid_gas = -1.f;
bp->velocity_gas[0] = 0.f;
bp->velocity_gas[1] = 0.f;
bp->velocity_gas[2] = 0.f;
......@@ -223,6 +227,7 @@ __attribute__((always_inline)) INLINE static void black_holes_end_density(
/* For the following, we also have to undo the mass smoothing
* (N.B.: bp->velocity_gas is in BH frame, in internal units). */
bp->sound_speed_gas *= h_inv_dim * rho_inv;
bp->internal_energy_gas *= h_inv_dim * rho_inv;
bp->velocity_gas[0] *= h_inv_dim * rho_inv;
bp->velocity_gas[1] *= h_inv_dim * rho_inv;
bp->velocity_gas[2] *= h_inv_dim * rho_inv;
......@@ -417,6 +422,8 @@ __attribute__((always_inline)) INLINE static void black_holes_swallow_bpart(
* @param props The properties of the black hole scheme.
* @param constants The physical constants (in internal units).
* @param cosmo The cosmological model.
* @param cooling Properties of the cooling model.
* @param floor_props Properties of the entropy fllor.
* @param time Time since the start of the simulation (non-cosmo mode).
* @param with_cosmology Are we running with cosmology?
* @param dt The time-step size (in physical internal units).
......@@ -424,7 +431,9 @@ __attribute__((always_inline)) INLINE static void black_holes_swallow_bpart(
__attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
struct bpart* restrict bp, const struct black_holes_props* props,
const struct phys_const* constants, const struct cosmology* cosmo,
const double time, const int with_cosmology, const double dt) {
const struct cooling_function_data* cooling,
const struct entropy_floor_properties* floor_props, const double time,
const int with_cosmology, const double dt) {
/* Record that the black hole has another active time step */
bp->number_of_time_steps++;
......@@ -453,8 +462,8 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
/* Convert the quantities we gathered to physical frame (all internal units).
* Note: for the velocities this means peculiar velocities */
const double gas_c_phys = bp->sound_speed_gas * cosmo->a_factor_sound_speed;
const double gas_c_phys2 = gas_c_phys * gas_c_phys;
double gas_c_phys = bp->sound_speed_gas * cosmo->a_factor_sound_speed;
double gas_c_phys2 = gas_c_phys * gas_c_phys;
const double gas_v_circular[3] = {
bp->circular_velocity_gas[0] * cosmo->a_inv,
bp->circular_velocity_gas[1] * cosmo->a_inv,
......@@ -480,11 +489,10 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
(4. * M_PI * G * G * BH_mass * BH_mass * hi_inv_dim);
} else {
/* Standard approach: compute accretion rate for all gas simultaneously.
*
* Convert the quantities we gathered to physical frame (all internal
* units). Note: velocities are already in black hole frame. */
const double gas_rho_phys = bp->rho_gas * cosmo->a3_inv;
/* Standard approach: compute accretion rate for all gas simultaneously */
/* Convert velocities to physical frame
* Note: velocities are already in black hole frame. */
const double gas_v_phys[3] = {bp->velocity_gas[0] * cosmo->a_inv,
bp->velocity_gas[1] * cosmo->a_inv,
bp->velocity_gas[2] * cosmo->a_inv};
......@@ -493,18 +501,89 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
gas_v_phys[1] * gas_v_phys[1] +
gas_v_phys[2] * gas_v_phys[2];
const double denominator2 = gas_v_norm2 + gas_c_phys2;
if (props->subgrid_bondi) {
/* Use subgrid rho and c for Bondi model */
/* Construct basic properties of the gas around the BH in
physical coordinates */
const double gas_rho_phys = bp->rho_gas * cosmo->a3_inv;
const double gas_u_phys =
bp->internal_energy_gas * cosmo->a_factor_internal_energy;
const double gas_P_phys =
gas_pressure_from_internal_energy(gas_rho_phys, gas_u_phys);
/* Assume primordial abundance and solar metallicity
* (Yes, that is inconsitent but makes no difference) */
const double logZZsol = 0.;
const double XH = 0.75;
/* Get the gas temperature */
const float gas_T = cooling_get_temperature_from_gas(
constants, cosmo, cooling, gas_rho_phys, logZZsol, XH, gas_u_phys);
const float log10_gas_T = log10f(gas_T);
/* Get the temperature on the EOS at this physical density */
const float T_EOS = entropy_floor_gas_temperature(
gas_rho_phys, bp->rho_gas, cosmo, floor_props);
/* Add the allowed offset */
const float log10_T_EOS_max =
log10f(max(T_EOS, FLT_MIN)) + cooling->dlogT_EOS;
/* Compute the subgrid density assuming pressure
* equilibirum if on the entropy floor */
const double rho_sub = compute_subgrid_density(
cooling, constants, floor_props, cosmo, gas_rho_phys, logZZsol, XH,
gas_P_phys, log10_gas_T, log10_T_EOS_max);
/* Record what we used */
bp->rho_subgrid_gas = rho_sub;
/* And the subgrid sound-speed */
const float c_sub = gas_soundspeed_from_pressure(rho_sub, gas_P_phys);
/* Also update the sound-speed to use in the angular momentum limiter */
gas_c_phys = c_sub;
gas_c_phys2 = c_sub * c_sub;
/* Record what we used */
bp->sound_speed_subgrid_gas = c_sub;
/* Now, compute the Bondi rate based on the normal velocities and
* the subgrid density and sound-speed */
const double denominator2 = gas_v_norm2 + gas_c_phys2;
#ifdef SWIFT_DEBUG_CHECKS
/* Make sure that the denominator is strictly positive */
if (denominator2 <= 0)
error(
"Invalid denominator for black hole particle %lld in Bondi rate "
"calculation.",
bp->id);
#endif
const double denominator_inv = 1. / sqrt(denominator2);
Bondi_rate = 4. * M_PI * G * G * BH_mass * BH_mass * rho_sub *
denominator_inv * denominator_inv * denominator_inv;
} else {
/* Use dynamical rho and c for Bondi model */
const double gas_rho_phys = bp->rho_gas * cosmo->a3_inv;
const double denominator2 = gas_v_norm2 + gas_c_phys2;
#ifdef SWIFT_DEBUG_CHECKS
/* Make sure that the denominator is strictly positive */
if (denominator2 <= 0)
error(
"Invalid denominator for black hole particle %lld in Bondi rate "
"calculation.",
bp->id);
/* Make sure that the denominator is strictly positive */
if (denominator2 <= 0)
error(
"Invalid denominator for black hole particle %lld in Bondi rate "
"calculation.",
bp->id);
#endif
const double denominator_inv = 1. / sqrt(denominator2);
Bondi_rate = 4. * M_PI * G * G * BH_mass * BH_mass * gas_rho_phys *
denominator_inv * denominator_inv * denominator_inv;
const double denominator_inv = 1. / sqrt(denominator2);
Bondi_rate = 4. * M_PI * G * G * BH_mass * BH_mass * gas_rho_phys *
denominator_inv * denominator_inv * denominator_inv;
}
}
/* Compute the reduction factor from Rosas-Guevara et al. (2015) */
......
......@@ -81,12 +81,18 @@ runner_iact_nonsym_bh_gas_density(
/* Contribution to the total neighbour mass */
bi->ngb_mass += mj;
/* Neighbour sounds speed */
/* Neighbour sound speed */
const float cj = hydro_get_comoving_soundspeed(pj);
/* Contribution to the smoothed sound speed */
bi->sound_speed_gas += mj * cj * wi;
/* Neighbour internal energy */
const float uj = hydro_get_drifted_comoving_internal_energy(pj);
/* Contribution to the smoothed internal energy */
bi->internal_energy_gas += mj * uj * wi;
/* Neighbour's (drifted) velocity in the frame of the black hole
* (we don't include a Hubble term since we are interested in the
* velocity contribution at the location of the black hole) */
......
......@@ -139,7 +139,7 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
int with_cosmology) {
/* Say how much we want to write */
*num_fields = 27;
*num_fields = 29;
/* List what we want to write */
list[0] = io_make_output_field_convert_bpart(
......@@ -178,7 +178,7 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
}
list[7] = io_make_output_field(
"GasDensities", FLOAT, 1, UNIT_CONV_DENSITY, 0.f, bparts, rho_gas,
"GasDensities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f, bparts, rho_gas,
"Co-moving densities of the gas around the particles");
list[8] = io_make_output_field(
......@@ -318,6 +318,16 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
"rates have been suppressed by the Rosas-Guevara et al. (2015) "
"accretion disc model.");
list[27] = io_make_output_field(
"SubgridDensities", FLOAT, 1, UNIT_CONV_DENSITY, 0.f, bparts,
rho_subgrid_gas,
"Physical subgrid densities used in the subgrid-Bondi model.");
list[28] = io_make_output_field(
"SubgridSoundSpeeds", FLOAT, 1, UNIT_CONV_SPEED, 0.f, bparts,
sound_speed_subgrid_gas,
"Physical subgrid sound-speeds used in the subgrid-Bondi model.");
#ifdef DEBUG_INTERACTIONS_BLACK_HOLES
list += *num_fields;
......
......@@ -90,9 +90,20 @@ struct bpart {
/*! Density of the gas surrounding the black hole. */
float rho_gas;
/*! Internal energy of the gas surrounding the black hole. */
float internal_energy_gas;
/*! Smoothed sound speed of the gas surrounding the black hole. */
float sound_speed_gas;
/*! Subgrid physical density of the gas (updated when using the subgrid Bondi
* model) */
float rho_subgrid_gas;
/*! Subgrid physical sound speed of the gas (updated when using the subgrid
* Bondi model) */
float sound_speed_subgrid_gas;
/*! Smoothed velocity of the gas surrounding the black hole,
* in the frame of the black hole (internal units) */
float velocity_gas[3];
......
......@@ -63,6 +63,9 @@ struct black_holes_props {
/*! Calculate Bondi accretion rate for individual neighbours? */
int multi_phase_bondi;
/*! Are we using the subgrid gas properties in the Bondi model? */
int subgrid_bondi;
/*! Are we applying the angular-momentum-based multiplicative term from
* Rosas-Guevara et al. (2015)? */
int with_angmom_limiter;
......@@ -209,6 +212,13 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
bp->multi_phase_bondi =
parser_get_param_int(params, "EAGLEAGN:multi_phase_bondi");
bp->subgrid_bondi = parser_get_param_int(params, "EAGLEAGN:subgrid_bondi");
if (bp->multi_phase_bondi && bp->subgrid_bondi)
error(
"Cannot run with both the multi-phase Bondi and subgrid Bondi models "
"at the same time!");
/* Rosas-Guevara et al. (2015) model */
bp->with_angmom_limiter =
parser_get_param_int(params, "EAGLEAGN:with_angmom_limiter");
......
......@@ -156,6 +156,52 @@ float cooling_get_internalenergy_for_temperature(
return exp10(log_10_U);
}
/**
* @brief Compute the temperature based on gas properties.
*
* The temperature returned is consistent with the cooling rates.
*
* @param phys_const #phys_const data structure.
* @param cosmo #cosmology data structure.
* @param cooling #cooling_function_data struct.
* @param rho_phys Density of the gas in internal physical units.
* @param logZZsol Logarithm base 10 of the gas' metallicity in units of solar
* metallicity.
* @param XH The Hydrogen abundance of the gas.
* @param u_phys Internal energy of the gas in internal physical units.
*/
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) {
/* Convert to CGS */
const double u_cgs = u_phys * cooling->internal_energy_to_cgs;
const double n_H = rho_phys * XH / phys_const->const_proton_mass;
const double n_H_cgs = n_H * cooling->number_density_to_cgs;
/* 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;
get_index_1d(cooling->Redshifts, colibre_cooling_N_redshifts, cosmo->z,
&red_index, &d_red);
get_index_1d(cooling->Metallicity, colibre_cooling_N_metallicity, logZZsol,
&met_index, &d_met);
get_index_1d(cooling->nH, colibre_cooling_N_density, log10(n_H_cgs),
&n_H_index, &d_n_H);
/* Compute the log10 of the temperature by interpolating the table */
const double log_10_T =
colibre_convert_u_to_temp(log10(u_cgs), cosmo->z, n_H_index, d_n_H,
met_index, d_met, red_index, d_red, cooling);
/* Undo the log! */
return exp10(log_10_T);
}
/**
* @brief Compute the temperature of a #part based on the cooling function.
*
......@@ -184,8 +230,7 @@ float cooling_get_temperature(const struct phys_const *phys_const,
#endif
/* Get physical internal energy */
const float u = hydro_get_physical_internal_energy(p, xp, cosmo);
const double u_cgs = u * cooling->internal_energy_to_cgs;
const float u_phys = hydro_get_physical_internal_energy(p, xp, cosmo);
/* Get the Hydrogen mass fraction */
float const *metal_fraction =
......@@ -193,9 +238,7 @@ float cooling_get_temperature(const struct phys_const *phys_const,
const float XH = metal_fraction[chemistry_element_H];
/* 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;
const double n_H_cgs = n_H * cooling->number_density_to_cgs;
const float rho_phys = hydro_get_physical_density(p, cosmo);
/* Get this particle's metallicity ratio to solar.
*
......@@ -204,26 +247,8 @@ 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);
/* 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;
get_index_1d(cooling->Redshifts, colibre_cooling_N_redshifts, cosmo->z,
&red_index, &d_red);
get_index_1d(cooling->Metallicity, colibre_cooling_N_metallicity, logZZsol,
&met_index, &d_met);
get_index_1d(cooling->nH, colibre_cooling_N_density, log10(n_H_cgs),
&n_H_index, &d_n_H);
/* Compute the log10 of the temperature by interpolating the table */
const double log_10_T =
colibre_convert_u_to_temp(log10(u_cgs), cosmo->z, n_H_index, d_n_H,
met_index, d_met, red_index, d_red, cooling);
/* Undo the log! */
return exp10(log_10_T);
return cooling_get_temperature_from_gas(phys_const, cosmo, cooling, rho_phys,
logZZsol, XH, u_phys);
}
/**
......@@ -692,11 +717,6 @@ float cooling_get_particle_subgrid_HI_fraction(
const float log10_T_EOS_max =
log10f(max(T_EOS, FLT_MIN)) + cooling->dlogT_EOS;
/* Get the 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);
......@@ -712,6 +732,14 @@ float cooling_get_particle_subgrid_HI_fraction(
/* 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 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);
......@@ -745,11 +773,6 @@ float cooling_get_particle_subgrid_HII_fraction(
const float log10_T_EOS_max =
log10f(max(T_EOS, FLT_MIN)) + cooling->dlogT_EOS;
/* Get the 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);
......@@ -765,6 +788,14 @@ float cooling_get_particle_subgrid_HII_fraction(
/* 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 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);
......@@ -798,11 +829,6 @@ float cooling_get_particle_subgrid_H2_fraction(
const float log10_T_EOS_max =
log10f(max(T_EOS, FLT_MIN)) + cooling->dlogT_EOS;
/* Get the 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);
......@@ -818,6 +844,14 @@ float cooling_get_particle_subgrid_H2_fraction(
/* 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 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);
......@@ -850,11 +884,6 @@ float cooling_get_particle_subgrid_temperature(
const float log10_T_EOS_max =
log10f(max(T_EOS, FLT_MIN)) + cooling->dlogT_EOS;
/* Get the 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);
......@@ -870,6 +899,14 @@ float cooling_get_particle_subgrid_temperature(
/* 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 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);
......@@ -904,11 +941,6 @@ float cooling_get_particle_subgrid_density(
const float log10_T_EOS_max =
log10f(max(T_EOS, FLT_MIN)) + cooling->dlogT_EOS;
/* Get the 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);
......@@ -924,6 +956,14 @@ 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 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);
......@@ -953,11 +993,6 @@ void cooling_set_particle_subgrid_properties(
const float log10_T_EOS_max =
log10f(max(T_EOS, FLT_MIN)) + cooling->dlogT_EOS;
/* Get the 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);
......@@ -973,6 +1008,14 @@ void cooling_set_particle_subgrid_properties(
/* 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 log10_T = log10f(T);
p->cooling_data.subgrid_temp = compute_subgrid_temperature(
cooling, phys_const, floor_props, cosmo, rho_phys, logZZsol, XH, P_phys,
log10_T, log10_T_EOS_max);
......
......@@ -61,6 +61,11 @@ void cooling_first_init_part(const struct phys_const *phys_const,
const struct cooling_function_data *cooling,
struct part *p, struct xpart *xp);
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);
float cooling_get_temperature(const struct phys_const *phys_const,
const struct hydro_props *hydro_props,
const struct unit_system *us,
......@@ -110,6 +115,14 @@ void cooling_set_particle_subgrid_properties(
const struct cooling_function_data *cooling, struct part *p,
struct xpart *xp);
double compute_subgrid_density(
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);
float cooling_get_radiated_energy(const struct xpart *xp);
void cooling_split_part(struct part *p, struct xpart *xp, double n);
......
......@@ -568,6 +568,29 @@ __attribute__((always_inline)) INLINE void cooling_first_init_part(
xp->cooling_data.radiated_energy = 0.f;
}
/**
* @brief Compute the temperature based on gas properties.
*
* Dummy function!
*
* @param phys_const #phys_const data structure.
* @param cosmo #cosmology data structure.
* @param cooling #cooling_function_data struct.
* @param rho_phys Density of the gas in internal physical units.
* @param logZZsol Logarithm base 10 of the gas' metallicity in units of solar
* metallicity.
* @param XH The Hydrogen abundance of the gas.
* @param u_phys Internal energy of the gas in internal physical units.
*/