From cf8a9bb029dc544f0edb1d106504c264cc33569c Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Fri, 13 Aug 2021 09:39:21 +0100
Subject: [PATCH] add artificial accretion rate on a per-stellar type basis

fix so that the mass gain timestep multipliers work properly
---
 src/binary_c_structures.h                     |  1 +
 .../artificial_mass_accretion_rates.c         |  9 +++
 src/evolution/evolution_difflog.c             |  2 +-
 src/logging/log_flags.def                     |  1 +
 src/setup/cmd_line_args_list.h                |  9 +++
 src/stellar_structure/stellar_structure_AGB.c |  9 +++
 src/stellar_structure/stellar_structure_WD.c  |  3 +
 src/timestep/timestep_RLOF.c                  |  1 +
 src/timestep/timestep_limits.c                | 67 +++++++++++--------
 .../timestep_set_default_multipliers.c        |  2 +-
 tbse                                          | 12 +++-
 11 files changed, 84 insertions(+), 32 deletions(-)

diff --git a/src/binary_c_structures.h b/src/binary_c_structures.h
index 1739c0aed..d3739984b 100644
--- a/src/binary_c_structures.h
+++ b/src/binary_c_structures.h
@@ -744,6 +744,7 @@ struct preferences_t {
     double artificial_orbital_angular_momentum_accretion_rate;
     double artificial_accretion_start_time;
     double artificial_accretion_end_time;
+    double artificial_mass_accretion_rate_by_stellar_type[NUMBER_OF_STELLAR_TYPES];
     double minimum_donor_menv_for_comenv;
     double magnetic_braking_factor;
     double magnetic_braking_gamma;
diff --git a/src/binary_star_functions/artificial_mass_accretion_rates.c b/src/binary_star_functions/artificial_mass_accretion_rates.c
index 64dee80da..61b87c40a 100644
--- a/src/binary_star_functions/artificial_mass_accretion_rates.c
+++ b/src/binary_star_functions/artificial_mass_accretion_rates.c
@@ -12,8 +12,17 @@ void artificial_mass_accretion_rates(struct stardata_t * const stardata)
     {
         Foreach_star(star)
         {
+            /*
+             * Accrete at a fixed rate independent of stellar type
+             */
             star->derivative[DERIVATIVE_STELLAR_MASS_ARTIFICIAL] =
                 stardata->preferences->artificial_mass_accretion_rate[star->starnum];
+
+            /*
+             * Accrete at a fixed rate depending on stellar type
+             */
+            star->derivative[DERIVATIVE_STELLAR_MASS_ARTIFICIAL] =
+                stardata->preferences->artificial_mass_accretion_rate_by_stellar_type[star->stellar_type];
         }
     }
 }
diff --git a/src/evolution/evolution_difflog.c b/src/evolution/evolution_difflog.c
index c584dc6a2..801204cde 100644
--- a/src/evolution/evolution_difflog.c
+++ b/src/evolution/evolution_difflog.c
@@ -775,7 +775,7 @@ static Boolean cf_diffstats(struct stardata_t * RESTRICT const stardata,
         }
 
         /* core collapse supernovae */
-        if(stardata->model.supernova)//now->SN_type[k] != SN_NONE)
+        if(stardata->model.supernova)
         {
             flags[SN_LABEL] = TRUE;
         }
diff --git a/src/logging/log_flags.def b/src/logging/log_flags.def
index e54491e8b..e3069df7e 100644
--- a/src/logging/log_flags.def
+++ b/src/logging/log_flags.def
@@ -9,6 +9,7 @@
         X( PROBABILITY )                        \
         X( PERIASTRON_COLLISION )               \
         X( SN_AND_KICK )                        \
+        X( MCHAND )                             \
         X( OFF_MS )                             \
         X( STELLAR_TYPE_CHANGE )                \
         X( COMENV_POST )                        \
diff --git a/src/setup/cmd_line_args_list.h b/src/setup/cmd_line_args_list.h
index 1a7aa9ae0..e7932170c 100644
--- a/src/setup/cmd_line_args_list.h
+++ b/src/setup/cmd_line_args_list.h
@@ -1301,6 +1301,15 @@ BATCHMODE_ARGS
         Var(stardata->preferences->artificial_mass_accretion_rate),
         1.0
         CMD_LINE_ARG_T_OPTIONAL_ARGS },
+{
+    ARG_SECTION_STARS,
+        "artificial_mass_accretion_rate_by_stellar_type%d",
+        "Constant mass accretion rate for stellar type <n>.",
+        ARG_FLOAT_SCANF ,
+        WTTS_USE_DEFAULT,
+        Var(stardata->preferences->artificial_mass_accretion_rate_by_stellar_type),
+        1.0
+        CMD_LINE_ARG_T_OPTIONAL_ARGS },
 {
     ARG_SECTION_STARS,
         "artificial_angular_momentum_accretion_rate%d",
diff --git a/src/stellar_structure/stellar_structure_AGB.c b/src/stellar_structure/stellar_structure_AGB.c
index 3867874e4..27f0508c0 100644
--- a/src/stellar_structure/stellar_structure_AGB.c
+++ b/src/stellar_structure/stellar_structure_AGB.c
@@ -253,6 +253,15 @@ Stellar_type stellar_structure_AGB(struct star_t * const oldstar,
         (carbon_ignition != CARBON_IGNITION_NONE &&
          Boolean_(newstar->core_mass[CORE_CO] > mcmax));
 
+
+    if(supernova == TRUE &&
+       newstar->stellar_type == TPAGB)
+    {
+        Append_logstring(LOG_MCHAND,
+                         "0TPAGB core MCh reached");
+    }
+
+
     Dprint("supernova? %d\n",supernova);
 
     const Boolean envelope_lost =
diff --git a/src/stellar_structure/stellar_structure_WD.c b/src/stellar_structure/stellar_structure_WD.c
index 13c2b91b9..2ae4f8331 100644
--- a/src/stellar_structure/stellar_structure_WD.c
+++ b/src/stellar_structure/stellar_structure_WD.c
@@ -125,6 +125,9 @@ void stellar_structure_WD(struct star_t * const newstar,
                newstar->mass,
                MCh);
 
+        Append_logstring(LOG_MCHAND,
+                         "0WD MCh reached");
+
         struct star_t * news;
 
         if(newstar->stellar_type==ONeWD)
diff --git a/src/timestep/timestep_RLOF.c b/src/timestep/timestep_RLOF.c
index a433c60d2..6a58c1e5f 100644
--- a/src/timestep/timestep_RLOF.c
+++ b/src/timestep/timestep_RLOF.c
@@ -109,6 +109,7 @@ void timestep_RLOF(Timestep_prototype_args)
     if(Is_not_zero(stardata->model.RLOF_recommended_timestep))
     {
         Limit_timestep(*dt,
+                       stardata->preferences->timestep_multipliers[DT_LIMIT_RLOF]*
                        stardata->model.RLOF_recommended_timestep,
                        star,
                        DT_LIMIT_RLOF);
diff --git a/src/timestep/timestep_limits.c b/src/timestep/timestep_limits.c
index dee831040..368f40356 100644
--- a/src/timestep/timestep_limits.c
+++ b/src/timestep/timestep_limits.c
@@ -145,6 +145,45 @@ void timestep_limits(Timestep_prototype_args)
         }
     }
 
+
+    /*
+     * Mass gain: this usually requires that we are
+     * in a binary, but artificial accretion could
+     * also be set, so check it here.
+     */
+
+    {
+        const double accretion_rate = Mdot_gain(star);
+
+        if(Is_not_zero(accretion_rate))
+        {
+            if(Use_timestep(DT_LIMIT_MASS_GAIN))
+            {
+                const double dtlim = fsolver * 1e-6 *
+                    stardata->preferences->timestep_multipliers[DT_LIMIT_MASS_GAIN] *
+                    star->mass / accretion_rate;
+                Limit_timestep(*dt,
+                               dtlim,
+                               star,
+                               DT_LIMIT_MASS_GAIN);
+            }
+
+#ifdef SECOND_DERIVATIVES
+            if(Use_timestep(DT_LIMIT_MASS_GAIN2))
+            {
+                const double dtlim2 =
+                    stardata->preferences->timestep_multipliers[DT_LIMIT_MASS_GAIN2] *
+                    SECOND_DERIVATIVE_F * Second_derivative(DT_LIMIT_MASS_GAIN) * 1e-6;
+                Limit_timestep(*dt,
+                               dtlim2,
+                               star,
+                               DT_LIMIT_MASS_GAIN2);
+            }
+#endif//SECOND_DERIVATIVES
+        }
+    }
+
+
     if(System_is_binary)
     {
         /*
@@ -250,34 +289,6 @@ void timestep_limits(Timestep_prototype_args)
 
         timestep_RLOF(Timestep_call_args);
 
-        {
-            const double accretion_rate = Mdot_gain(star);
-
-            if(Is_not_zero(accretion_rate))
-            {
-                if(Use_timestep(DT_LIMIT_MASS_GAIN))
-                {
-                    const double dtlim = fsolver * 0.1 * star->mass / accretion_rate;
-                    Limit_timestep(*dt,
-                                   dtlim,
-                                   star,
-                                   DT_LIMIT_MASS_GAIN);
-                }
-
-#ifdef SECOND_DERIVATIVES
-                if(Use_timestep(DT_LIMIT_MASS_GAIN2))
-                {
-                    const double dtlim2 =  SECOND_DERIVATIVE_F * Second_derivative(DT_LIMIT_MASS_GAIN) * 1e-6;
-                    Limit_timestep(*dt,
-                                   dtlim2,
-                                   star,
-                                   DT_LIMIT_MASS_GAIN2);
-                }
-#endif//SECOND_DERIVATIVES
-            }
-        }
-
-
         if(stardata->preferences->individual_novae == TRUE &&
            WHITE_DWARF(star->stellar_type) &&
            (star->nova != NOVA_STATE_NONE ||
diff --git a/src/timestep/timestep_set_default_multipliers.c b/src/timestep/timestep_set_default_multipliers.c
index a00fc77c1..ad7b8f4fb 100644
--- a/src/timestep/timestep_set_default_multipliers.c
+++ b/src/timestep/timestep_set_default_multipliers.c
@@ -62,7 +62,7 @@ void timestep_set_default_multipliers(struct preferences_t * preferences)
     preferences->timestep_multipliers[DT_LIMIT_NOVAE] = 0.1; /* 39 */
     preferences->timestep_multipliers[DT_LIMIT_ARTIFICIAL_ACCRETION] = 1.0; /* 40 */
     preferences->timestep_multipliers[DT_LIMIT_SN] = 1.0; /* 41 */
-    preferences->timestep_multipliers[DT_LIMIT_MASS_GAIN] = 1.0; /* 42 */
+    preferences->timestep_multipliers[DT_LIMIT_MASS_GAIN] = 0.01; /* 42 */
     preferences->timestep_multipliers[DT_LIMIT_MASS_LOSS] = 1.0; /* 43 */
 #ifdef EXPONENTIAL_TIDES
     /* exponential tides are more stable, so can have a longer multiplier */
diff --git a/tbse b/tbse
index 219fdd075..e1c57ba23 100755
--- a/tbse
+++ b/tbse
@@ -382,10 +382,17 @@ ARTIFICIAL_ANGULAR_MOMENTUM_ACCRETION_RATE1=0.0
 ARTIFICIAL_ANGULAR_MOMENTUM_ACCRETION_RATE2=0.0
 ARTIFICIAL_ORBITAL_ANGULAR_MOMENTUM_ACCRETION_RATE=0.0
 # Start and end times for the accretion, ignored if set to
-# ARTIFICIAL_ACCRETION_IGNORE = -1
+# ARTIFICIAL_ACCRETION_IGNORE == -1
 ARTIFICIAL_ACCRETION_START_TIME=ARTIFICIAL_ACCRETION_IGNORE
 ARTIFICIAL_ACCRETION_END_TIME=ARTIFICIAL_ACCRETION_IGNORE
 
+# You can also accrete by stellar type, e.g.
+# to accrete on COWDs (type 11) at a rate 1e-10 Msun/yr
+# do this. This still respects the start and end time
+# set above.
+#ARTIFICIAL_MASS_ACCRETION_BY_STELLAR_TYPE="--artificial_mass_accretion_rate_by_stellar_type11 1e-10"
+ARTIFICIAL_MASS_ACCRETION_BY_STELLAR_TYPE=""
+
 # Mass loss at the helium flash
 # If >=0, this is an amount in Msun
 # If <0, this is an algorithm 
@@ -2519,6 +2526,7 @@ $WARMUP \
 --individual_novae $INDIVIDUAL_NOVAE \
 --artificial_mass_accretion_rate1 $ARTIFICIAL_MASS_ACCRETION_RATE1 \
 --artificial_mass_accretion_rate2 $ARTIFICIAL_MASS_ACCRETION_RATE2 \
+$ARTIFICIAL_MASS_ACCRETION_BY_STELLAR_TYPE \
 --artificial_angular_momentum_accretion_rate1 $ARTIFICIAL_ANGULAR_MOMENTUM_ACCRETION_RATE1 \
 --artificial_angular_momentum_accretion_rate2 $ARTIFICIAL_ANGULAR_MOMENTUM_ACCRETION_RATE2 \
 --artificial_accretion_start_time $ARTIFICIAL_ACCRETION_START_TIME \
@@ -2980,7 +2988,7 @@ else
     ############################################################
     
     eval "echo $ARGS_PREFIX $ARGS $ARGS_POSTFIX | xargs $EXECUTABLE"
-
+    
     STATUS=$?
     if [ "$VERBOSE" -eq "1" ];
     then
-- 
GitLab