diff --git a/src/feedback/EAGLE/feedback.c b/src/feedback/EAGLE/feedback.c index 7a86785b8371b1e3bad9f65e6c0a8ecc1188ea49..8d81eae0ac8932fcae39c8456710992048f73a80 100644 --- a/src/feedback/EAGLE/feedback.c +++ b/src/feedback/EAGLE/feedback.c @@ -76,9 +76,9 @@ double eagle_feedback_number_of_SNIa(const struct spart* sp, const double t0, * eq. 3 of Schaye 2015 paper. */ const double tau = props->SNIa_timescale_Gyr_inv; const double nu = props->SNIa_efficiency; - const double num_SNe_per_Msun = nu * (exp(-t0 * tau) - exp(-t1 * tau)); + const double num_SNIa_per_Msun = nu * (exp(-t0 * tau) - exp(-t1 * tau)); - return num_SNe_per_Msun * sp->mass_init * props->mass_to_solar_mass; + return num_SNIa_per_Msun * sp->mass_init * props->mass_to_solar_mass; } /** @@ -323,6 +323,11 @@ INLINE static void evolve_SNIa(const float log10_min_mass, star_age_Gyr = lifetime_Gyr; } +#ifdef SWIFT_DEBUG_CHECKS + if (dt_Gyr < 0.) error("Negative time-step length!"); + if (star_age_Gyr < 0.) error("Negative age!"); +#endif + /* Compute the number of SNIa */ const float num_SNIa = eagle_feedback_number_of_SNIa( sp, star_age_Gyr, star_age_Gyr + dt_Gyr, props); diff --git a/src/feedback/EAGLE/feedback_iact.h b/src/feedback/EAGLE/feedback_iact.h index 8d55820520efcb81ea664870fb630e6fefbd3b87..06fbd98a3cf5768a670506d689afffe9d95390a2 100644 --- a/src/feedback/EAGLE/feedback_iact.h +++ b/src/feedback/EAGLE/feedback_iact.h @@ -228,10 +228,9 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx, 0.5 * cosmo->a2_inv * current_mass * current_v2; /* Update velocity following injection of momentum */ - const double delta_m = si->feedback_data.to_distribute.mass * Omega_frac; - xpj->v_full[0] += delta_m * si->v[0] * new_mass_inv; - xpj->v_full[1] += delta_m * si->v[1] * new_mass_inv; - xpj->v_full[2] += delta_m * si->v[2] * new_mass_inv; + 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 new kinetic energy */ const double new_v2 = xpj->v_full[0] * xpj->v_full[0] + @@ -239,11 +238,13 @@ 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; - /* Total energy of that particle */ - const double new_total_energy = - current_kinetic_energy_gas + + /* Injection energy */ + 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; + /* 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; @@ -253,7 +254,8 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx, hydro_get_physical_internal_energy(pj, xpj, cosmo); #ifdef SWIFT_DEBUG_CHECKS - if (delta_u_enrich < 0.) error("Removing energy from the system."); + if (delta_u_enrich < -1e-3 * u_init_enrich) + error("Removing energy from the system."); #endif /* Do the energy injection.