Improve the logic in the stars ghost

src/feedback/EAGLE/feedback.c

@@ -651,38 +651,6 @@ INLINE static void evolve_AGB(const float log10_min_mass, float log10_max_mass,
   }
 }
 
-/**
- * @brief Prepares a star's feedback field before computing what
- * needs to be distributed.
- */
-INLINE static void feedback_init_to_distribute(struct spart* sp) {
-
-  /* Zero the amount of mass that is distributed */
-  sp->feedback_data.to_distribute.mass = 0.f;
-
-  /* Zero the metal enrichment quantities */
-  for (int i = 0; i < chemistry_element_count; i++) {
-    sp->feedback_data.to_distribute.metal_mass[i] = 0.f;
-  }
-  sp->feedback_data.to_distribute.total_metal_mass = 0.f;
-  sp->feedback_data.to_distribute.mass_from_AGB = 0.f;
-  sp->feedback_data.to_distribute.metal_mass_from_AGB = 0.f;
-  sp->feedback_data.to_distribute.mass_from_SNII = 0.f;
-  sp->feedback_data.to_distribute.metal_mass_from_SNII = 0.f;
-  sp->feedback_data.to_distribute.mass_from_SNIa = 0.f;
-  sp->feedback_data.to_distribute.metal_mass_from_SNIa = 0.f;
-  sp->feedback_data.to_distribute.Fe_mass_from_SNIa = 0.f;
-
-  /* Zero the energy to inject */
-  sp->feedback_data.to_distribute.d_energy = 0.f;
-
-  /* Zero the SNII feedback probability */
-  sp->feedback_data.to_distribute.SNII_heating_probability = 0.f;
-
-  /* Zero the SNII feedback energy */
-  sp->feedback_data.to_distribute.SNII_delta_u = 0.f;
-}
-
 /**
  * @brief calculates stellar mass in spart that died over the timestep, calls
  * functions to calculate feedback due to SNIa, SNII and AGB
@@ -721,10 +689,11 @@ void compute_stellar_evolution(const struct feedback_props* feedback_props,
    * particles. Do _NOT_ read from the to_collect substructure any more. */
 
   /* Zero all the output fields */
-  feedback_init_to_distribute(sp);
+  feedback_reset_feedback(sp, feedback_props);
 
   /* Update the weights used for distribution */
-  const float enrichment_weight = 1.f / enrichment_weight_inv;
+  const float enrichment_weight =
+      (enrichment_weight_inv != 0.f) ? 1.f / enrichment_weight_inv : 0.f;
   sp->feedback_data.to_distribute.enrichment_weight = enrichment_weight;
 
   /* Compute properties of the stochastic SNe feedback model. */

src/feedback/EAGLE/feedback.h

@@ -33,6 +33,17 @@ void compute_stellar_evolution(const struct feedback_props* feedback_props,
                                const struct unit_system* us, const float age,
                                const float dt);
 
+/**
+ * @brief Should we do feedback for this star?
+ *
+ * @param sp The star to consider.
+ */
+__attribute__((always_inline)) INLINE static int feedback_do_feedback(
+    const struct spart* sp) {
+
+  return (sp->birth_time != -1.);
+}
+
 /**
  * @brief Prepares a s-particle for its feedback interactions
  *
@@ -45,6 +56,42 @@ __attribute__((always_inline)) INLINE static void feedback_init_spart(
   sp->feedback_data.to_collect.ngb_mass = 0.f;
 }
 
+/**
+ * @brief Prepares a star's feedback field before computing what
+ * needs to be distributed.
+ */
+__attribute__((always_inline)) INLINE static void feedback_reset_feedback(
+    struct spart* sp, const struct feedback_props* feedback_props) {
+
+  /* Zero the distribution weights */
+  sp->feedback_data.to_distribute.enrichment_weight = 0.f;
+
+  /* Zero the amount of mass that is distributed */
+  sp->feedback_data.to_distribute.mass = 0.f;
+
+  /* Zero the metal enrichment quantities */
+  for (int i = 0; i < chemistry_element_count; i++) {
+    sp->feedback_data.to_distribute.metal_mass[i] = 0.f;
+  }
+  sp->feedback_data.to_distribute.total_metal_mass = 0.f;
+  sp->feedback_data.to_distribute.mass_from_AGB = 0.f;
+  sp->feedback_data.to_distribute.metal_mass_from_AGB = 0.f;
+  sp->feedback_data.to_distribute.mass_from_SNII = 0.f;
+  sp->feedback_data.to_distribute.metal_mass_from_SNII = 0.f;
+  sp->feedback_data.to_distribute.mass_from_SNIa = 0.f;
+  sp->feedback_data.to_distribute.metal_mass_from_SNIa = 0.f;
+  sp->feedback_data.to_distribute.Fe_mass_from_SNIa = 0.f;
+
+  /* Zero the energy to inject */
+  sp->feedback_data.to_distribute.d_energy = 0.f;
+
+  /* Zero the SNII feedback probability */
+  sp->feedback_data.to_distribute.SNII_heating_probability = 0.f;
+
+  /* Zero the SNII feedback energy */
+  sp->feedback_data.to_distribute.SNII_delta_u = 0.f;
+}
+
 /**
  * @brief Initialises the s-particles feedback props for the first time
  *
@@ -91,6 +138,10 @@ __attribute__((always_inline)) INLINE static void feedback_evolve_spart(
     const struct cosmology* cosmo, const struct unit_system* us,
     const double star_age_beg_step, const double dt) {
 
+#ifdef SWIFT_DEBUG_CHECKS
+  if (sp->birth_time == -1.) error("Evolving a star particle that shoul not!");
+#endif
+
   /* Compute amount of enrichment and feedback that needs to be done in this
    * step */
   compute_stellar_evolution(feedback_props, cosmo, sp, us, star_age_beg_step,

src/feedback/none/feedback.h

@@ -34,6 +34,24 @@
 __attribute__((always_inline)) INLINE static void feedback_init_spart(
     struct spart* sp) {}
 
+/**
+ * @brief Should we do feedback for this star?
+ *
+ * @param sp The star to consider.
+ */
+__attribute__((always_inline)) INLINE static int feedback_do_feedback(
+    const struct spart* sp) {
+
+  return 0;
+}
+
+/**
+ * @brief Prepares a star's feedback field before computing what
+ * needs to be distributed.
+ */
+__attribute__((always_inline)) INLINE static void feedback_reset_feedback(
+    struct spart* sp, const struct feedback_props* feedback_props) {}
+
 /**
  * @brief Initialises the s-particles feedback props for the first time
  *

src/runner.c

@@ -260,6 +260,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
               ((sp->h <= stars_h_min) && (f > 0.f))) {
 
             stars_reset_feedback(sp);
+            feedback_reset_feedback(sp, feedback_props);
 
             /* Ok, we are done with this particle */
             continue;
@@ -356,7 +357,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
         stars_reset_feedback(sp);
 
         /* Only do feedback if stars have a reasonable birth time */
-        if (sp->birth_time != -1.) {
+        if (feedback_do_feedback(sp)) {
 
           const integertime_t ti_step = get_integer_timestep(sp->time_bin);
           const integertime_t ti_begin =
@@ -380,9 +381,6 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
             star_age_end_of_step = e->time - sp->birth_time;
           }
 
-          /* Reset the feedback fields of the star particle */
-          feedback_prepare_spart(sp, feedback_props);
-
           /* Has this star been around for a while ? */
           if (star_age_end_of_step > 0.) {
 
@@ -393,7 +391,15 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
             /* Compute the stellar evolution  */
             feedback_evolve_spart(sp, feedback_props, cosmo, us,
                                   star_age_beg_of_step, dt);
+          } else {
+
+            /* Reset the feedback fields of the star particle */
+            feedback_reset_feedback(sp, feedback_props);
           }
+        } else {
+
+          /* Reset the feedback fields of the star particle */
+          feedback_reset_feedback(sp, feedback_props);
         }
       }
 
@@ -725,8 +731,6 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
             /* Did we get a star? (Or did we run out of spare ones?) */
             if (sp != NULL) {
 
-              message("Formed a star ID=%lld", sp->id);
-
               /* Copy the properties of the gas particle to the star particle */
               star_formation_copy_properties(p, xp, sp, e, sf_props, cosmo,
                                              with_cosmology);
@@ -1601,12 +1605,14 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
         const float h_init = h_0[i];
         const float h_old = p->h;
         const float h_old_dim = pow_dimension(h_old);
+        const float h_old_inv_dim = pow_dimension(1.f / h_old);
         const float h_old_dim_minus_one = pow_dimension_minus_one(h_old);
 
         float h_new;
         int has_no_neighbours = 0;
 
-        if (p->density.wcount == 0.f) { /* No neighbours case */
+        if (p->density.wcount <
+            1e-5 * kernel_root * h_old_inv_dim) { /* No neighbours case */
 
           /* Flag that there were no neighbours */
           has_no_neighbours = 1;