diff --git a/src/feedback/EAGLE/feedback_iact.h b/src/feedback/EAGLE/feedback_iact.h index 06fbd98a3cf5768a670506d689afffe9d95390a2..e6a3a9ec9340fd72eb707fab016c692256c234f2 100644 --- a/src/feedback/EAGLE/feedback_iact.h +++ b/src/feedback/EAGLE/feedback_iact.h @@ -227,10 +227,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 +249,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;