From a7ce844d03c1e61d30d598369241f108494a232b Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Mon, 23 Nov 2020 14:20:05 +0100
Subject: [PATCH] Fix rounding issues in the new calculation of the enrichment
 dt in the time-step code

---
 .../CosmologicalStellarEvolution/run.sh       | 20 +++++++++----------
 .../stellar_evolution.yml                     |  4 +++-
 src/feedback/EAGLE/feedback.h                 | 18 +++++++++++------
 src/runner_time_integration.c                 | 10 +++++++---
 4 files changed, 32 insertions(+), 20 deletions(-)

diff --git a/examples/SubgridTests/CosmologicalStellarEvolution/run.sh b/examples/SubgridTests/CosmologicalStellarEvolution/run.sh
index 637a100567..ec267259e9 100755
--- a/examples/SubgridTests/CosmologicalStellarEvolution/run.sh
+++ b/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
diff --git a/examples/SubgridTests/CosmologicalStellarEvolution/stellar_evolution.yml b/examples/SubgridTests/CosmologicalStellarEvolution/stellar_evolution.yml
index e6b24b5194..032ac26df0 100644
--- a/examples/SubgridTests/CosmologicalStellarEvolution/stellar_evolution.yml
+++ b/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.
diff --git a/src/feedback/EAGLE/feedback.h b/src/feedback/EAGLE/feedback.h
index b8a905cac1..0a42bf94c7 100644
--- a/src/feedback/EAGLE/feedback.h
+++ b/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? */
diff --git a/src/runner_time_integration.c b/src/runner_time_integration.c
index b7ab9c6879..b7b014b731 100644
--- a/src/runner_time_integration.c
+++ b/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 */
-- 
GitLab