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

Moved the feedback related functions to a separate directory.

parent 98c285d0
......@@ -769,7 +769,8 @@ int main(int argc, char *argv[]) {
if (myrank == 0) chemistry_print(&chemistry);
/* Initialise stellar evolution */
if (with_feedback) stars_evolve_init(params, &stars_properties);
// MATTHIEU
// if (with_feedback) stars_evolve_init(params, &stars_properties);
/* Be verbose about what happens next */
if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
......
This diff is collapsed.
/**
* @brief Density interaction between two particles (non-symmetric).
* MATTHIEU: check with RGB about concerns with comment (in this and other
* subgrid schemes)
*
* @param r2 Comoving square distance between the two particles.
* @param dx Comoving vector separating both particles (pi - pj).
* @param hi Comoving smoothing-length of particle i.
* @param hj Comoving smoothing-length of particle j.
* @param si First sparticle.
* @param pj Second particle (not updated).
* @param a Current scale factor.
* @param H Current Hubble parameter.
* @param xp Extra particle data
* @param ti_current Current integer time value
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_stars_density(
float r2, const float *dx, float hi, float hj, struct spart *restrict si,
const struct part *restrict pj, const struct cosmology *restrict cosmo,
const struct stars_props *restrict stars_properties,
struct xpart *restrict xp, integertime_t ti_current) {
/* Get the gas mass. */
const float mj = hydro_get_mass(pj);
float wi, wi_dx;
/* Get r and 1/r. */
const float r_inv = 1.0f / sqrtf(r2);
const float r = r2 * r_inv;
/* Compute the kernel function */
const float hi_inv = 1.0f / hi;
const float ui = r * hi_inv;
kernel_deval(ui, &wi, &wi_dx);
float wj, wj_dx;
const float hj_inv = 1.0f / hj;
const float uj = r * hj_inv;
kernel_deval(uj, &wj, &wj_dx);
/* Add mass of pj to neighbour mass of si */
si->ngb_mass += hydro_get_mass(pj);
/* Add contribution of pj to normalisation of density weighted fraction
* which determines how much mass to distribute to neighbouring
* gas particles */
const float rho = hydro_get_comoving_density(pj);
if (rho != 0) si->density_weighted_frac_normalisation_inv += wj / rho;
/* Compute contribution to the density */
si->rho_gas += mj * wi;
}
/**
* @brief Feedback interaction between two particles (non-symmetric).
* Used for updating properties of gas particles neighbouring a star particle
*
* @param r2 Comoving square distance between the two particles.
* @param dx Comoving vector separating both particles (si - pj).
* @param hi Comoving smoothing-length of particle i.
* @param hj Comoving smoothing-length of particle j.
* @param si First (star) particle.
* @param pj Second (gas) particle.
* @param a Current scale factor.
* @param H Current Hubble parameter.
* @param xp Extra particle data
* @param ti_current Current integer time used value for seeding random number
* generator
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_stars_feedback(
float r2, const float *dx, float hi, float hj,
const struct spart *restrict si, struct part *restrict pj,
const struct cosmology *restrict cosmo,
const struct stars_props *restrict stars_properties,
struct xpart *restrict xp, integertime_t ti_current) {
// MATTHIEU
/* float wj; */
/* /\* Get r and 1/r. *\/ */
/* const float r_inv = 1.0f / sqrtf(r2); */
/* const float r = r2 * r_inv; */
/* /\* Compute the kernel function *\/ */
/* const float hj_inv = 1.0f / hj; */
/* const float uj = r * hj_inv; */
/* kernel_eval(uj, &wj); */
/* /\* Compute weighting for distributing feedback quantities *\/ */
/* float density_weighted_frac; */
/* float rho = hydro_get_comoving_density(pj); */
/* if (rho * si->density_weighted_frac_normalisation_inv != 0) { */
/* density_weighted_frac = */
/* wj / (rho * si->density_weighted_frac_normalisation_inv); */
/* } else { */
/* density_weighted_frac = 0.f; */
/* } */
/* /\* Update particle mass *\/ */
/* const float current_mass = hydro_get_mass(pj); */
/* float new_mass = */
/* current_mass + si->to_distribute.mass * density_weighted_frac; */
/* // for testing energy injection */
/* // new_mass = current_mass; */
/* hydro_set_mass(pj, new_mass); */
/* /\* Update total metallicity *\/ */
/* const float current_metal_mass_total = */
/* pj->chemistry_data.metal_mass_fraction_total * current_mass; */
/* const float new_metal_mass_total = */
/* current_metal_mass_total + */
/* si->to_distribute.total_metal_mass * density_weighted_frac; */
/* pj->chemistry_data.metal_mass_fraction_total = */
/* new_metal_mass_total / new_mass; */
/* /\* Update mass fraction of each tracked element *\/ */
/* for (int elem = 0; elem < chemistry_element_count; elem++) { */
/* const float current_metal_mass = */
/* pj->chemistry_data.metal_mass_fraction[elem] * current_mass; */
/* const float new_metal_mass = */
/* current_metal_mass + */
/* si->to_distribute.metal_mass[elem] * density_weighted_frac; */
/* pj->chemistry_data.metal_mass_fraction[elem] = new_metal_mass / new_mass;
*/
/* } */
/* /\* Update iron mass fraction from SNIa *\/ */
/* const float current_iron_from_SNIa_mass = */
/* pj->chemistry_data.iron_mass_fraction_from_SNIa * current_mass; */
/* const float new_iron_from_SNIa_mass = */
/* current_iron_from_SNIa_mass + */
/* si->to_distribute.Fe_mass_from_SNIa * density_weighted_frac; */
/* pj->chemistry_data.iron_mass_fraction_from_SNIa = */
/* new_iron_from_SNIa_mass / new_mass; */
/* /\* Update mass fraction from SNIa *\/ */
/* const float current_mass_from_SNIa = */
/* pj->chemistry_data.mass_from_SNIa * current_mass; */
/* const float new_mass_from_SNIa = */
/* current_mass_from_SNIa + */
/* si->to_distribute.mass_from_SNIa * density_weighted_frac; */
/* pj->chemistry_data.mass_from_SNIa = new_mass_from_SNIa / new_mass; */
/* /\* Update metal mass fraction from SNIa *\/ */
/* const float current_metal_mass_fraction_from_SNIa = */
/* pj->chemistry_data.metal_mass_fraction_from_SNIa * current_mass; */
/* const float new_metal_mass_fraction_from_SNIa = */
/* current_metal_mass_fraction_from_SNIa + */
/* si->to_distribute.metal_mass_from_SNIa * density_weighted_frac; */
/* pj->chemistry_data.metal_mass_fraction_from_SNIa = */
/* new_metal_mass_fraction_from_SNIa / new_mass; */
/* /\* Update mass fraction from SNII *\/ */
/* const float current_mass_from_SNII = */
/* pj->chemistry_data.mass_from_SNII * current_mass; */
/* const float new_mass_from_SNII = */
/* current_mass_from_SNII + */
/* si->to_distribute.mass_from_SNII * density_weighted_frac; */
/* pj->chemistry_data.mass_from_SNII = new_mass_from_SNII / new_mass; */
/* /\* Update metal mass fraction from SNII *\/ */
/* const float current_metal_mass_fraction_from_SNII = */
/* pj->chemistry_data.metal_mass_fraction_from_SNII * current_mass; */
/* const float new_metal_mass_fraction_from_SNII = */
/* current_metal_mass_fraction_from_SNII + */
/* si->to_distribute.metal_mass_from_SNII * density_weighted_frac; */
/* pj->chemistry_data.metal_mass_fraction_from_SNII = */
/* new_metal_mass_fraction_from_SNII / new_mass; */
/* /\* Update mass fraction from AGB *\/ */
/* const float current_mass_from_AGB = */
/* pj->chemistry_data.mass_from_AGB * current_mass; */
/* const float new_mass_from_AGB = */
/* current_mass_from_AGB + */
/* si->to_distribute.mass_from_AGB * density_weighted_frac; */
/* pj->chemistry_data.mass_from_AGB = new_mass_from_AGB / new_mass; */
/* /\* Update metal mass fraction from AGB *\/ */
/* const float current_metal_mass_fraction_from_AGB = */
/* pj->chemistry_data.metal_mass_fraction_from_AGB * current_mass; */
/* const float new_metal_mass_fraction_from_AGB = */
/* current_metal_mass_fraction_from_AGB + */
/* si->to_distribute.metal_mass_from_AGB * density_weighted_frac; */
/* pj->chemistry_data.metal_mass_fraction_from_AGB = */
/* new_metal_mass_fraction_from_AGB / new_mass; */
/* /\* Update momentum *\/ */
/* for (int i = 0; i < 3; i++) { */
/* pj->v[i] += */
/* si->to_distribute.mass * density_weighted_frac * (si->v[i] -
* pj->v[i]); */
/* } */
/* /\* Energy feedback *\/ */
/* float u_init = hydro_get_physical_internal_energy(pj, xp, cosmo); */
/* float heating_probability = -1.f, du = 0.f, d_energy = 0.f; */
/* d_energy += si->to_distribute.d_energy * density_weighted_frac; */
/* if (stars_properties->feedback.continuous_heating) { */
/* // We're doing ONLY continuous heating */
/* d_energy += si->to_distribute.num_SNIa * */
/* stars_properties->feedback.total_energy_SNe * */
/* density_weighted_frac * si->mass_init; */
/* } else { */
/* // We're doing stochastic heating */
/* heating_probability = si->to_distribute.heating_probability; */
/* du = stars_properties->feedback.SNe_deltaT_desired * */
/* stars_properties->feedback.temp_to_u_factor; */
/* if (heating_probability >= 1) { */
/* du = stars_properties->feedback.total_energy_SNe * */
/* si->to_distribute.num_SNIa / si->ngb_mass; */
/* heating_probability = 1; */
/* } */
/* double random_num = random_unit_interval(pj->id, ti_current, */
/* random_number_stellar_feedback);
*/
/* if (random_num < heating_probability) { */
/* message( */
/* "we did some heating! id %llu star id %llu probability %.5e " */
/* "random_num %.5e du %.5e " */
/* "du/ini %.5e", */
/* pj->id, si->id, heating_probability, random_num, du, */
/* du / hydro_get_physical_internal_energy(pj, xp, cosmo)); */
/* hydro_set_physical_internal_energy(pj, xp, cosmo, u_init + du); */
/* hydro_set_drifted_physical_internal_energy(pj, cosmo, u_init + du); */
/* } */
/* } */
/* /\* Add contribution from thermal and kinetic energy of ejected material
* (and */
/* * continuous SNIa feedback) *\/ */
/* u_init = hydro_get_physical_internal_energy(pj, xp, cosmo); */
/* du = d_energy / hydro_get_mass(pj); */
/* hydro_set_physical_internal_energy(pj, xp, cosmo, u_init + du); */
/* hydro_set_drifted_physical_internal_energy(pj, cosmo, u_init + du); */
}
/**
* @brief Stores AGB and SNII yield tables
*/
struct yield_table {
/* Yield table mass bins */
double *mass;
/* Yield table metallicity bins */
double *metallicity;
/* Array to store yield table resampled by IMF mass bins */
double *yield_IMF_resampled;
/* Array to store yield table being read in */
double *yield;
/* Array to store table of ejecta resampled by IMF mass bins */
double *ejecta_IMF_resampled;
/* Array to store table of ejecta being read in */
double *ejecta;
/* Array to store table of total mass released resampled by IMF mass bins */
double *total_metals_IMF_resampled;
/* Array to store table of total mass released being read in */
double *total_metals;
};
/**
* @brief Stores tables to determine stellar lifetimes. Used for calculation of
* IMF
*/
struct lifetime_table {
/* number of elements, mass, and initial metallicity bins */
int n_mass;
int n_z;
/* table of masses */
double *mass;
/* table of metallicities */
double *metallicity;
/* table of lifetimes depending on mass an metallicity */
double **dyingtime;
};
struct feedback_properties {
/* Flag to switch between continuous and stochastic heating */
int continuous_heating;
/* Desired temperature increase due to supernovae */
float SNe_deltaT_desired;
/* Conversion factor from temperature to internal energy */
float temp_to_u_factor;
/* Energy released by one supernova */
float total_energy_SNe;
/* Kinetic energy of SN ejecta per unit mass (check name with Richard)*/
float ejecta_specific_thermal_energy;
/* Solar mass */
float const_solar_mass;
/* Flag for testing energy injection */
int const_feedback_energy_testing;
/* Yield tables for AGB and SNII */
struct yield_table yield_AGB;
struct yield_table yield_SNII;
/* Array of adjustment factors for SNII */
double *typeII_factor;
/* Arrays of yield tables for SNIa */
double *yield_SNIa_IMF_resampled;
double yield_SNIa_total_metals_IMF_resampled;
double *yields_SNIa;
/* Parameters to SNIa enrichment model */
int SNIa_mode;
float SNIa_efficiency;
float SNIa_timescale_Gyr;
/* Arrays for names of elements being tracked for each enrichment channel */
char **SNIa_element_names;
char **SNII_element_names;
char **AGB_element_names;
/* Element name string length */
int element_name_length;
/* Sizes of dimensions of arrays in yield tables for each enrichment channel
*/
int SNIa_n_elements;
int SNII_n_mass;
int SNII_n_elements;
int SNII_n_z;
int AGB_n_mass;
int AGB_n_elements;
int AGB_n_z;
/* log10 of max and min allowable masses for SNII and SNIa in msun */
float log10_SNII_min_mass_msun;
float log10_SNII_max_mass_msun;
float log10_SNIa_max_mass_msun;
/* Array of mass bins for yield calculations */
double *yield_mass_bins;
/* Parameters for IMF */
/* IMF model name */
char IMF_Model[10];
/* Exponent for IMF if using power law */
float IMF_Exponent;
/* Array to store calculated IMF */
float *imf;
/* Arrays to store IMF mass bins */
float *imf_mass_bin;
float *imf_mass_bin_log10;
/* Number of IMF mass bins, maximum and minimum bins */
int n_imf_mass_bins;
float imf_max_mass_msun;
float imf_min_mass_msun;
float log10_imf_min_mass_msun;
float log10_imf_max_mass_msun;
/* Table of lifetime values */
struct lifetime_table lifetimes;
/* Location of yield tables */
char yield_table_path[50];
/* number of type II supernovae per solar mass */
float num_SNII_per_msun;
/* wind delay time for SNII */
float SNII_wind_delay;
/* Value to set birth time of stars read from ICs */
float spart_first_init_birth_time;
};
/**
* @brief Initialize the global properties of the feedback scheme.
*
* By default, takes the values provided by the hydro.
*
* @param sp The #feedback_properties.
* @param phys_const The physical constants in the internal unit system.
* @param us The internal unit system.
* @param params The parsed parameters.
* @param p The already read-in properties of the hydro scheme.
*/
INLINE static void feedbakc_props_init(struct feedback_props *fp,
const struct phys_const *phys_const,
const struct unit_system *us,
struct swift_params *params,
const struct hydro_props *p,
const struct cosmology *cosmo) {
/* Read SNIa timscale */
fp->SNIa_timescale_Gyr =
parser_get_param_float(params, "EAGLEFeedback:SNIa_timescale_Gyr");
/* Read the efficiency of producing SNIa */
fp->SNIa_efficiency =
parser_get_param_float(params, "EAGLEFeedback:SNIa_efficiency");
/* Are we doing continuous heating? */
fp->continuous_heating =
parser_get_param_int(params, "EAGLEFeedback:continuous_heating_switch");
/* Set the delay time before SNII occur */
const float Gyr_in_cgs = 1e9 * 365 * 24 * 3600;
fp->SNII_wind_delay =
parser_get_param_float(params, "EAGLEFeedback:SNII_wind_delay_Gyr") *
Gyr_in_cgs / units_cgs_conversion_factor(us, UNIT_CONV_TIME);
/* Read the temperature change to use in stochastic heating */
fp->SNe_deltaT_desired =
parser_get_param_float(params, "EAGLEFeedback:SNe_heating_temperature_K");
fp->SNe_deltaT_desired /=
units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
/* Set ejecta thermal energy */
const float ejecta_velocity =
1.0e6 / units_cgs_conversion_factor(
us, UNIT_CONV_SPEED); // EAGLE parameter is 10 km/s
fp->ejecta_specific_thermal_energy = 0.5 * ejecta_velocity * ejecta_velocity;
/* Energy released by supernova */
fp->total_energy_SNe =
1.0e51 / units_cgs_conversion_factor(us, UNIT_CONV_ENERGY);
/* Calculate temperature to internal energy conversion factor */
fp->temp_to_u_factor =
phys_const->const_boltzmann_k /
(p->mu_ionised * (hydro_gamma_minus_one)*phys_const->const_proton_mass);
/* Read birth time to set all stars in ICs to (defaults to -1 to indicate star
* present in ICs) */
fp->spart_first_init_birth_time = parser_get_opt_param_float(
params, "EAGLEFeedback:birth_time_override", -1);
/* Copy over solar mass */
fp->const_solar_mass = phys_const->const_solar_mass;
}
......@@ -206,7 +206,7 @@ inline static void init_imf(struct stars_props *restrict star_properties) {
const float imf_log10_mass_bin_size =
(star_properties->feedback.log10_imf_max_mass_msun -
star_properties->feedback.log10_imf_min_mass_msun) /
(float)(star_properties->feedback.n_imf_mass_bins - 1);
(float)(star_properties->feedback.n_imf_mass_bins - 1);
/* Allocate IMF array */
if (swift_memalign(
......@@ -233,9 +233,10 @@ inline static void init_imf(struct stars_props *restrict star_properties) {
if (strcmp(star_properties->feedback.IMF_Model, "Chabrier") == 0) {
for (int i = 0; i < star_properties->feedback.n_imf_mass_bins; i++) {
const float log10_mass_msun = star_properties->feedback.log10_imf_min_mass_msun +
i * imf_log10_mass_bin_size;
const float log10_mass_msun =
star_properties->feedback.log10_imf_min_mass_msun +
i * imf_log10_mass_bin_size;
const float mass_msun = exp10f(log10_mass_msun);
if (mass_msun > 1.0)
......@@ -256,10 +257,11 @@ inline static void init_imf(struct stars_props *restrict star_properties) {
}
/* Normalize the IMF */
const float norm = integrate_imf(star_properties->feedback.log10_imf_min_mass_msun,
star_properties->feedback.log10_imf_max_mass_msun,
eagle_imf_integration_mass_weight,
/*(stellar_yields=)*/ NULL, star_properties);
const float norm =
integrate_imf(star_properties->feedback.log10_imf_min_mass_msun,
star_properties->feedback.log10_imf_max_mass_msun,
eagle_imf_integration_mass_weight,
/*(stellar_yields=)*/ NULL, star_properties);
for (int i = 0; i < star_properties->feedback.n_imf_mass_bins; i++)
star_properties->feedback.imf[i] /= norm;
......@@ -275,15 +277,16 @@ inline static void init_imf(struct stars_props *restrict star_properties) {
* @param metallicity Star's metallicity
* @param star_properties the #stars_props data structure
*/
inline static double dying_mass_msun(const float age_Gyr, const float metallicity,
inline static double dying_mass_msun(const float age_Gyr,
const float metallicity,
const struct stars_props *star_props) {
float metallicity_offset, time_offset_low_metallicity = 0,
time_offset_high_metallicity = 0, mass_low_metallicity,
mass_high_metallicity;
time_offset_high_metallicity = 0,
mass_low_metallicity, mass_high_metallicity;
int metallicity_index, time_index_low_metallicity = -1,
time_index_high_metallicity = -1;
time_index_high_metallicity = -1;
if (age_Gyr <= 0) {
return star_props->feedback.imf_max_mass_msun;
......@@ -303,15 +306,14 @@ inline static double dying_mass_msun(const float age_Gyr, const float metallicit
for (metallicity_index = 0;
metallicity_index < star_props->feedback.lifetimes.n_z - 1;
metallicity_index++)
if (star_props->feedback.lifetimes
.metallicity[metallicity_index + 1] > metallicity)
if (star_props->feedback.lifetimes.metallicity[metallicity_index + 1] >
metallicity)
break;
metallicity_offset =
(metallicity -
star_props->feedback.lifetimes.metallicity[metallicity_index]) /
(star_props->feedback.lifetimes
.metallicity[metallicity_index + 1] -
(star_props->feedback.lifetimes.metallicity[metallicity_index + 1] -
star_props->feedback.lifetimes.metallicity[metallicity_index]);
}
......@@ -335,12 +337,10 @@ inline static double dying_mass_msun(const float age_Gyr, const float metallicit
star_props->feedback.lifetimes
.dyingtime[metallicity_index + 1]
[star_props->feedback.lifetimes.n_mass - 1]) {
time_index_high_metallicity =
star_props->feedback.lifetimes.n_mass - 2;
time_index_high_metallicity = star_props->feedback.lifetimes.n_mass - 2;
time_offset_high_metallicity = 1.0;
}
int i = star_props->feedback.lifetimes.n_mass;
while (i >= 0 && (time_index_low_metallicity == -1 ||
time_index_high_metallicity == -1)) {
......@@ -358,8 +358,8 @@ inline static double dying_mass_msun(const float age_Gyr, const float metallicit
star_props->feedback.lifetimes
.dyingtime[metallicity_index][time_index_low_metallicity]);
}
if (star_props->feedback.lifetimes
.dyingtime[metallicity_index + 1][i] >= log10_age_yr &&
if (star_props->feedback.lifetimes.dyingtime[metallicity_index + 1][i] >=
log10_age_yr &&
time_index_high_metallicity == -1) {
time_index_high_metallicity = i;
time_offset_high_metallicity =
......@@ -382,7 +382,7 @@ inline static double dying_mass_msun(const float age_Gyr, const float metallicit
time_index_high_metallicity, time_offset_high_metallicity);
float mass = (1.0 - metallicity_offset) * mass_low_metallicity +
metallicity_offset * mass_high_metallicity;
metallicity_offset * mass_high_metallicity;
/* Check that we haven't killed too many stars */
if (mass > star_props->feedback.imf_max_mass_msun)
......
......@@ -136,8 +136,8 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
struct spart *restrict sparts = c->stars.parts;
const struct engine *e = r->e;
const struct unit_system *us = e->internal_units;
const int with_cosmology = (e->policy & engine_policy_cosmology);
// const struct unit_system *us = e->internal_units;
// const int with_cosmology = (e->policy & engine_policy_cosmology);
const struct cosmology *cosmo = e->cosmology;
const float stars_h_max = e->hydro_properties->h_max;
const float stars_h_min = e->hydro_properties->h_min;
......@@ -204,16 +204,17 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
struct spart *sp = &sparts[sid[i]];
/* get particle timestep */
double dt;
if (with_cosmology) {
const integertime_t ti_step = get_integer_timestep(sp->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(e->ti_current - 1, sp->time_bin);
dt = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
ti_begin + ti_step);
} else {
dt = get_timestep(sp->time_bin, e->time_base);
}
/* double dt; */
/* if (with_cosmology) { */
/* const integertime_t ti_step = get_integer_timestep(sp->time_bin);
*/
/* const integertime_t ti_begin = */
/* get_integer_time_begin(e->ti_current - 1, sp->time_bin); */
/* dt = cosmology_get_therm_kick_factor(e->cosmology, ti_begin, */
/* ti_begin + ti_step); */