Better debugging checks in the EAGLE feedback.

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);

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.