GEAR: fix the stellar feedback

src/feedback/GEAR/initial_mass_function.c

@@ -107,7 +107,7 @@ void initial_mass_function_integrate(const struct initial_mass_function *imf,
     }
 
     /* Integrate the data */
-    while (m < imf->mass_limits[i + 1] && j < count) {
+    while ((m < imf->mass_limits[i + 1] || i == imf->n_parts - 1) && j < count) {
 
       /* Compute the masses */
       const float log_m1 = log_mass_min + (j - 1) * step_size;

src/feedback/GEAR/interpolation.h

@@ -28,6 +28,9 @@ enum interpolate_boundary_condition {
 
   /* Zero (left boundary) and constant (right boundary) boundary conditions */
   boundary_condition_zero_const,
+
+  /* constant boundary conditions */
+  boundary_condition_const,
 };
 
 struct interpolation_1d {
@@ -96,6 +99,9 @@ __attribute__((always_inline)) static INLINE void interpolate_1d_init(
         case boundary_condition_zero_const:
           interp->data[i] = 0;
           break;
+        case boundary_condition_const:
+          interp->data[i] = data[0];
+          break;
         default:
           error("Interpolation type not implemented");
       }
@@ -109,6 +115,7 @@ __attribute__((always_inline)) static INLINE void interpolate_1d_init(
           interp->data[i] = 0;
           break;
         case boundary_condition_zero_const:
+        case boundary_condition_const:
           interp->data[i] = interp->data[i - 1];
           break;
         default:
@@ -149,6 +156,8 @@ __attribute__((always_inline)) static INLINE float interpolate_1d(
       case boundary_condition_zero:
       case boundary_condition_zero_const:
         return 0;
+      case boundary_condition_const:
+	return interp->data[0];
       default:
         error("Interpolation type not implemented");
     }
@@ -160,6 +169,7 @@ __attribute__((always_inline)) static INLINE float interpolate_1d(
       case boundary_condition_zero:
         return 0;
       case boundary_condition_zero_const:
+      case boundary_condition_const:
         return interp->data[interp->N - 1];
       default:
         error("Interpolation type not implemented");

src/feedback/GEAR/stellar_evolution.c

@@ -94,27 +94,28 @@ int stellar_evolution_compute_integer_number_supernovae(
  * (solMass)
  * @param m_end_step Mass of a star ending its life at the end of the step
  * (solMass)
- * @param number_snia Number of SNIa produced by the stellar particle.
- * @param number_snii Number of SNII produced by the stellar particle.
+ * @param number_snia_f (Floating) Number of SNIa produced by the stellar particle.
+ * @param number_snii_f (Floating) Number of SNII produced by the stellar particle.
  *
  */
 void stellar_evolution_compute_continuous_feedback_properties(
     struct spart* restrict sp, const struct stellar_model* sm,
     const struct phys_const* phys_const, const float log_m_beg_step,
     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) {
+    const float m_init, const float number_snia_f, const float number_snii_f) {
 
   /* Compute the mass ejected */
   /* SNIa */
   const float mass_snia =
-    supernovae_ia_get_ejected_mass_processed(&sm->snia) * number_snia;
+    supernovae_ia_get_ejected_mass_processed(&sm->snia) * number_snia_f;
 
   /* SNII */
   const float mass_frac_snii = supernovae_ii_get_ejected_mass_fraction_processed(
           &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 + mass_snia;
+  sp->feedback_data.mass_ejected = mass_frac_snii * sp->birth.mass
+    + mass_snia * phys_const->const_solar_mass;
 
   if (sp->mass <= sp->feedback_data.mass_ejected) {
     error("Stars cannot have negative mass. (%g <= %g). Initial mass = %g",
@@ -154,7 +155,7 @@ void stellar_evolution_compute_continuous_feedback_properties(
 
     /* Add the Supernovae Ia */
     sp->feedback_data.metal_mass_ejected[i] +=
-      snia_yields[i] * number_snia * phys_const->const_solar_mass;
+      snia_yields[i] * number_snia_f * phys_const->const_solar_mass;
 
   }
 }
@@ -185,18 +186,15 @@ void stellar_evolution_compute_discrete_feedback_properties(
 
   /* 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);
+    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 mass ejected */
   /* SNIa */
   const float mass_snia =
-      (number_snia == 0)
-          ? 0
-          : (supernovae_ia_get_ejected_mass_processed(&sm->snia) * number_snia);
+    (number_snia == 0) ? 0 :
+    (supernovae_ia_get_ejected_mass_processed(&sm->snia) * number_snia);
 
   /* SNII */
   const float mass_snii =
@@ -305,50 +303,36 @@ void stellar_evolution_evolve_spart(
   const float m_init = sp->sf_data.birth_mass / phys_const->const_solar_mass;
 
   /* Compute number of SNIa */
-  int number_snia = 0;
-  if (can_produce_snia) {
-    /* Compute rates */
-    const float number_snia_f =
-        supernovae_ia_get_number_per_unit_mass(&sm->snia, m_end_step, m_beg_step) * m_init;
-
-    /* Get the integer number of supernovae */
-    number_snia = stellar_evolution_compute_integer_number_supernovae(
-        sp, number_snia_f, ti_begin, random_number_stellar_feedback_1);
-
-    if (number_snia != 0) {
-      message("Exploding %i SNIa", number_snia);
-    }
-  }
+  const float number_snia_f = can_produce_snia ?
+    supernovae_ia_get_number_per_unit_mass(&sm->snia, m_end_step, m_beg_step) * m_init : 0;
 
   /* Compute number of SNII */
-  int number_snii = 0;
-  if (can_produce_snii) {
-    /* Compute rates */
-    const float number_snii_f =
-        supernovae_ii_get_number_per_unit_mass(&sm->snii, m_end_step, m_beg_step) * m_init;
-
-    /* Get the integer number of supernovae */
-    number_snii = stellar_evolution_compute_integer_number_supernovae(
-        sp, number_snii_f, ti_begin, random_number_stellar_feedback_2);
-    if (number_snii != 0) {
-      message("Exploding %i SNII", number_snii);
-    }
-  }
+  const float number_snii_f = can_produce_snii ?
+    supernovae_ii_get_number_per_unit_mass(&sm->snii, m_end_step, m_beg_step) * m_init : 0;
 
   /* Does this star produce a supernovae? */
-  if (number_snia == 0 && number_snii == 0) return;
+  if (number_snia_f == 0 && number_snii_f == 0) return;
 
-  sp->feedback_data.number_sn = number_snia + number_snii;
+  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 */
+    const int number_snia = stellar_evolution_compute_integer_number_supernovae(
+      sp, number_snia_f, ti_begin, random_number_stellar_feedback_1);
+
+    /* Get the integer number of supernovae */
+    const int number_snii = stellar_evolution_compute_integer_number_supernovae(
+      sp, number_snii_f, ti_begin, random_number_stellar_feedback_2);
+
     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 {
     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, number_snii);
+        m_end_step, m_init, number_snia_f, number_snii_f);
   }
 }
 

src/feedback/GEAR/stellar_evolution.h

@@ -39,7 +39,7 @@ void stellar_evolution_compute_continuous_feedback_properties(
     struct spart* restrict sp, const struct stellar_model* sm,
     const struct phys_const* phys_const, const float log_m_beg_step,
     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);
+    const float m_init, const float number_snia_f, const float number_snii_f);
 void stellar_evolution_compute_discrete_feedback_properties(
     struct spart* restrict sp, const struct stellar_model* sm,
     const struct phys_const* phys_const, const float log_m_beg_step,

src/feedback/GEAR/supernovae_ii.c

@@ -191,14 +191,13 @@ void supernovae_ii_read_yields_array(
   /* Read the dataset */
   io_read_array_dataset(group_id, hdf5_dataset_name, FLOAT, data, count);
 
-  /* Integrate the yields */
   initial_mass_function_integrate(&sm->imf, data, count, log_mass_min, step_size);
   // TODO: decrease count in order to keep the same distance between points
 
   /* Initialize the interpolation */
   interpolate_1d_init(interp, log10(snii->mass_min), log10(snii->mass_max),
                       interpolation_size, log_mass_min, step_size, count, data,
-		      boundary_condition_zero_const);
+		      boundary_condition_const);
 
 
   /* Cleanup the memory */