diff --git a/examples/SubgridTests/StellarEvolution/stellar_evolution.yml b/examples/SubgridTests/StellarEvolution/stellar_evolution.yml
index 559e415f706ebfbd0aa20db7b85dc79b61c48151..9baa1b57be2b14564c50aec2796d15fb3364aa17 100644
--- a/examples/SubgridTests/StellarEvolution/stellar_evolution.yml
+++ b/examples/SubgridTests/StellarEvolution/stellar_evolution.yml
@@ -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.
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 26a85f4bc2b40a9bfe21b3c7eaa5b25520e1e1be..d077a216890b04ce05e73110b34802f09a032dd5 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -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.
diff --git a/src/feedback/EAGLE/feedback.c b/src/feedback/EAGLE/feedback.c
index 612fc8c798bf9c797a5e2b82694790c8f431ceb2..d0d05d35fa3840bf0f82a5da0ad334f90c1694b9 100644
--- a/src/feedback/EAGLE/feedback.c
+++ b/src/feedback/EAGLE/feedback.c
@@ -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);
diff --git a/src/feedback/EAGLE/feedback_iact.h b/src/feedback/EAGLE/feedback_iact.h
index 151308696ede4d67bf1e26f7cb5a7ca54cf7e9a0..c32524dceaea5a544df32b70464d95040e5e29a4 100644
--- a/src/feedback/EAGLE/feedback_iact.h
+++ b/src/feedback/EAGLE/feedback_iact.h
@@ -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);
diff --git a/src/feedback/EAGLE/feedback_properties.h b/src/feedback/EAGLE/feedback_properties.h
index 3ea28b82d9689695236b239b7c5a3f7c7ffa2e9c..69aee1124c88ccb21a37d22dea5b720a3af8d0aa 100644
--- a/src/feedback/EAGLE/feedback_properties.h
+++ b/src/feedback/EAGLE/feedback_properties.h
@@ -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 */