From bf840731f0949e6476b11ef4c8b878312011c5e7 Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Sat, 10 Jul 2021 09:26:30 +0100
Subject: [PATCH] add function supernova_should_split_evolution.c to check
whether evolution splitting should happen
add preference evolution_splitting_sn_eccentricity_threshold to set the eccentricity above which evolution splitting should happen
---
src/binary_c_function_macros.h | 2 +
src/binary_c_structures.h | 3 +
src/evolution/evolution_split.c | 43 +++++------
src/evolution/evolve_system_binary_c.c | 1 -
src/setup/cmd_line_args_list.h | 9 +++
src/setup/cmd_line_macros.h | 2 +
src/setup/set_default_preferences.c | 4 +-
src/supernovae/set_kick_velocity.c | 5 +-
.../supernova_should_split_evolution.c | 71 +++++++++++++++++++
src/supernovae/supernovae_prototypes.h | 2 +
tbse | 7 +-
11 files changed, 122 insertions(+), 27 deletions(-)
create mode 100644 src/supernovae/supernova_should_split_evolution.c
diff --git a/src/binary_c_function_macros.h b/src/binary_c_function_macros.h
index 11522bb72..7c35b432b 100644
--- a/src/binary_c_function_macros.h
+++ b/src/binary_c_function_macros.h
@@ -475,6 +475,8 @@
#define System_is_single _System_is_single(stardata)
#define System_is_binary _System_is_binary(stardata)
+#define System_was_single _System_is_single(stardata->previous_stardata)
+#define System_was_binary _System_is_binary(stardata->previous_stardata)
/*
* Define an observable binary as one in which
diff --git a/src/binary_c_structures.h b/src/binary_c_structures.h
index 5420c8318..8e2f3a58d 100644
--- a/src/binary_c_structures.h
+++ b/src/binary_c_structures.h
@@ -987,6 +987,9 @@ struct preferences_t {
double orbiting_objects_close_pc_threshold;
#endif // ORBITING_OBJECTS
+#ifdef EVOLUTION_SPLITTING
+ double evolution_splitting_sn_eccentricity_threshold;
+#endif
/*
* Strings (e.g. Filenames)
diff --git a/src/evolution/evolution_split.c b/src/evolution/evolution_split.c
index 2e30b946f..4682c1f8d 100644
--- a/src/evolution/evolution_split.c
+++ b/src/evolution/evolution_split.c
@@ -165,8 +165,9 @@ int evolution_split(struct stardata_t * RESTRICT const stardata,
stardata->star[1].stellar_type
);
Append_logstring(LOG_SPLIT_EVOLUTION,
- "Evolution restart @ t=%g, %d.%d",
+ "Evolution restart @ t=%g, prev e=%g, %d.%d",
stardata->model.time,
+ stardata->previous_stardata->common.orbit.eccentricity,
stardata->preferences->current_splitdepth-1,
p->count);
#ifdef LINUX
@@ -290,8 +291,9 @@ int evolution_split(struct stardata_t * RESTRICT const stardata,
stardata->model.doingsplit = TRUE;
Append_logstring(LOG_SPLIT_EVOLUTION,
- "Evolution split and restart at t=%g, %d.%d",
+ "Evolution split and restart at t=%g, prev e=%g, %d.%d",
stardata->model.time,
+ stardata->previous_stardata->common.orbit.eccentricity,
stardata->preferences->current_splitdepth-1,
p->count);
@@ -322,10 +324,6 @@ int evolution_split(struct stardata_t * RESTRICT const stardata,
static int splittest(struct stardata_t * RESTRICT const stardata)
{
- struct stardata_t * previous_stardata = stardata->previous_stardata;
- int split = 0;
- Star_number i;
-
/*
* Detect new split cases after stellar evolution has
* been done in this system. You can use stardata and
@@ -338,25 +336,28 @@ static int splittest(struct stardata_t * RESTRICT const stardata)
* Returns zero if there should be no split.
*/
- Starloop(i)
+ Foreach_star(star)
{
- if(!split)
+ /*
+ * Supernova check
+ */
+ SDEBUG("check SN %d : type %d\n",
+ stardata->model.supernova,
+ star->stellar_type);
+
+ const Boolean supernova =
+ stardata->model.supernova ||
+ ((star->stellar_type == BLACK_HOLE ||
+ star->stellar_type == NEUTRON_STAR) &&
+ (stardata->previous_stardata->star[star->starnum].stellar_type<=HeGB));
+
+ if(supernova == TRUE &&
+ supernova_should_split_evolution(stardata))
{
- // Supernova
- SDEBUG("check SN %d : type %d\n",
- stardata->model.supernova,
- stardata->star[i].stellar_type);
-
- if(stardata->model.supernova ||
- ((stardata->star[i].stellar_type == BLACK_HOLE ||
- stardata->star[i].stellar_type == NEUTRON_STAR) &&
- (previous_stardata->star[i].stellar_type<=HeGB)))
- {
- split = stardata->preferences->evolution_splitting_sn_n;
- }
+ return stardata->preferences->evolution_splitting_sn_n;
}
}
- return split;
+ return 0;
}
static void modulate_stardata(struct stardata_t * RESTRICT const stardata,
diff --git a/src/evolution/evolve_system_binary_c.c b/src/evolution/evolve_system_binary_c.c
index 45c10be27..836a42ea7 100644
--- a/src/evolution/evolve_system_binary_c.c
+++ b/src/evolution/evolve_system_binary_c.c
@@ -287,7 +287,6 @@ int evolve_system_binary_c(struct stardata_t * RESTRICT const stardata)
*/
Evolution_logging;
-
/*
* Update timestep triggers
*/
diff --git a/src/setup/cmd_line_args_list.h b/src/setup/cmd_line_args_list.h
index bff99ffe9..2f959e2d0 100644
--- a/src/setup/cmd_line_args_list.h
+++ b/src/setup/cmd_line_args_list.h
@@ -952,6 +952,15 @@ BATCHMODE_ARGS
Var(stardata->preferences->disable_events),
1.0
CMD_LINE_ARG_T_OPTIONAL_ARGS },
+{
+ ARG_SECTION_STARS,
+ "evolution_splitting_sn_eccentricity_threshold",
+ "Threshold eccentricity above which evolution splitting happens in a system with no SN kick. (0.01)",
+ ARG_FLOAT ,
+ WTTS_USE_DEFAULT,
+ EVOLUTION_SPLITTING_SN_ECCENTRICITY_THRESHOLD_VAR,
+ 1.0
+ CMD_LINE_ARG_T_OPTIONAL_ARGS },
{
ARG_SECTION_STARS,
diff --git a/src/setup/cmd_line_macros.h b/src/setup/cmd_line_macros.h
index 25f47570b..8837d55a4 100644
--- a/src/setup/cmd_line_macros.h
+++ b/src/setup/cmd_line_macros.h
@@ -298,10 +298,12 @@
#define EVOLUTION_SPLITTING_VAR Var(stardata->preferences->evolution_splitting)
#define EVOLUTION_SPLITTING_MAXDEPTH_VAR Var(stardata->preferences->evolution_splitting_maxdepth)
#define EVOLUTION_SPLITTING_SN_N_VAR Var(stardata->preferences->evolution_splitting_sn_n)
+#define EVOLUTION_SPLITTING_SN_ECCENTRICITY_THRESHOLD_VAR Var(stardata->preferences->evolution_splitting_sn_eccentricity_threshold)
#else
#define EVOLUTION_SPLITTING_VAR Not_used_arg
#define EVOLUTION_SPLITTING_SN_N_VAR Not_used_arg
#define EVOLUTION_SPLITTING_MAXDEPTH_VAR Not_used_arg
+#define EVOLUTION_SPLITTING_SN_ECCENTRICITY_THRESHOLD_VAR Not_used_arg
#endif
#ifdef DISCS
diff --git a/src/setup/set_default_preferences.c b/src/setup/set_default_preferences.c
index 55b634bcb..4fd7e5ab3 100644
--- a/src/setup/set_default_preferences.c
+++ b/src/setup/set_default_preferences.c
@@ -18,7 +18,8 @@ void set_default_preferences(struct preferences_t * RESTRICT const preferences)
*
* to list the default parameters.
*
- * Beware that these defaults are not
+ * The preferences struct is allocated using Calloc,
+ * so everything is zero by default.
*/
timestep_set_default_multipliers(preferences);
#ifdef TIMESTEP_MODULATION
@@ -427,6 +428,7 @@ preferences->WD_accretion_rate_novae_upper_limit_other_donor = ACCRETION_RATE_NO
#ifdef EVOLUTION_SPLITTING
preferences->evolution_splitting_sn_n = EVOLUTION_SPLITTING_SUPERNOVA_N;
preferences->evolution_splitting_maxdepth = EVOLUTION_SPLITTING_MAX_SPLITDEPTH_DEFAULT;
+ preferences->evolution_splitting_sn_eccentricity_threshold = 0.01;
#endif // EVOLUTION_SPLITTING
#ifdef DISCS
#if defined DISC_LOG || defined DISC_LOG_2D
diff --git a/src/supernovae/set_kick_velocity.c b/src/supernovae/set_kick_velocity.c
index 468c84e87..a48cacbe0 100644
--- a/src/supernovae/set_kick_velocity.c
+++ b/src/supernovae/set_kick_velocity.c
@@ -35,13 +35,12 @@ void set_kick_velocity(struct stardata_t * const stardata,
* Use Henon's method for pairwise components (Douglas Heggie 22/5/97).
*/
#if defined TEST_KICK_SPEED_DISTRIBUTION || defined TEST_KICK_DIRECTIONS
- int il;
- for(il=0;il<NTESTS;il++)
+ for(int il=0;il<NTESTS;il++)
{
#endif
double u1=0.0,v[5]; // v is the kick velocity vector
unsigned int k;
- for (k=1; k<=2; k++)
+ for (k=1; k<3; k++)
{
u1 = (double)random_number(stardata,random_seed);
double u2 = (double)random_number(stardata,random_seed);
diff --git a/src/supernovae/supernova_should_split_evolution.c b/src/supernovae/supernova_should_split_evolution.c
new file mode 100644
index 000000000..37eeadaa1
--- /dev/null
+++ b/src/supernovae/supernova_should_split_evolution.c
@@ -0,0 +1,71 @@
+#include "../binary_c.h"
+No_empty_translation_unit_warning;
+
+#include "supernovae.h"
+
+Boolean supernova_should_split_evolution(struct stardata_t * const stardata)
+{
+ /*
+ * This function returns TRUE if evolution splitting
+ * should happen. This is based on a possibly complicated
+ * set of logical decisions.
+ *
+ * Note: we must check both stars for supernovae,
+ * and we should use the previous_stardata to
+ * check the orbit (e.g. eccentricity, duplicity)
+ * because the stardata we have now has already
+ * had a kick if there was a suitable supernova.
+ */
+
+ /*
+ * We don't need to split single stars: their
+ * orbit cannot be sampled and their kick
+ * speed is (probably) irrelevant.
+ */
+ if(System_was_binary)
+ {
+ Foreach_star(star)
+ {
+ if(star->SN_type != SN_NONE)
+ {
+ /*
+ * If the orbit is eccentric, we should split to
+ * sample different anomalies regardless of the
+ * kick speed.
+ */
+ if(stardata->previous_stardata->common.orbit.eccentricity >
+ stardata->preferences->evolution_splitting_sn_eccentricity_threshold)
+ {
+ return TRUE;
+ }
+
+ /*
+ * If there is a Maxwellian kick, with non-zero dispersion,
+ * splitting should happen because the kick speed is random.
+ */
+ else if(
+ stardata->preferences->sn_kick_distribution[star->SN_type] == KICK_VELOCITY_MAXWELLIAN
+ &&
+ Is_not_zero(stardata->preferences->sn_kick_dispersion[star->SN_type])
+ )
+ {
+ return TRUE;
+ }
+
+ /*
+ * If the kick is a custom function, we should split just
+ * in case.
+ */
+ else if(stardata->preferences->sn_kick_distribution[star->SN_type] == KICK_VELOCITY_CUSTOM)
+ {
+ return TRUE;
+ }
+ }
+ }
+ }
+
+ /*
+ * ... no criteria are satisfied, so no splitting is required.
+ */
+ return FALSE;
+}
diff --git a/src/supernovae/supernovae_prototypes.h b/src/supernovae/supernovae_prototypes.h
index 6a3626bc2..b6e7c5b4f 100644
--- a/src/supernovae/supernovae_prototypes.h
+++ b/src/supernovae/supernovae_prototypes.h
@@ -97,4 +97,6 @@ Yield_source supernova_source_id(struct star_t * exploder,
struct star_t * pre_explosion_star,
struct star_t * companion);
+Boolean supernova_should_split_evolution(struct stardata_t * const stardata);
+
#endif // SUPERNOVA_PROTOTYPES_H
diff --git a/tbse b/tbse
index 6efea3001..d9afcb915 100755
--- a/tbse
+++ b/tbse
@@ -1316,10 +1316,14 @@ SN_KICK_COMPANION_IA_Hybrid_HeCOWD_subluminous=SN_IMPULSE_NONE
# when a SN occurs.
# (If this is 1, then the evolution is as normal.)
#
+# The evolution splitting algorithm has code to detect whether
+# a split should happen when there is a supernova. You have to choose
+# an eccentricity above which splitting happens in
+# evolution_splitting_sn_eccentricity_threshold
EVOLUTION_SPLITTING=False
EVOLUTION_SPLITTING_MAXDEPTH=2
EVOLUTION_SPLITTING_SN_N=10
-
+EVOLUTION_SPLITTING_SN_ECCENTRICITY_THRESHOLD=0.01
##############################
# WD kick distribution velocity dispersion
@@ -2355,6 +2359,7 @@ $WARMUP \
--evolution_splitting $EVOLUTION_SPLITTING \
--evolution_splitting_sn_n $EVOLUTION_SPLITTING_SN_N \
--evolution_splitting_maxdepth $EVOLUTION_SPLITTING_MAXDEPTH \
+--evolution_splitting_sn_eccentricity_threshold $EVOLUTION_SPLITTING_SN_ECCENTRICITY_THRESHOLD \
--RLOF_angular_momentum_transfer_model $RLOF_ANGULAR_MOMENTUM_TRANSFER_MODEL \
--gravitational_radiation_model $GRAVITATIONAL_RADIATION_MODEL \
--gravitational_radiation_modulator_J $GRAVITATIONAL_RADIATION_MODULATOR_J \
--
GitLab