GEAR: discrete yields work correctly

src/feedback/GEAR/initial_mass_function.c

@@ -193,6 +193,13 @@ float initial_mass_function_get_coefficient(
 float initial_mass_function_get_integral_xi(
     const struct initial_mass_function *imf, float m1, float m2) {
 
+  /* Ensure the masses to be withing the limits */
+  m1 = min(m1, imf->mass_max);
+  m1 = max(m1, imf->mass_min);
+
+  m2 = min(m2, imf->mass_max);
+  m2 = max(m2, imf->mass_min);
+
   int k = -1;
   /* Find the correct part */
   for (int i = 0; i < imf->n_parts; i++) {
@@ -259,16 +266,14 @@ float initial_mass_function_get_imf(const struct initial_mass_function *imf,
  * @return The integral of the mass fraction.
  */
 float initial_mass_function_get_integral_imf(
-    const struct initial_mass_function *imf, const float m1, const float m2) {
+    const struct initial_mass_function *imf, float m1, float m2) {
 
-#ifdef SWIFT_DEBUG_CHECKS
-  if (m1 > imf->mass_max || m1 < imf->mass_min)
-    error("Mass 1 below or above limits expecting %g < %g < %g.", imf->mass_min,
-          m1, imf->mass_max);
-  if (m2 > imf->mass_max || m2 < imf->mass_min)
-    error("Mass 2 below or above limits expecting %g < %g < %g.", imf->mass_min,
-          m2, imf->mass_max);
-#endif
+  /* Ensure the masses to be withing the limits */
+  m1 = min(m1, imf->mass_max);
+  m1 = max(m1, imf->mass_min);
+
+  m2 = min(m2, imf->mass_max);
+  m2 = max(m2, imf->mass_min);
 
   for (int i = 0; i < imf->n_parts; i++) {
     if (m1 <= imf->mass_limits[i + 1]) {

src/feedback/GEAR/stellar_evolution.c

@@ -77,7 +77,7 @@ int stellar_evolution_compute_integer_number_supernovae(
   const float frac_sn = number_supernovae_f - number_supernovae_i;
 
   /* Get the integer number of SN */
-  return number_supernovae_i + ((rand_sn < frac_sn) ? 1 : 0);
+  return number_supernovae_i + (rand_sn < frac_sn);
 }
 
 /**
@@ -110,11 +110,11 @@ void stellar_evolution_compute_continuous_feedback_properties(
     supernovae_ia_get_ejected_mass_processed(&sm->snia) * number_snia_f;
 
   /* SNII */
-  const float mass_frac_snii = supernovae_ii_get_ejected_mass_fraction_processed(
+  const float mass_frac_snii = supernovae_ii_get_ejected_mass_fraction_processed_from_integral(
           &sm->snii, log_m_end_step, log_m_beg_step); 
 
   /* Sum the contributions from SNIa and SNII */
-  sp->feedback_data.mass_ejected = mass_frac_snii * sp->birth.mass
+  sp->feedback_data.mass_ejected = mass_frac_snii * sp->sf_data.birth_mass
     + mass_snia * phys_const->const_solar_mass;
 
   if (sp->mass <= sp->feedback_data.mass_ejected) {
@@ -132,11 +132,11 @@ void stellar_evolution_compute_continuous_feedback_properties(
 
   /* Compute the SNII yields */
   float snii_yields[GEAR_CHEMISTRY_ELEMENT_COUNT];
-  supernovae_ii_get_yields(&sm->snii, log_m_end_step, log_m_beg_step,
+  supernovae_ii_get_yields_from_integral(&sm->snii, log_m_end_step, log_m_beg_step,
                            snii_yields);
 
   /* Compute the mass fraction of non processed elements */
-  const float non_processed = supernovae_ii_get_ejected_mass_fraction(
+  const float non_processed = supernovae_ii_get_ejected_mass_fraction_from_integral(
       &sm->snii, log_m_end_step, log_m_beg_step);
 
   /* Set the yields */
@@ -174,6 +174,7 @@ void stellar_evolution_compute_continuous_feedback_properties(
  * (solMass)
  * @param m_end_step Mass of a star ending its life at the end of the step
  * (solMass)
+ * @param m_init Birth mass in solMass.
  * @param number_snia Number of SNIa produced by the stellar particle.
  * @param number_snii Number of SNII produced by the stellar particle.
  *
@@ -184,22 +185,20 @@ 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) {
 
-  /* Get the normalization to the average */
-  const float normalization =
-    number_snii == 0 ? 0. :
-    number_snii / (supernovae_ii_get_number_per_unit_mass(&sm->snii, m_end_step, m_beg_step) *
-		   m_init);
+  /* 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 log_m_avg = log10(m_avg);
 
   /* Compute the mass ejected */
   /* SNIa */
   const float mass_snia =
-    (number_snia == 0) ? 0 :
-    (supernovae_ia_get_ejected_mass_processed(&sm->snia) * number_snia);
+    supernovae_ia_get_ejected_mass_processed(&sm->snia) * number_snia;
 
   /* SNII */
   const float mass_snii =
-      normalization * supernovae_ii_get_ejected_mass_fraction_processed(
-                          &sm->snii, log_m_end_step, log_m_beg_step);
+    supernovae_ii_get_ejected_mass_fraction_processed_from_raw(
+      &sm->snii, log_m_avg) * m_avg * number_snii;
 
   sp->feedback_data.mass_ejected = mass_snia + mass_snii;
 
@@ -217,29 +216,33 @@ void stellar_evolution_compute_discrete_feedback_properties(
   /* Get the SNIa yields */
   const float* snia_yields = supernovae_ia_get_yields(&sm->snia);
 
-  /* Compute the SNII yields (without the normalization) */
+  /* Compute the SNII yields */
   float snii_yields[GEAR_CHEMISTRY_ELEMENT_COUNT];
-  supernovae_ii_get_yields(&sm->snii, log_m_end_step, log_m_beg_step,
-                           snii_yields);
+  supernovae_ii_get_yields_from_raw(&sm->snii, log_m_avg, snii_yields);
 
   /* Compute the mass fraction of non processed elements */
   const float non_processed =
-      normalization * supernovae_ii_get_ejected_mass_fraction(
-                          &sm->snii, log_m_end_step, log_m_beg_step);
+    supernovae_ii_get_ejected_mass_fraction_from_raw(
+        &sm->snii, log_m_avg);
 
   /* Set the yields */
   for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
     /* Compute the mass fraction of metals */
     sp->feedback_data.metal_mass_ejected[i] =
-        /* Supernovae Ia yields */
-        snia_yields[i] * number_snia +
-        /* Supernovae II yields */
-        normalization * snii_yields[i] +
-        /* Gas contained in stars initial metallicity */
-        chemistry_get_metal_mass_fraction_for_feedback(sp)[i] * non_processed;
+      /* Supernovae II yields */
+      snii_yields[i] +
+      /* Gas contained in stars initial metallicity */
+      chemistry_get_metal_mass_fraction_for_feedback(sp)[i] * non_processed;
 
     /* Convert it to total mass */
-    sp->feedback_data.metal_mass_ejected[i] *= sp->sf_data.birth_mass;
+    sp->feedback_data.metal_mass_ejected[i] *= m_avg * number_snii;
+
+    /* Supernovae Ia yields */
+    sp->feedback_data.metal_mass_ejected[i] +=
+      snia_yields[i] * number_snia;
+
+    /* Convert everything in code units */
+    sp->feedback_data.metal_mass_ejected[i] *= phys_const->const_solar_mass;
   }
 }
 
@@ -313,8 +316,6 @@ void stellar_evolution_evolve_spart(
   /* Does this star produce a supernovae? */
   if (number_snia_f == 0 && number_snii_f == 0) return;
 
-  sp->feedback_data.number_sn = number_snia_f + number_snii_f;
-
   /* Compute the properties of the feedback (e.g. yields) */
   if (sm->discrete_yields) {
     /* Get the integer number of supernovae */
@@ -325,11 +326,22 @@ void stellar_evolution_evolve_spart(
     const int number_snii = stellar_evolution_compute_integer_number_supernovae(
       sp, number_snii_f, ti_begin, random_number_stellar_feedback_2);
 
+    /* Do we have a supernovae? */
+    if (number_snia == 0 && number_snii == 0) return;
+
+    /* Save the number of supernovae */
+    sp->feedback_data.number_sn = number_snia + number_snii;
+
+    /* Compute the yields */
     stellar_evolution_compute_discrete_feedback_properties(
         sp, sm, phys_const, log_m_beg_step, log_m_end_step, m_beg_step,
         m_end_step, m_init, number_snia, number_snii);
 
   } else {
+    /* Save the number of supernovae */
+    sp->feedback_data.number_sn = number_snia_f + number_snii_f;
+
+    /* Compute the yields */
     stellar_evolution_compute_continuous_feedback_properties(
         sp, sm, phys_const, log_m_beg_step, log_m_end_step, m_beg_step,
         m_end_step, m_init, number_snia_f, number_snii_f);

src/feedback/GEAR/stellar_evolution_struct.h

@@ -113,15 +113,30 @@ struct supernovae_ia {
  */
 struct supernovae_ii {
 
-  /*! Integrated (over the IMF) mass fraction of metals ejected by a supernovae
-   */
-  struct interpolation_1d integrated_yields[GEAR_CHEMISTRY_ELEMENT_COUNT];
+  /*! Yields not integrated */
+  struct {
+    /*! Mass fraction of metals ejected by a supernovae. */
+    struct interpolation_1d yields[GEAR_CHEMISTRY_ELEMENT_COUNT];
+    
+    /*! Total mass fraction ejected. */
+    struct interpolation_1d ejected_mass_processed;
 
-  /*! Total mass fraction ejected (integrated over the IMF) */
-  struct interpolation_1d integrated_ejected_mass_processed;
+    /*! Mass fraction ejected and not processed (=> with the star metallicity). */
+    struct interpolation_1d ejected_mass;
+  } raw;
 
-  /*! Mass fraction ejected and not processed (=> with the star metallicity) */
-  struct interpolation_1d integrated_ejected_mass;
+  /*! Yields integrated */
+  struct {
+    /*! Integrated (over the IMF) mass fraction of metals ejected by a supernovae
+     */
+    struct interpolation_1d yields[GEAR_CHEMISTRY_ELEMENT_COUNT];
+    
+    /*! Total mass fraction ejected (integrated over the IMF) */
+    struct interpolation_1d ejected_mass_processed;
+
+    /*! Mass fraction ejected and not processed (=> with the star metallicity) */
+    struct interpolation_1d ejected_mass;
+  } integrated;
 
   /*! Minimal mass for a SNII */
   float mass_min;

src/feedback/GEAR/supernovae_ii.h

@@ -33,15 +33,22 @@ int supernovae_ii_can_explode(const struct supernovae_ii *snii, float m_low,
                               float m_high);
 float supernovae_ii_get_number_per_unit_mass(const struct supernovae_ii *snii, float m1,
 					     float m2);
-void supernovae_ii_get_yields(const struct supernovae_ii *snii, float log_m1,
-                              float log_m2, float *yields);
-float supernovae_ii_get_ejected_mass_fraction(const struct supernovae_ii *snii,
-                                              float log_m1, float log_m2);
+void supernovae_ii_get_yields_from_integral(const struct supernovae_ii *snii, float log_m1,
+					    float log_m2, float *yields);
+void supernovae_ii_get_yields_from_raw(const struct supernovae_ii *snii, float log_m,
+				       float *yields);
+float supernovae_ii_get_ejected_mass_fraction_from_integral(const struct supernovae_ii *snii,
+							    float log_m1, float log_m2);
+float supernovae_ii_get_ejected_mass_fraction_from_raw(const struct supernovae_ii *snii,
+						       float log_m);
 
-float supernovae_ii_get_ejected_mass_fraction_processed(
+float supernovae_ii_get_ejected_mass_fraction_processed_from_integral(
     const struct supernovae_ii *snii, float log_m1, float log_m2);
+float supernovae_ii_get_ejected_mass_fraction_processed_from_raw(
+    const struct supernovae_ii *snii, float log_m);
 void supernovae_ii_read_yields_array(
-    struct supernovae_ii *snii, struct interpolation_1d *interp,
+    struct supernovae_ii *snii, struct interpolation_1d *interp_raw,
+    struct interpolation_1d *interp_int,
     const struct phys_const *phys_const, const struct stellar_model *sm,
     hid_t group_id, const char *hdf5_dataset_name, hsize_t *previous_count,
     int interpolation_size);