From 4e6987ec70b60aced3d2820c176fc63a1ba4748f Mon Sep 17 00:00:00 2001
From: Loic Hausammann <loic.hausammann@protonmail.ch>
Date: Tue, 4 Aug 2020 09:17:52 +0100
Subject: [PATCH] Update GEAR feedback
---
src/feedback/GEAR/feedback.c | 9 ++++++++-
src/feedback/GEAR/stellar_evolution.c | 8 +++++---
2 files changed, 13 insertions(+), 4 deletions(-)
diff --git a/src/feedback/GEAR/feedback.c b/src/feedback/GEAR/feedback.c
index f4f8f30961..e32c72dd49 100644
--- a/src/feedback/GEAR/feedback.c
+++ b/src/feedback/GEAR/feedback.c
@@ -225,7 +225,13 @@ void feedback_evolve_spart(struct spart* restrict sp,
#ifdef SWIFT_DEBUG_CHECKS
if (sp->birth_time == -1.) error("Evolving a star particle that should not!");
+
+ if (star_age_beg_step < -1e-6) {
+ error("Negative age for a star");
+ }
#endif
+ const double star_age_beg_step_safe =
+ star_age_beg_step < 0 ? 0 : star_age_beg_step;
/* Reset the feedback */
feedback_reset_feedback(sp, feedback_props);
@@ -238,7 +244,8 @@ void feedback_evolve_spart(struct spart* restrict sp,
/* Compute the stellar evolution */
stellar_evolution_evolve_spart(sp, &feedback_props->stellar_model, cosmo, us,
- phys_const, ti_begin, star_age_beg_step, dt);
+ phys_const, ti_begin, star_age_beg_step_safe,
+ dt);
/* Transform the number of SN to the energy */
sp->feedback_data.energy_ejected =
diff --git a/src/feedback/GEAR/stellar_evolution.c b/src/feedback/GEAR/stellar_evolution.c
index cae77ee33e..2dcbc0b5ee 100644
--- a/src/feedback/GEAR/stellar_evolution.c
+++ b/src/feedback/GEAR/stellar_evolution.c
@@ -187,10 +187,12 @@ void stellar_evolution_compute_discrete_feedback_properties(
const float log_m_end_step, const float m_beg_step, const float m_end_step,
const float m_init, const int number_snia, const int number_snii) {
+ /* Limit the mass within the imf limits */
+ const float m_beg_lim = min(m_beg_step, sm->imf.mass_max);
+ const float m_end_lim = max(m_end_step, sm->imf.mass_min);
+
/* Compute the average mass */
- const float m_avg =
- initial_mass_function_get_integral_imf(&sm->imf, m_end_step, m_beg_step) /
- initial_mass_function_get_integral_xi(&sm->imf, m_end_step, m_beg_step);
+ const float m_avg = 0.5 * (m_beg_lim + m_end_lim);
const float log_m_avg = log10(m_avg);
/* Compute the mass ejected */
--
GitLab