diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml
index c3b8690564c6490cb1fae9ee09c73edfc9e70e27..4773205de440ca70b2f9ffcf334044009f08c501 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml
@@ -1,6 +1,6 @@
# 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
diff --git a/src/cell.c b/src/cell.c
index b64d802ab5e8e20800f54a940f1ea4f1265d9187..c9abe72cd62751e3c5348e4493ca4f922328ca12 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -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);
}
}
diff --git a/src/feedback/EAGLE/feedback.c b/src/feedback/EAGLE/feedback.c
index afdc23f80fb54878907436c0a677d2493b6fec32..3d4ff89791db2a471d2c44a5d6991908d9b6e23f 100644
--- a/src/feedback/EAGLE/feedback.c
+++ b/src/feedback/EAGLE/feedback.c
@@ -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");
}
diff --git a/src/feedback/EAGLE/feedback.h b/src/feedback/EAGLE/feedback.h
index 6e5aab1182b60c2e25fafb19b24de64daaccb1cf..ed1183a53c2196264bdae8c7a976ae001bfcf04e 100644
--- a/src/feedback/EAGLE/feedback.h
+++ b/src/feedback/EAGLE/feedback.h
@@ -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;
}
/**
diff --git a/src/feedback/EAGLE/feedback_iact.h b/src/feedback/EAGLE/feedback_iact.h
index 28b37e48a1b6b63af87317bffca4ae8e7242462d..08e82c06cefb7be653b9ef8536b3ee7a694cbd3c 100644
--- a/src/feedback/EAGLE/feedback_iact.h
+++ b/src/feedback/EAGLE/feedback_iact.h
@@ -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 */
diff --git a/src/feedback/EAGLE/feedback_properties.h b/src/feedback/EAGLE/feedback_properties.h
index fdc60739ef6936bc07692f143a6f38aeb28823c4..960abc5c5cdea90ccb3d16da419ffa3fd539fc09 100644
--- a/src/feedback/EAGLE/feedback_properties.h
+++ b/src/feedback/EAGLE/feedback_properties.h
@@ -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,
diff --git a/src/feedback/EAGLE/feedback_struct.h b/src/feedback/EAGLE/feedback_struct.h
index a1e2c5ba151501de858651b81030d91eb097f995..cd9c923224bef3c2ea75b86992759502e56c068b 100644
--- a/src/feedback/EAGLE/feedback_struct.h
+++ b/src/feedback/EAGLE/feedback_struct.h
@@ -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;
diff --git a/src/feedback/EAGLE/imf.h b/src/feedback/EAGLE/imf.h
index 84929b78a8ce306783c7b7e9400c2fa0e84308b2..c48c86d2ed4affcd70e906f77c94e14e4fe9a02b 100644
--- a/src/feedback/EAGLE/imf.h
+++ b/src/feedback/EAGLE/imf.h
@@ -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;
diff --git a/src/runner.c b/src/runner.c
index 0273c7d00f5ada379c16e8d37dd62535a3433318..4b393914f51171657999516d481abc5c68171158 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -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);