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

Added calculation of the energy to inject from the enrichment channels.

parent 0dd11a70
Branches
Tags
1 merge request!787Eagle stellar evolution matthieu
......@@ -78,33 +78,35 @@ EAGLEChemistry:
# Properties of the EAGLE feedback and enrichment model.
EAGLEFeedback:
use_SNe_feedback: 0 # No SNe feedback in this test.
use_AGB_enrichment: 1
use_SNII_enrichment: 1
use_SNIa_enrichment: 1
filename: ./yieldtables/
IMF_min_mass_Msun: 0.1
IMF_max_mass_Msun: 100.0
SNII_min_mass_Msun: 6.0
SNII_max_mass_Msun: 100.0
SNII_wind_delay_Gyr: 0.03
SNII_delta_T_K: 3.16228e7
SNII_Energy_erg: 1.0e51
SNII_Energy_fraction_min: 3.0
SNII_Energy_fraction_max: 0.3
SNII_Energy_fraction_Z_0: 0.0012663729
SNII_Energy_fraction_n_0_H_p_cm3: 0.67
SNII_Energy_fraction_n_n: 0.8686
SNII_Energy_fraction_n_Z: 0.8686
SNIa_max_mass_Msun: 8.0
SNIa_timescale_Gyr: 2.0
SNIa_efficiency_p_Msun: 0.002
SNII_yield_factor_Hydrogen: 1.0 # Correction to the yields following the EAGLE REF values
SNII_yield_factor_Helium: 1.0
SNII_yield_factor_Carbon: 0.5
SNII_yield_factor_Nitrogen: 1.0
SNII_yield_factor_Oxygen: 1.0
SNII_yield_factor_Neon: 1.0
SNII_yield_factor_Magnesium: 2.0
SNII_yield_factor_Silicon: 1.0
SNII_yield_factor_Iron: 0.5
use_SNe_feedback: 0 # Global switch for SNe thermal feedback.
use_AGB_enrichment: 1 # Global switch for enrichement from AGB stars.
use_SNII_enrichment: 1 # Global switch for enrichement from SNII stars.
use_SNIa_enrichment: 1 # Global switch for enrichement from SNIa stars.
filename: ./yieldtables/ # Path to the directory containing the EAGLE yield tables.
IMF_min_mass_Msun: 0.1 # Minimal stellar mass considered for the Chabrier IMF in solar masses.
IMF_max_mass_Msun: 100.0 # Maximal stellar mass considered for the Chabrier IMF in solar masses.
SNII_min_mass_Msun: 6.0 # Minimal mass considered for SNII feedback (not SNII enrichment!) in solar masses.
SNII_max_mass_Msun: 100.0 # Maximal mass considered for SNII feedback (not SNII enrichment!) in solar masses.
SNII_wind_delay_Gyr: 0.03 # Time in Gyr between a star's birth and the SNII thermal feedback event.
SNII_delta_T_K: 3.16228e7 # Change in temperature to apply to the gas particle in a SNII thermal feedback event in Kelvin.
SNII_Energy_erg: 1.0e51 # Energy of one SNII explosion in ergs.
SNII_Energy_fraction_min: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
SNII_Energy_fraction_max: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
SNII_Energy_fraction_Z_0: 0.0012663729 # Pivot point for the metallicity dependance of the SNII energy fraction (metal mass fraction).
SNII_Energy_fraction_n_0_H_p_cm3: 0.67 # Pivot point for the birth density dependance of the SNII energy fraction in cm^-3.
SNII_Energy_fraction_n_Z: 0.8686 # Power-law for the metallicity dependance of the SNII energy fraction.
SNII_Energy_fraction_n_n: 0.8686 # Power-law for the birth density dependance of the SNII energy fraction.
SNIa_max_mass_Msun: 8.0 # Maximal mass considered for SNIa feedback and enrichment in solar masses.
SNIa_timescale_Gyr: 2.0 # Time-scale of the exponential decay of the SNIa rates in Gyr.
SNIa_efficiency_p_Msun: 0.002 # Normalisation of the SNIa rates in inverse solar masses.
SNIa_Energy_erg: 1.0e51 # Energy of one SNIa explosion in ergs.
AGB_ejecta_velocity_km_p_s: 10.0 # Velocity of the AGB ejectas in km/s.
SNII_yield_factor_Hydrogen: 1.0 # (Optional) Correction factor to apply to the Hydrogen yield from the SNII channel.
SNII_yield_factor_Helium: 1.0 # (Optional) Correction factor to apply to the Helium yield from the SNII channel.
SNII_yield_factor_Carbon: 0.5 # (Optional) Correction factor to apply to the Carbon yield from the SNII channel.
SNII_yield_factor_Nitrogen: 1.0 # (Optional) Correction factor to apply to the Nitrogen yield from the SNII channel.
SNII_yield_factor_Oxygen: 1.0 # (Optional) Correction factor to apply to the Oxygen yield from the SNII channel.
SNII_yield_factor_Neon: 1.0 # (Optional) Correction factor to apply to the Neon yield from the SNII channel.
SNII_yield_factor_Magnesium: 2.0 # (Optional) Correction factor to apply to the Magnesium yield from the SNII channel.
SNII_yield_factor_Silicon: 1.0 # (Optional) Correction factor to apply to the Silicon yield from the SNII channel.
SNII_yield_factor_Iron: 0.5 # (Optional) Correction factor to apply to the Iron yield from the SNII channel.
......@@ -360,6 +360,8 @@ EAGLEFeedback:
SNIa_max_mass_Msun: 8.0 # Maximal mass considered for SNIa feedback and enrichment in solar masses.
SNIa_timescale_Gyr: 2.0 # Time-scale of the exponential decay of the SNIa rates in Gyr.
SNIa_efficiency_p_Msun: 0.002 # Normalisation of the SNIa rates in inverse solar masses.
SNIa_Energy_erg: 1.0e51 # Energy of one SNIa explosion in ergs.
AGB_ejecta_velocity_km_p_s: 10.0 # Velocity of the AGB ejectas in km/s.
SNII_yield_factor_Hydrogen: 1.0 # (Optional) Correction factor to apply to the Hydrogen yield from the SNII channel.
SNII_yield_factor_Helium: 1.0 # (Optional) Correction factor to apply to the Helium yield from the SNII channel.
SNII_yield_factor_Carbon: 0.5 # (Optional) Correction factor to apply to the Carbon yield from the SNII channel.
......
......@@ -344,6 +344,9 @@ INLINE static void evolve_SNIa(const float log10_min_mass,
sp->feedback_data.to_distribute.Fe_mass_from_SNIa +=
num_SNIa * props->yield_SNIa_IMF_resampled[chemistry_element_Fe] *
props->solar_mass_to_mass;
/* Compute the energy to be injected */
sp->feedback_data.to_distribute.d_energy += num_SNIa * props->E_SNIa;
}
/**
......@@ -717,7 +720,7 @@ void compute_stellar_evolution(const struct feedback_props* feedback_props,
/* Integration interval is zero - this can happen if minimum and maximum
* dying masses are above imf_max_mass_Msun. Return without doing any
* feedback. */
* enrichment. */
if (min_dying_mass_Msun == max_dying_mass_Msun) return;
/* Life is better in log */
......@@ -739,20 +742,27 @@ void compute_stellar_evolution(const struct feedback_props* feedback_props,
stellar_yields, feedback_props, sp);
}
#ifdef SWIFT_DEBUG_CHECKS
if (sp->feedback_data.to_distribute.mass != 0.f)
error("Injected mass will be lost");
#endif
/* 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 energy change due to kinetic energy of ejectas */
sp->feedback_data.to_distribute.d_energy +=
sp->feedback_data.to_distribute.mass *
feedback_props->AGB_ejecta_specific_kinetic_energy;
/* Compute energy change due to kinetic energy of the star */
sp->feedback_data.to_distribute.d_energy +=
sp->feedback_data.to_distribute.mass * 0.5f *
(sp->v[0] * sp->v[0] + sp->v[1] * sp->v[1] + sp->v[2] * sp->v[2]) *
cosmo->a2_inv;
}
/**
......@@ -867,6 +877,27 @@ void feedback_props_init(struct feedback_props* fp,
parser_get_param_float(params, "EAGLEFeedback:SNIa_timescale_Gyr");
fp->SNIa_timescale_Gyr_inv = 1.f / fp->SNIa_timescale_Gyr;
/* Energy released by supernova type Ia */
fp->E_SNIa_cgs =
parser_get_param_double(params, "EAGLEFeedback:SNIa_Energy_erg");
fp->E_SNIa =
fp->E_SNIa_cgs / units_cgs_conversion_factor(us, UNIT_CONV_ENERGY);
/* Properties of the SNIa enrichment model -------------------------------- */
/* Read AGB ejecta velocity */
const float ejecta_velocity_km_p_s = parser_get_param_float(
params, "EAGLEFeedback:AGB_ejecta_velocity_km_p_s");
/* Convert to internal units */
const float ejecta_velocity_cgs = ejecta_velocity_km_p_s * 1e5;
const float ejecta_velocity =
ejecta_velocity_cgs / units_cgs_conversion_factor(us, UNIT_CONV_SPEED);
/* Convert to specific thermal energy */
fp->AGB_ejecta_specific_kinetic_energy =
0.5f * ejecta_velocity * ejecta_velocity;
/* Gather common conversion factors --------------------------------------- */
/* Calculate internal mass to solar mass conversion factor */
......@@ -889,14 +920,6 @@ void feedback_props_init(struct feedback_props* fp,
fp->rho_to_n_cgs =
(X_H / m_p) * units_cgs_conversion_factor(us, UNIT_CONV_NUMBER_DENSITY);
// 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;
/* Initialise the IMF ------------------------------------------------- */
init_imf(fp);
......
......@@ -216,15 +216,14 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx,
pj->chemistry_data.metal_mass_fraction_from_AGB =
new_metal_mass_from_AGB / new_mass;
/* /\* Update momentum *\/ */
/* for (int i = 0; i < 3; i++) { */
/* pj->v[i] += si->feedback_data.to_distribute.mass * Omega_frac * */
/* (si->v[i] - pj->v[i]); */
/* } */
/* Update momentum */
for (int i = 0; i < 3; i++) {
pj->v[i] += si->feedback_data.to_distribute.mass * Omega_frac *
(si->v[i] - pj->v[i]);
}
/* Energy feedback */
// d_energy += si->feedback_data.to_distribute.d_energy *
// Omega_frac;
// d_energy += ;
/* if (feedback_props->continuous_heating) { */
/* // We're doing ONLY continuous heating */
......@@ -233,16 +232,18 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx,
/* 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); */
/* Add contribution from thermal and kinetic energy of ejected material
(and continuous SNIa feedback) */
/* const double u_init = hydro_get_physical_internal_energy(pj, xp, cosmo); */
/* const double delta_energy = si->feedback_data.to_distribute.d_energy *
* Omega_frac; */
/* const double delta_u = delta_energy / new_mass; */
/* const double u_new = u_init + delta_u; */
/* hydro_set_physical_internal_energy(pj, xp, cosmo, u_new); */
/* hydro_set_drifted_physical_internal_energy(pj, cosmo, u_new); */
/* 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) {
......@@ -254,8 +255,9 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx,
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;
const double u_init = hydro_get_physical_internal_energy(pj, xp, cosmo);
const float delta_u = si->feedback_data.to_distribute.SNII_delta_u;
const double u_new = u_init + delta_u;
/* Inject energy into the particle */
hydro_set_physical_internal_energy(pj, xp, cosmo, u_new);
......
......@@ -87,11 +87,6 @@ struct feedback_props {
/*! Are we doing SNIa enrichment? */
int with_SNIa_enrichment;
/* ------------ Energy injection ----------------- */
/* Kinetic energy of SN ejecta per unit mass (check name with Richard)*/
float ejecta_specific_thermal_energy;
/* ------------ Yield tables ----------------- */
/* Yield tables for AGB and SNII */
......@@ -138,6 +133,17 @@ struct feedback_props {
/*! Log 10 of the maximal mass used for SNIa feedback (in solar masses) */
float log10_SNIa_max_mass_msun;
/*! Energy released by one supernova type II in cgs units */
double E_SNIa_cgs;
/*! Energy released by one supernova type II in internal units */
float E_SNIa;
/* ------------- AGB parameters ---------------- */
/*! Specific kinetic energy injected from AGB ejectas (in internal units). */
float AGB_ejecta_specific_kinetic_energy;
/* ------------- Conversion factors --------------- */
/*! Conversion factor from internal mass unit to solar mass */
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment