diff --git a/src/feedback/EAGLE/feedback.c b/src/feedback/EAGLE/feedback.c
index 8d81eae0ac8932fcae39c8456710992048f73a80..76370156e4cda9c7abab8f538bd9d2517be1ebf8 100644
--- a/src/feedback/EAGLE/feedback.c
+++ b/src/feedback/EAGLE/feedback.c
@@ -690,6 +690,10 @@ void compute_stellar_evolution(const struct feedback_props* feedback_props,
TIMER_TIC;
+#ifdef SWIFT_DEBUG_CHECKS
+ if (age < 0.f) error("Negative age for a star.");
+#endif
+
/* Allocate temporary array for calculating imf weights */
float stellar_yields[eagle_feedback_N_imf_bins];
diff --git a/src/feedback/EAGLE/feedback.h b/src/feedback/EAGLE/feedback.h
index 5fa32fc2fa8053dbeaa1c1f96b5abb8526942b68..14f6821b065db5d05e408c1438af28d7492a5d94 100644
--- a/src/feedback/EAGLE/feedback.h
+++ b/src/feedback/EAGLE/feedback.h
@@ -44,6 +44,27 @@ __attribute__((always_inline)) INLINE static int feedback_do_feedback(
return (sp->birth_time != -1.);
}
+/**
+ * @brief Should this particle be doing any feedback-related operation?
+ *
+ * @param sp The #spart.
+ * @param time The current simulation time (Non-cosmological runs).
+ * @param cosmo The cosmological model (cosmological runs).
+ * @param with_cosmology Are we doing a cosmological run?
+ */
+__attribute__((always_inline)) INLINE static int feedback_is_active(
+ const struct spart* sp, const float time, const struct cosmology* cosmo,
+ const int with_cosmology) {
+
+ if (sp->birth_time == -1.) return 0;
+
+ if (with_cosmology) {
+ return ((float)cosmo->a) > sp->birth_scale_factor;
+ } else {
+ return time > sp->birth_time;
+ }
+}
+
/**
* @brief Prepares a s-particle for its feedback interactions
*
diff --git a/src/feedback/EAGLE/feedback_iact.h b/src/feedback/EAGLE/feedback_iact.h
index 06fbd98a3cf5768a670506d689afffe9d95390a2..63d5332f4d7eed815d7e32119d3ce6db67eb3af6 100644
--- a/src/feedback/EAGLE/feedback_iact.h
+++ b/src/feedback/EAGLE/feedback_iact.h
@@ -114,6 +114,11 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx,
Omega_frac = 0.f;
}
+#ifdef SWIFT_DEBUG_CHECKS
+ if (Omega_frac < 0. || Omega_frac > 1.)
+ error("Invalid fraction of material to dsitribute.");
+#endif
+
/* Update particle mass */
const double current_mass = hydro_get_mass(pj);
const double delta_mass = si->feedback_data.to_distribute.mass * Omega_frac;
@@ -227,10 +232,21 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx,
const double current_kinetic_energy_gas =
0.5 * cosmo->a2_inv * current_mass * current_v2;
- /* Update velocity following injection of momentum */
- xpj->v_full[0] += delta_mass * si->v[0] * new_mass_inv;
- xpj->v_full[1] += delta_mass * si->v[1] * new_mass_inv;
- xpj->v_full[2] += delta_mass * si->v[2] * new_mass_inv;
+ /* Compute the current thermal energy */
+ const double current_thermal_energy =
+ current_mass * hydro_get_physical_internal_energy(pj, xpj, cosmo);
+
+ /* Apply conservation of momentum */
+
+ /* Update velocity following change in gas mass */
+ xpj->v_full[0] *= current_mass * new_mass_inv;
+ xpj->v_full[1] *= current_mass * new_mass_inv;
+ xpj->v_full[2] *= current_mass * new_mass_inv;
+
+ /* Update velocity following addition of mass with different momentum */
+ xpj->v_full[0] += delta_mass * new_mass_inv * si->v[0];
+ xpj->v_full[1] += delta_mass * new_mass_inv * si->v[1];
+ xpj->v_full[2] += delta_mass * new_mass_inv * si->v[2];
/* Compute the new kinetic energy */
const double new_v2 = xpj->v_full[0] * xpj->v_full[0] +
@@ -238,32 +254,30 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx,
xpj->v_full[2] * xpj->v_full[2];
const double new_kinetic_energy_gas = 0.5 * cosmo->a2_inv * new_mass * new_v2;
- /* Injection energy */
+ /* Energy injected
+ * (thermal SNIa + kinetic energy of ejecta + kinetic energy of star) */
const double injected_energy =
si->feedback_data.to_distribute.energy * Omega_frac;
- /* Total energy of that particle */
- const double new_total_energy = current_kinetic_energy_gas + injected_energy;
+ /* Apply energy conservation to recover the new thermal energy of the gas */
+ const double new_thermal_energy = current_kinetic_energy_gas +
+ current_thermal_energy + injected_energy -
+ new_kinetic_energy_gas;
- /* Thermal energy of the particle */
- const double thermal_energy = new_total_energy - new_kinetic_energy_gas;
- const double delta_u_enrich = thermal_energy / new_mass;
-
- /* Energy feedback (ejecta energy + SNIa)*/
- const double u_init_enrich =
- hydro_get_physical_internal_energy(pj, xpj, cosmo);
+ /* Convert this to a specific thermal energy */
+ const double u_new_enrich = new_thermal_energy * new_mass_inv;
#ifdef SWIFT_DEBUG_CHECKS
- if (delta_u_enrich < -1e-3 * u_init_enrich)
- error("Removing energy from the system.");
+ if (new_thermal_energy < 0.99 * current_thermal_energy)
+ error("Enrichment is cooling the gas");
#endif
- /* Do the energy injection.
- * Note: We take a max() here just to be safe of rounding errors. */
- const double u_new_enrich = u_init_enrich + max(delta_u_enrich, 0.);
+ /* Do the energy injection. */
hydro_set_physical_internal_energy(pj, xpj, cosmo, u_new_enrich);
hydro_set_drifted_physical_internal_energy(pj, cosmo, u_new_enrich);
+ /* Finally, SNII stochastic feedback */
+
/* Get the SNII feedback properties */
const float prob = si->feedback_data.to_distribute.SNII_heating_probability;
diff --git a/src/feedback/none/feedback.h b/src/feedback/none/feedback.h
index de37015a9974df93d7dee81f7e386213f83b59bb..84a1cf36df9d95c2a672e59f8204f272b29a9b72 100644
--- a/src/feedback/none/feedback.h
+++ b/src/feedback/none/feedback.h
@@ -45,6 +45,21 @@ __attribute__((always_inline)) INLINE static int feedback_do_feedback(
return 0;
}
+/**
+ * @brief Should this particle be doing any feedback-related operation?
+ *
+ * @param sp The #spart.
+ * @param time The current simulation time (Non-cosmological runs).
+ * @param cosmo The cosmological model (cosmological runs).
+ * @param with_cosmology Are we doing a cosmological run?
+ */
+__attribute__((always_inline)) INLINE static int feedback_is_active(
+ const struct spart* sp, const float time, const struct cosmology* cosmo,
+ const int with_cosmology) {
+
+ return 0;
+}
+
/**
* @brief Prepares a star's feedback field before computing what
* needs to be distributed.
diff --git a/src/runner.c b/src/runner.c
index 183204faad02d0e255820f4c9fe470ffa8bd8f53..e56b7f05132a64606249490b8a5c6816cf291951 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -185,7 +185,8 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
if ((right = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
error("Can't allocate memory for right.");
for (int k = 0; k < c->stars.count; k++)
- if (spart_is_active(&sparts[k], e)) {
+ if (spart_is_active(&sparts[k], e) &&
+ feedback_is_active(&sparts[k], e->time, cosmo, with_cosmology)) {
sid[scount] = k;
h_0[scount] = sparts[k].h;
left[scount] = 0.f;
@@ -378,7 +379,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
star_age_end_of_step = cosmology_get_delta_time_from_scale_factors(
cosmo, sp->birth_scale_factor, (float)cosmo->a);
} else {
- star_age_end_of_step = e->time - sp->birth_time;
+ star_age_end_of_step = (float)e->time - sp->birth_time;
}
/* Has this star been around for a while ? */
@@ -1608,14 +1609,12 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
const float h_init = h_0[i];
const float h_old = p->h;
const float h_old_dim = pow_dimension(h_old);
- // const float h_old_inv_dim = pow_dimension(1.f / h_old);
const float h_old_dim_minus_one = pow_dimension_minus_one(h_old);
float h_new;
int has_no_neighbours = 0;
- if (p->density.wcount == 0.f) {
- // 1e-5 * kernel_root * h_old_inv_dim) { /* No neighbours case */
+ if (p->density.wcount == 0.f) { /* No neighbours case */
/* Flag that there were no neighbours */
has_no_neighbours = 1;
diff --git a/src/runner_doiact_stars.h b/src/runner_doiact_stars.h
index 50cdad07a5b25ab124c1b7cbfca83a81aed85efd..da279fcd01a7f2a3fe900863f7564451f6e88406 100644
--- a/src/runner_doiact_stars.h
+++ b/src/runner_doiact_stars.h
@@ -102,6 +102,7 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
TIMER_TIC;
const struct engine *e = r->e;
+ const int with_cosmology = e->policy & engine_policy_cosmology;
const integertime_t ti_current = e->ti_current;
const struct cosmology *cosmo = e->cosmology;
@@ -124,8 +125,13 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
/* Get a hold of the ith spart in ci. */
struct spart *restrict si = &sparts[sid];
+
+ /* Skip inactive particles */
if (!spart_is_active(si, e)) continue;
+ /* Skip inactive particles */
+ if (!feedback_is_active(si, e->time, cosmo, with_cosmology)) continue;
+
const float hi = si->h;
const float hig2 = hi * hi * kernel_gamma2;
const float six[3] = {(float)(si->x[0] - c->loc[0]),
@@ -191,6 +197,7 @@ void DO_NONSYM_PAIR1_STARS_NAIVE(struct runner *r, struct cell *restrict ci,
#endif
const struct engine *e = r->e;
+ const int with_cosmology = e->policy & engine_policy_cosmology;
const integertime_t ti_current = e->ti_current;
const struct cosmology *cosmo = e->cosmology;
@@ -222,8 +229,13 @@ void DO_NONSYM_PAIR1_STARS_NAIVE(struct runner *r, struct cell *restrict ci,
/* Get a hold of the ith spart in ci. */
struct spart *restrict si = &sparts_i[sid];
+
+ /* Skip inactive particles */
if (!spart_is_active(si, e)) continue;
+ /* Skip inactive particles */
+ if (!feedback_is_active(si, e->time, cosmo, with_cosmology)) continue;
+
const float hi = si->h;
const float hig2 = hi * hi * kernel_gamma2;
const float six[3] = {(float)(si->x[0] - (cj->loc[0] + shift[0])),
@@ -284,6 +296,7 @@ void DO_SYM_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
TIMER_TIC;
const struct engine *e = r->e;
+ const int with_cosmology = e->policy & engine_policy_cosmology;
const integertime_t ti_current = e->ti_current;
const struct cosmology *cosmo = e->cosmology;
@@ -350,6 +363,9 @@ void DO_SYM_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
/* Skip inactive particles */
if (!spart_is_active(spi, e)) continue;
+ /* Skip inactive particles */
+ if (!feedback_is_active(spi, e->time, cosmo, with_cosmology)) continue;
+
/* Compute distance from the other cell. */
const double px[3] = {spi->x[0], spi->x[1], spi->x[2]};
float dist = px[0] * runner_shift[sid][0] + px[1] * runner_shift[sid][1] +
@@ -475,6 +491,9 @@ void DO_SYM_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
/* Skip inactive particles */
if (!spart_is_active(spj, e)) continue;
+ /* Skip inactive particles */
+ if (!feedback_is_active(spj, e->time, cosmo, with_cosmology)) continue;
+
/* Compute distance from the other cell. */
const double px[3] = {spj->x[0], spj->x[1], spj->x[2]};
float dist = px[0] * runner_shift[sid][0] + px[1] * runner_shift[sid][1] +