Skip to content
Snippets Groups Projects
Commit 36c89901 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Can run the isolated galaxy with SNII feedback on.

parent 36656159
No related branches found
No related tags found
1 merge request!787Eagle stellar evolution matthieu
Hide whitespace changes
Inline Side-by-side
Showing
with 440 additions and 310 deletions
examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml
+ 6
− 8
View file @ 36c89901
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.9891E43 # 10^10 solar masses
UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun in grams
UnitLength_in_cgs: 3.08567758E21 # 1 kpc
UnitVelocity_in_cgs: 1E5 # km/s
UnitCurrent_in_cgs: 1 # Amperes
......@@ -45,6 +45,7 @@ SPH:
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
h_min_ratio: 0.1 # Minimal smoothing in units of softening.
h_max: 10.
minimal_temperature: 100.
# Standard EAGLE cooling options
EAGLECooling:
......@@ -110,10 +111,7 @@ EAGLEEntropyFloor:
Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
EAGLEFeedback:
filename: ./yieldtables/
imf_model: Chabrier
continuous_heating_switch: 0
SNIa_timescale_Gyr: 2.0
SNIa_efficiency: 2.e-3
SNII_wind_delay_Gyr: 0.03
SNe_heating_temperature_K: 3.16228e7
filename: ./yieldtables/
SNII_wind_delay_Gyr: 0.0001
SNII_delta_T_K: 3.16228e7
SNII_Energy_erg: 1.0e51
src/cell.c
+ 2
− 0
View file @ 36c89901
......@@ -53,6 +53,7 @@
#include "drift.h"
#include "engine.h"
#include "error.h"
#include "feedback.h"
#include "gravity.h"
#include "hydro.h"
#include "hydro_properties.h"
......@@ -4366,6 +4367,7 @@ void cell_drift_spart(struct cell *c, const struct engine *e, int force) {
/* Get ready for a density calculation */
if (spart_is_active(sp, e)) {
stars_init_spart(sp);
feedback_init_spart(sp);
}
}
......
src/feedback/EAGLE/feedback.c
+ 232
− 128
View file @ 36c89901
......@@ -27,25 +27,113 @@
#include "yield_tables.h"
/**
* @brief Compute number of SN that should go off given the age of the spart
* @brief Return the change in temperature (in internal units) to apply to a
* gas particle affected by SNe feedback.
*
* @param sp spart we're evolving
* @param stars_properties stars_props data structure
* @param age age of star at the beginning of the step
* @param dt length of step
* @param sp The #spart.
* @param props The properties of the feedback model.
*/
float compute_SNe(const struct spart* sp,
const struct feedback_props* feedback_props, const float age,
const float dt) {
double eagle_feedback_temperature_change(const struct spart* sp,
const struct feedback_props* props) {
/* In the EAGLE REF model, the change of temperature is constant */
return props->SNe_deltaT_desired;
}
/**
* @brief Computes the number of super-novae exploding for a given
* star particle.
*
* @param sp The #spart.
* @param props The properties of the stellar model.
*/
double eagle_feedback_number_of_SNe(const struct spart* sp,
const struct feedback_props* props) {
message("init_mass=%e conv=%e N/Msun=%e", sp->mass_init,
props->mass_to_solar_mass, props->num_SNII_per_msun);
/* Note: For a Chabrier 2003 IMF and SNII going off between 6 and 100
* M_sun, the first term is 0.017362 M_sun^-1 */
return props->num_SNII_per_msun * sp->mass_init * props->mass_to_solar_mass;
}
/**
* @brief Computes the fraction of the available super-novae energy to
* inject for a given event.
*
* Note that the fraction can be > 1.
*
* @param sp The #spart.
* @param props The properties of the feedback model.
*/
double eagle_feedback_energy_fraction(const struct spart* sp,
const struct feedback_props* props) {
return 1.;
}
/**
* @brief Compute the properties of the SNe feedback energy injection.
*
* Only does something if the particle reached the SNe age during this time
* step.
*
* @param sp The star particle.
* @param feedback_props The properties of the feedback model.
* @param star_age Age of star at the beginning of the step in internal units.
* @param dt Length of time-step in internal units.
*/
inline static void compute_SNe_feedback(
struct spart* sp, const double star_age, const double dt,
const struct feedback_props* feedback_props) {
/* Time after birth considered for SNII feedback (internal units) */
const float SNII_wind_delay = feedback_props->SNII_wind_delay;
if (age <= SNII_wind_delay && (age + dt) > SNII_wind_delay) {
/* Are we doing feedback this step? */
if (star_age <= SNII_wind_delay && (star_age + dt) > SNII_wind_delay) {
return feedback_props->num_SNII_per_msun * sp->mass_init /
feedback_props->const_solar_mass;
} else {
return 0.f;
/* Properties of the model (all in internal units) */
const double delta_T =
eagle_feedback_temperature_change(sp, feedback_props);
const double N_SNe = eagle_feedback_number_of_SNe(sp, feedback_props);
const double E_SNe = feedback_props->E_SNII;
const double f_E = eagle_feedback_energy_fraction(sp, feedback_props);
/* Conversion factor from T to internal energy */
const double conv_factor = feedback_props->temp_to_u_factor;
/* Calculate the default heating probability */
double prob = f_E * E_SNe * N_SNe /
(conv_factor * delta_T * sp->feedback_data.ngb_mass);
/* Calculate the change in internal energy of the gas particles that get
* heated */
double delta_u;
if (prob <= 1.) {
/* Normal case */
delta_u = delta_T * conv_factor;
} else {
/* Special case: we need to adjust the energy irrespective of the
desired deltaT to ensure we inject all the available energy. */
prob = 1.;
delta_u = f_E * E_SNe * N_SNe / sp->feedback_data.ngb_mass;
}
message(
"ID=%lld E_SNe=%e N_SNe=%e deltaT= %e f_E=%e prob=%f ngb_mass=%f "
"fac=%e delta_u=%e",
sp->id, E_SNe, N_SNe, delta_T, f_E, prob, sp->feedback_data.ngb_mass,
conv_factor, delta_u);
/* Store all of this in the star for delivery onto the gas */
sp->feedback_data.to_distribute.SNII_heating_probability = prob;
sp->feedback_data.to_distribute.SNII_delta_u = delta_u;
}
}
......@@ -182,53 +270,59 @@ inline static void evolve_SNIa(float log10_min_mass, float log10_max_mass,
struct spart* restrict sp, float star_age_Gyr,
float dt_Gyr) {
/* Check if we're outside the mass range for SNIa */
if (log10_min_mass >= feedback_props->log10_SNIa_max_mass_msun) return;
/* If the max mass is outside the mass range update it to be the maximum and
* use updated values for the star's age and timestep in this function */
if (log10_max_mass > feedback_props->log10_SNIa_max_mass_msun) {
log10_max_mass = feedback_props->log10_SNIa_max_mass_msun;
float lifetime_Gyr = lifetime_in_Gyr(
exp(M_LN10 * feedback_props->log10_SNIa_max_mass_msun),
sp->chemistry_data.metal_mass_fraction_total, feedback_props);
dt_Gyr = star_age_Gyr + dt_Gyr - lifetime_Gyr;
star_age_Gyr = lifetime_Gyr;
}
/* /\* Check if we're outside the mass range for SNIa *\/ */
/* if (log10_min_mass >= feedback_props->log10_SNIa_max_mass_msun) return; */
/* compute the number of SNIa */
/* Efolding (Forster 2006) */
float num_SNIa_per_msun =
feedback_props->SNIa_efficiency *
(exp(-star_age_Gyr / feedback_props->SNIa_timescale_Gyr) -
exp(-(star_age_Gyr + dt_Gyr) / feedback_props->SNIa_timescale_Gyr)) *
sp->mass_init;
sp->feedback_data.to_distribute.num_SNIa =
num_SNIa_per_msun / feedback_props->const_solar_mass;
/* compute mass fractions of each metal */
for (int i = 0; i < chemistry_element_count; i++) {
sp->feedback_data.to_distribute.metal_mass[i] +=
num_SNIa_per_msun * feedback_props->yield_SNIa_IMF_resampled[i];
}
/* /\* If the max mass is outside the mass range update it to be the maximum
* and */
/* * use updated values for the star's age and timestep in this function *\/
*/
/* if (log10_max_mass > feedback_props->log10_SNIa_max_mass_msun) { */
/* log10_max_mass = feedback_props->log10_SNIa_max_mass_msun; */
/* float lifetime_Gyr = lifetime_in_Gyr( */
/* exp(M_LN10 * feedback_props->log10_SNIa_max_mass_msun), */
/* sp->chemistry_data.metal_mass_fraction_total, feedback_props); */
/* dt_Gyr = star_age_Gyr + dt_Gyr - lifetime_Gyr; */
/* star_age_Gyr = lifetime_Gyr; */
/* } */
/* Update the metallicity of the material released */
sp->feedback_data.to_distribute.metal_mass_from_SNIa +=
num_SNIa_per_msun * feedback_props->yield_SNIa_total_metals_IMF_resampled;
/* /\* compute the number of SNIa *\/ */
/* /\* Efolding (Forster 2006) *\/ */
/* float num_SNIa_per_msun = */
/* feedback_props->SNIa_efficiency * */
/* (exp(-star_age_Gyr / feedback_props->SNIa_timescale_Gyr) - */
/* exp(-(star_age_Gyr + dt_Gyr) / feedback_props->SNIa_timescale_Gyr)) *
*/
/* sp->mass_init; */
/* Update the metal mass produced */
sp->feedback_data.to_distribute.total_metal_mass +=
num_SNIa_per_msun * feedback_props->yield_SNIa_total_metals_IMF_resampled;
/* sp->feedback_data.to_distribute.num_SNIa = */
/* num_SNIa_per_msun / feedback_props->const_solar_mass; */
/* Compute the mass produced by SNIa */
sp->feedback_data.to_distribute.mass_from_SNIa +=
num_SNIa_per_msun * feedback_props->yield_SNIa_total_metals_IMF_resampled;
/* /\* compute mass fractions of each metal *\/ */
/* for (int i = 0; i < chemistry_element_count; i++) { */
/* sp->feedback_data.to_distribute.metal_mass[i] += */
/* num_SNIa_per_msun * feedback_props->yield_SNIa_IMF_resampled[i]; */
/* } */
/* Compute the iron mass produced */
sp->feedback_data.to_distribute.Fe_mass_from_SNIa +=
num_SNIa_per_msun *
feedback_props->yield_SNIa_IMF_resampled[chemistry_element_Fe];
/* /\* Update the metallicity of the material released *\/ */
/* sp->feedback_data.to_distribute.metal_mass_from_SNIa += */
/* num_SNIa_per_msun *
* feedback_props->yield_SNIa_total_metals_IMF_resampled; */
/* /\* Update the metal mass produced *\/ */
/* sp->feedback_data.to_distribute.total_metal_mass += */
/* num_SNIa_per_msun *
* feedback_props->yield_SNIa_total_metals_IMF_resampled; */
/* /\* Compute the mass produced by SNIa *\/ */
/* sp->feedback_data.to_distribute.mass_from_SNIa += */
/* num_SNIa_per_msun *
* feedback_props->yield_SNIa_total_metals_IMF_resampled; */
/* /\* Compute the iron mass produced *\/ */
/* sp->feedback_data.to_distribute.Fe_mass_from_SNIa += */
/* num_SNIa_per_msun * */
/* feedback_props->yield_SNIa_IMF_resampled[chemistry_element_Fe]; */
}
/**
......@@ -560,7 +654,7 @@ void compute_stellar_evolution(const struct feedback_props* feedback_props,
const float dt) {
/* Allocate temporary array for calculating imf weights */
float stellar_yields[eagle_feedback_N_imf_bins];
// float stellar_yields[eagle_feedback_N_imf_bins];
/* Convert dt and stellar age from internal units to Gyr. */
const double Gyr_in_cgs = 1e9 * 365. * 24. * 3600.;
......@@ -572,7 +666,10 @@ void compute_stellar_evolution(const struct feedback_props* feedback_props,
/* Get the metallicity */
const float Z = sp->chemistry_data.metal_mass_fraction_total;
/* calculate mass of stars that has died from the star's birth up to the
/* Compute properties of the stochastic SNe feedback model. */
compute_SNe_feedback(sp, age, dt, feedback_props);
/* Calculate mass of stars that has died from the star's birth up to the
* beginning and end of timestep */
const float max_dying_mass_Msun =
dying_mass_msun(star_age_Gyr, Z, feedback_props);
......@@ -592,42 +689,45 @@ void compute_stellar_evolution(const struct feedback_props* feedback_props,
if (min_dying_mass_Msun == max_dying_mass_Msun) return;
/* Life is better in log */
const float log10_max_dying_mass_Msun = log10f(max_dying_mass_Msun);
const float log10_min_dying_mass_Msun = log10f(min_dying_mass_Msun);
/* Compute elements, energy and momentum to distribute from the
* three channels SNIa, SNII, AGB */
evolve_SNIa(log10_min_dying_mass_Msun, log10_max_dying_mass_Msun,
feedback_props, sp, star_age_Gyr, dt_Gyr);
evolve_SNII(log10_min_dying_mass_Msun, log10_max_dying_mass_Msun,
stellar_yields, feedback_props, sp);
evolve_AGB(log10_min_dying_mass_Msun, log10_max_dying_mass_Msun,
stellar_yields, feedback_props, sp);
/* Compute the total mass to distribute (H + He metals) */
sp->feedback_data.to_distribute.mass =
sp->feedback_data.to_distribute.total_metal_mass +
sp->feedback_data.to_distribute.metal_mass[chemistry_element_H] +
sp->feedback_data.to_distribute.metal_mass[chemistry_element_He];
/* Compute the number of type II SNe that went off */
sp->feedback_data.to_distribute.num_SNe =
compute_SNe(sp, feedback_props, age, dt);
/* Compute energy change due to thermal and kinetic energy of ejecta */
sp->feedback_data.to_distribute.d_energy =
sp->feedback_data.to_distribute.mass *
(feedback_props->ejecta_specific_thermal_energy +
0.5 * (sp->v[0] * sp->v[0] + sp->v[1] * sp->v[1] + sp->v[2] * sp->v[2]) *
cosmo->a2_inv);
/* Compute probability of heating neighbouring particles */
if (dt > 0 && sp->feedback_data.ngb_mass > 0)
sp->feedback_data.to_distribute.heating_probability =
feedback_props->total_energy_SNe *
sp->feedback_data.to_distribute.num_SNe /
(feedback_props->temp_to_u_factor * feedback_props->SNe_deltaT_desired *
sp->feedback_data.ngb_mass);
/* const float log10_max_dying_mass_Msun = log10f(max_dying_mass_Msun); */
/* const float log10_min_dying_mass_Msun = log10f(min_dying_mass_Msun); */
/* /\* Compute elements, energy and momentum to distribute from the */
/* * three channels SNIa, SNII, AGB *\/ */
/* evolve_SNIa(log10_min_dying_mass_Msun, log10_max_dying_mass_Msun, */
/* feedback_props, sp, star_age_Gyr, dt_Gyr); */
/* evolve_SNII(log10_min_dying_mass_Msun, log10_max_dying_mass_Msun, */
/* stellar_yields, feedback_props, sp); */
/* evolve_AGB(log10_min_dying_mass_Msun, log10_max_dying_mass_Msun, */
/* stellar_yields, feedback_props, sp); */
/* /\* Compute the total mass to distribute (H + He metals) *\/ */
/* sp->feedback_data.to_distribute.mass = */
/* sp->feedback_data.to_distribute.total_metal_mass + */
/* sp->feedback_data.to_distribute.metal_mass[chemistry_element_H] + */
/* sp->feedback_data.to_distribute.metal_mass[chemistry_element_He]; */
/* /\* Compute energy change due to thermal and kinetic energy of ejecta *\/
*/
/* sp->feedback_data.to_distribute.d_energy = */
/* sp->feedback_data.to_distribute.mass * */
/* (feedback_props->ejecta_specific_thermal_energy + */
/* 0.5 * (sp->v[0] * sp->v[0] + sp->v[1] * sp->v[1] + sp->v[2] *
* sp->v[2]) * */
/* cosmo->a2_inv); */
/* /\* Compute the number of type II SNe that went off *\/ */
/* sp->feedback_data.to_distribute.num_SNe = */
/* compute_SNe(sp, feedback_props, age, dt); */
/* /\* Compute probability of heating neighbouring particles *\/ */
/* if (dt > 0 && sp->feedback_data.ngb_mass > 0) */
/* sp->feedback_data.to_distribute.heating_probability = */
/* feedback_props->total_energy_SNe * */
/* sp->feedback_data.to_distribute.num_SNe / */
/* (feedback_props->temp_to_u_factor *
* feedback_props->SNe_deltaT_desired * */
/* sp->feedback_data.ngb_mass); */
}
/**
......@@ -648,17 +748,13 @@ void feedback_props_init(struct feedback_props* fp,
const struct hydro_props* hydro_props,
const struct cosmology* cosmo) {
/* Read SNIa timscale */
fp->SNIa_timescale_Gyr =
parser_get_param_float(params, "EAGLEFeedback:SNIa_timescale_Gyr");
/* /\* 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");
/* /\* Read the efficiency of producing SNIa *\/ */
/* fp->SNIa_efficiency = */
/* parser_get_param_float(params, "EAGLEFeedback:SNIa_efficiency"); */
/* Set the delay time before SNII occur */
const float Gyr_in_cgs = 1e9 * 365 * 24 * 3600;
......@@ -666,36 +762,43 @@ void feedback_props_init(struct feedback_props* fp,
parser_get_param_float(params, "EAGLEFeedback:SNII_wind_delay_Gyr") *
Gyr_in_cgs / units_cgs_conversion_factor(us, UNIT_CONV_TIME);
message("Wind delay: %e", fp->SNII_wind_delay);
/* Read the temperature change to use in stochastic heating */
fp->SNe_deltaT_desired =
parser_get_param_float(params, "EAGLEFeedback:SNe_heating_temperature_K");
parser_get_param_float(params, "EAGLEFeedback:SNII_delta_T_K");
fp->SNe_deltaT_desired /=
units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
/* Energy released by supernova type II */
fp->E_SNII_cgs =
parser_get_param_double(params, "EAGLEFeedback:SNII_Energy_erg");
fp->E_SNII =
fp->E_SNII_cgs / units_cgs_conversion_factor(us, UNIT_CONV_ENERGY);
/* Gather common conversion factors */
/* Calculate internal mass to solar mass conversion factor */
const double Msun_cgs = phys_const->const_solar_mass *
units_cgs_conversion_factor(us, UNIT_CONV_MASS);
const double unit_mass_cgs = units_cgs_conversion_factor(us, UNIT_CONV_MASS);
fp->mass_to_solar_mass = unit_mass_cgs / Msun_cgs;
/* Calculate temperature to internal energy conversion factor (all internal
* units) */
const double k_B = phys_const->const_boltzmann_k;
const double m_p = phys_const->const_proton_mass;
const double mu = hydro_props->mu_ionised;
fp->temp_to_u_factor = k_B / (mu * hydro_gamma_minus_one * m_p);
// MATTHIEU up to here
/* 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 /
(hydro_props->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;
/* Set number of elements found in yield tables */
fp->lifetimes.n_mass = 30;
fp->lifetimes.n_z = 6;
......@@ -745,12 +848,13 @@ void feedback_props_init(struct feedback_props* fp,
* mass bins used in IMF */
compute_ejecta(fp);
/* Calculate number of type II SN per solar mass
* Note: since we are integrating the IMF without weighting it by the yields
* pass NULL pointer for stellar_yields array */
fp->num_SNII_per_msun = integrate_imf(fp->log10_SNII_min_mass_msun,
fp->log10_SNII_max_mass_msun, 0,
/*(stellar_yields=)*/ NULL, fp);
/* Calculate number of type II SN per unit solar mass based on our choice
* of IMF and integration limits for type II SNe.
* Note: No weighting by yields here. */
fp->num_SNII_per_msun =
integrate_imf(fp->log10_SNII_min_mass_msun, fp->log10_SNII_max_mass_msun,
eagle_imf_integration_no_weight,
/*(stellar_yields=)*/ NULL, fp);
message("initialized stellar feedback");
}
src/feedback/EAGLE/feedback.h
+ 6
− 2
View file @ 36c89901
......@@ -20,6 +20,7 @@
#define SWIFT_FEEDBACK_EAGLE_H
#include "cosmology.h"
#include "error.h"
#include "feedback_properties.h"
#include "hydro_properties.h"
#include "part.h"
......@@ -94,8 +95,11 @@ __attribute__((always_inline)) INLINE static void feedback_prepare_spart(
/* Zero the energy to inject */
sp->feedback_data.to_distribute.d_energy = 0.f;
/* Zero the feedback probability */
sp->feedback_data.to_distribute.heating_probability = 0.f;
/* Zero the SNII feedback probability */
sp->feedback_data.to_distribute.SNII_heating_probability = 0.f;
/* Zero the SNII feedback energy */
sp->feedback_data.to_distribute.SNII_delta_u = 0.f;
}
/**
......
src/feedback/EAGLE/feedback_iact.h
+ 156
− 146
View file @ 36c89901
......@@ -56,18 +56,14 @@ runner_iact_nonsym_feedback_density(float r2, const float *dx, float hi,
kernel_eval(ui, &wi);
/* Add mass of pj to neighbour mass of si */
si->feedback_data.ngb_mass += hydro_get_mass(pj);
si->feedback_data.ngb_mass += mj;
/* 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->feedback_data.density_weighted_frac_normalisation_inv += wi / rho;
/* Compute contribution to the density */
si->rho_gas += mj * wi;
// const float rho = hydro_get_comoving_density(pj);
// si->feedback_data.density_weighted_frac_normalisation_inv += wi / rho;
}
/**
......@@ -104,153 +100,167 @@ runner_iact_nonsym_feedback_apply(
float wi;
kernel_eval(ui, &wi);
/* Compute weighting for distributing feedback quantities */
float density_weighted_frac;
float rho = hydro_get_comoving_density(pj);
if (rho * si->feedback_data.density_weighted_frac_normalisation_inv != 0) {
density_weighted_frac =
wi / (rho * si->feedback_data.density_weighted_frac_normalisation_inv);
} else {
density_weighted_frac = 0.f;
}
/* /\* Compute weighting for distributing feedback quantities *\/ */
/* float density_weighted_frac; */
/* float rho = hydro_get_comoving_density(pj); */
/* if (rho * si->feedback_data.density_weighted_frac_normalisation_inv != 0) {
*/
/* density_weighted_frac = */
/* wi / (rho *
* si->feedback_data.density_weighted_frac_normalisation_inv); */
/* } else { */
/* density_weighted_frac = 0.f; */
/* } */
/* Update particle mass */
const float current_mass = hydro_get_mass(pj);
const float new_mass = current_mass + si->feedback_data.to_distribute.mass *
density_weighted_frac;
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->feedback_data.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->feedback_data.to_distribute.metal_mass[elem] *
density_weighted_frac;
pj->chemistry_data.metal_mass_fraction[elem] = new_metal_mass / new_mass;
}
/* /\* Update particle mass *\/ */
/* const float current_mass = hydro_get_mass(pj); */
/* const float new_mass = current_mass + si->feedback_data.to_distribute.mass
* * */
/* density_weighted_frac; */
/* 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->feedback_data.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->feedback_data.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->feedback_data.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->feedback_data.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->feedback_data.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->feedback_data.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->feedback_data.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->feedback_data.to_distribute.mass * density_weighted_frac *
(si->v[i] - pj->v[i]);
}
/* 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->feedback_data.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->feedback_data.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->feedback_data.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->feedback_data.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->feedback_data.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->feedback_data.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->feedback_data.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->feedback_data.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->feedback_data.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->feedback_data.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->feedback_data.to_distribute.d_energy * density_weighted_frac;
if (feedback_props->continuous_heating) {
// We're doing ONLY continuous heating
d_energy += si->feedback_data.to_distribute.num_SNIa *
feedback_props->total_energy_SNe * density_weighted_frac *
si->mass_init;
} else {
// We're doing stochastic heating
heating_probability = si->feedback_data.to_distribute.heating_probability;
du = feedback_props->SNe_deltaT_desired * feedback_props->temp_to_u_factor;
if (heating_probability >= 1) {
du = feedback_props->total_energy_SNe *
si->feedback_data.to_distribute.num_SNIa /
si->feedback_data.ngb_mass;
heating_probability = 1;
}
// d_energy += si->feedback_data.to_distribute.d_energy *
// density_weighted_frac;
/* if (feedback_props->continuous_heating) { */
/* // We're doing ONLY continuous heating */
/* d_energy += si->feedback_data.to_distribute.num_SNIa * */
/* feedback_props->total_energy_SNe * density_weighted_frac * */
/* si->mass_init; */
/* } else { */
/* /\* 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); */
/* Get the SNII feedback properties */
const float prob = si->feedback_data.to_distribute.SNII_heating_probability;
const float delta_u = si->feedback_data.to_distribute.SNII_delta_u;
/* Are we doing some SNII feedback? */
if (prob > 0.f) {
/* Draw a random number (Note mixing both IDs) */
const float rand = random_unit_interval(si->id + pj->id, ti_current,
random_number_stellar_feedback);
/* Are we lucky? */
if (rand < prob) {
/* Compute new energy of this particle */
const float u_init = hydro_get_physical_internal_energy(pj, xp, cosmo);
const float u_new = u_init + delta_u;
/* Inject energy into the particle */
hydro_set_physical_internal_energy(pj, xp, cosmo, u_new);
hydro_set_drifted_physical_internal_energy(pj, cosmo, u_new);
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);
"We did some heating! id %llu star id %llu probability %.5e "
"random_num %.5e du %.5e du/ini %.5e",
pj->id, si->id, prob, rand, delta_u, delta_u / u_init);
}
}
/* 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);
}
#endif /* SWIFT_EAGLE_FEEDBACK_IACT_H */
src/feedback/EAGLE/feedback_properties.h
+ 27
− 23
View file @ 36c89901
......@@ -73,24 +73,9 @@ struct lifetime_table {
struct feedback_props {
/* 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;
/* Yield tables for AGB and SNII */
struct yield_table yield_AGB;
struct yield_table yield_SNII;
......@@ -120,7 +105,21 @@ struct feedback_props {
/* Array of mass bins for yield calculations */
double *yield_mass_bins;
/* Parameters for IMF */
/* Table of lifetime values */
struct lifetime_table lifetimes;
/* Location of yield tables */
char yield_table_path[200];
/* ------------- Conversion factors --------------- */
/*! Conversion factor from internal mass unit to solar mass */
double mass_to_solar_mass;
/*! Conversion factor from temperature to internal energy */
float temp_to_u_factor;
/* ------------- Parameters for IMF --------------- */
/*! Array to store calculated IMF */
float *imf;
......@@ -137,17 +136,22 @@ struct feedback_props {
float log10_imf_min_mass_msun;
float log10_imf_max_mass_msun;
/* Table of lifetime values */
struct lifetime_table lifetimes;
/* ------------ SNe feedback properties ------------ */
/* Location of yield tables */
char yield_table_path[200];
/* number of type II supernovae per solar mass */
/*! Number of type II supernovae per solar mass */
float num_SNII_per_msun;
/* wind delay time for SNII */
/*! Wind delay time for SNII */
float SNII_wind_delay;
/*! Temperature increase induced by SNe feedback */
float SNe_deltaT_desired;
/* Energy released by one supernova type II in cgs units */
double E_SNII_cgs;
/* Energy released by one supernova type II in internal units */
float E_SNII;
};
void feedback_props_init(struct feedback_props *fp,
......
src/feedback/EAGLE/feedback_struct.h
+ 5
− 2
View file @ 36c89901
......@@ -65,8 +65,11 @@ struct feedback_spart_data {
/* Energy change due to thermal and kinetic energy of ejecta */
float d_energy;
/* Probability for heating neighbouring gas particles */
float heating_probability;
/*! Probability to heating neighbouring gas particle for SNII feedback */
float SNII_heating_probability;
/*! Change in energy from SNII feedback energy injection */
float SNII_delta_u;
} to_distribute;
......
src/feedback/EAGLE/imf.h
+ 3
− 1
View file @ 36c89901
......@@ -193,7 +193,7 @@ inline static float integrate_imf(const float log10_min_mass,
* table.
*
* @param star_properties #stars_props data structure */
inline static void init_imf(struct feedback_props *restrict feedback_props) {
inline static void init_imf(struct feedback_props *feedback_props) {
/* Define max and min imf masses */
feedback_props->imf_max_mass_msun = 100.f;
......@@ -280,6 +280,8 @@ inline static float dying_mass_msun(
const float age_Gyr, const float Z,
const struct feedback_props *feedback_props) {
// MATTHIEU check this!!!
/* Pull out some common terms */
const double *lifetime_Z = feedback_props->lifetimes.metallicity;
const double *lifetime_m = feedback_props->lifetimes.mass;
......
src/runner.c
+ 3
− 0
View file @ 36c89901
......@@ -326,6 +326,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
/* Re-initialise everything */
stars_init_spart(sp);
feedback_init_spart(sp);
/* Off we go ! */
continue;
......@@ -713,6 +714,8 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
/* Did we get a star? (Or did we run out of spare ones?) */
if (sp != NULL) {
message("Formed a star ID=%lld", sp->id);
/* Copy the properties of the gas particle to the star particle */
star_formation_copy_properties(p, xp, sp, e, sf_props, cosmo,
with_cosmology);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment