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] +