Fix rounding issues in the new calculation of the enrichment dt in the time-step code

examples/SubgridTests/CosmologicalStellarEvolution/run.sh

@@ -26,22 +26,22 @@ then
     ./getSolutions.sh
 fi
 
-../../swift  --feedback --stars --hydro --cosmology --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.08 2>&1 | tee output_0p08.log
+../../swift  --temperature --feedback --stars --hydro --cosmology --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.08 -P EAGLEChemistry:init_abundance_Hydrogen:0.71 -P EAGLEChemistry:init_abundance_Helium:0.21 2>&1 | tee output_0p08.log
 
-python plot_box_evolution.py
+python3 plot_box_evolution.py
 
-../../swift  --feedback --stars --hydro --cosmology --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.04 2>&1 | tee output_0p04.log
+../../swift  --temperature --feedback --stars --hydro --cosmology --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.04 -P EAGLEChemistry:init_abundance_Hydrogen:0.74 -P EAGLEChemistry:init_abundance_Helium:0.23 2>&1 | tee output_0p04.log
 
-python plot_box_evolution.py
+python3 plot_box_evolution.py
 
-../../swift  --feedback --stars --hydro --cosmology --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.01 2>&1 | tee output_0p01.log
+../../swift  --temperature --feedback --stars --hydro --cosmology --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.01 2>&1 | tee output_0p01.log
 
-python plot_box_evolution.py
+python3 plot_box_evolution.py
 
-../../swift  --feedback --stars --hydro --cosmology --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.001 2>&1 | tee output_0p001.log
+../../swift  --temperature --feedback --stars --hydro --cosmology --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.001 2>&1 | tee output_0p001.log
 
-python plot_box_evolution.py
+python3 plot_box_evolution.py
 
-../../swift  --feedback --stars --hydro --cosmology --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.0001 2>&1 | tee output_0p0001.log
+../../swift  --temperature --feedback --stars --hydro --cosmology --external-gravity --threads=4 stellar_evolution.yml -P EAGLEChemistry:init_abundance_metal:0.0001 2>&1 | tee output_0p0001.log
 
-python plot_box_evolution.py
+python3 plot_box_evolution.py

examples/SubgridTests/CosmologicalStellarEvolution/stellar_evolution.yml

@@ -96,6 +96,7 @@ EAGLEFeedback:
   IMF_max_mass_Msun:                  100.0             # Maximal stellar mass considered for the Chabrier IMF in solar masses.
   SNII_min_mass_Msun:                   6.0             # Minimal mass considered for SNII stars in solar masses.
   SNII_max_mass_Msun:                 100.0             # Maximal mass considered for SNII stars in solar masses.
+  SNII_feedback_model:                  Random          # Feedback modes: Random, Isotropic, MinimumDistance, MinimumDensity
   SNII_sampled_delay:                   0               # Sample the SNII lifetimes to do feedback.
   SNII_wind_delay_Gyr:                  0.03            # Time in Gyr between a star's birth and the SNII thermal feedback event when not sampling.
   SNII_delta_T_K:                       3.16228e7       # Change in temperature to apply to the gas particle in a SNII thermal feedback event in Kelvin.
@@ -106,7 +107,8 @@ EAGLEFeedback:
   SNII_energy_fraction_n_0_H_p_cm3:     0.67            # Pivot point for the birth density dependance of the SNII energy fraction in cm^-3.
   SNII_energy_fraction_n_Z:             0.8686          # Power-law for the metallicity dependance of the SNII energy fraction.
   SNII_energy_fraction_n_n:             0.8686          # Power-law for the birth density dependance of the SNII energy fraction.
-  SNII_energy_fraction_use_birth_props: 0               # Are we using the density and metallicity at birth to compute f_E or at feedback time?
+  SNII_energy_fraction_use_birth_density: 0             # Are we using the density and metallicity at birth to compute f_E or at feedback time?
+  SNII_energy_fraction_use_birth_metallicity: 0         # Are we using the density and metallicity at birth to compute f_E or at feedback time?
   SNIa_DTD:                             Exponential     # Use the EAGLE-Ref SNIa DTD.
   SNIa_DTD_delay_Gyr:                   0.04            # Age of the most massive SNIa in Gyr.
   SNIa_DTD_exp_timescale_Gyr:           2.0             # Time-scale of the exponential decay of the SNIa rates in Gyr.

src/feedback/EAGLE/feedback.h

@@ -99,10 +99,13 @@ INLINE static double feedback_get_enrichment_timestep(
     const struct cosmology* cosmo, const double time, const double dt_star) {
 
   if (with_cosmology) {
-    return cosmology_get_delta_time_from_scale_factors(
-        cosmo, (double)sp->last_enrichment_time, cosmo->a);
+    if (cosmo->a > (double)sp->last_enrichment_time)
+      return cosmology_get_delta_time_from_scale_factors(
+          cosmo, (double)sp->last_enrichment_time, cosmo->a);
+    else
+      return 0.;
   } else {
-    return time - sp->last_enrichment_time;
+    return max(time - sp->last_enrichment_time, 0.);
   }
 }
 
@@ -264,10 +267,13 @@ __attribute__((always_inline)) INLINE static void feedback_will_do_feedback(
   /* Calculate age of the star at current time */
   double age_of_star;
   if (with_cosmology) {
-    age_of_star = cosmology_get_delta_time_from_scale_factors(
-        cosmo, (double)sp->birth_scale_factor, cosmo->a);
+    if (cosmo->a > (double)sp->birth_scale_factor)
+      age_of_star = cosmology_get_delta_time_from_scale_factors(
+          cosmo, (double)sp->birth_scale_factor, cosmo->a);
+    else
+      age_of_star = 0.;
   } else {
-    age_of_star = time - (double)sp->birth_time;
+    age_of_star = max(time - (double)sp->birth_time, 0.);
   }
 
   /* Is the star still young? */

src/runner_time_integration.c

@@ -860,10 +860,14 @@ void runner_do_timestep(struct runner *r, struct cell *c, const int timer) {
           /* Calculate age of the star at current time */
           double star_age_end_of_step;
           if (with_cosmology) {
-            star_age_end_of_step = cosmology_get_delta_time_from_scale_factors(
-                cosmo, (double)sp->birth_scale_factor, cosmo->a);
+            if (cosmo->a > (double)sp->birth_scale_factor)
+              star_age_end_of_step =
+                  cosmology_get_delta_time_from_scale_factors(
+                      cosmo, (double)sp->birth_scale_factor, cosmo->a);
+            else
+              star_age_end_of_step = 0.;
           } else {
-            star_age_end_of_step = e->time - (double)sp->birth_time;
+            star_age_end_of_step = max(e->time - (double)sp->birth_time, 0.);
           }
 
           /* Get the length of the enrichment time-step */

[Matthieu Schaller]

@lhausammann I had to add this.

[Loic Hausammann]

Thanks for the info and the work!

[Matthieu Schaller]

We create a lot more unnecessary calculation in the non-GEAR case so I am actually not super happy with this change...

[Loic Hausammann]

Ok so let's rework it. I guess that you are not in an hurry to fix it, it can wait a week?

[Matthieu Schaller]

No hurry.