From dbbae12819a7142782d16b60107d95636b315224 Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Thu, 7 Nov 2019 14:41:59 +0000
Subject: [PATCH] add individual nova timestep acceleration
---
src/binary_c_derivatives.h | 7 ++--
src/binary_c_parameters.h | 2 +-
src/binary_c_structures.h | 7 ++++
...ar_mass_and_angular_momentum_derivatives.c | 10 +++++
.../limit_accretion_rates.c | 24 ++++++++---
src/logging/log_every_timestep.c | 6 ++-
src/perl/scripts2/template.pl | 2 +-
src/perl/scripts2/yield_vs_time.pl | 6 +--
src/setup/cmd_line_args.h | 32 +++++++++++++-
src/setup/set_default_preferences.c | 5 +++
src/timestep/timestep_limits.c | 42 ++++++++++++++++++-
tbse | 35 ++++++++++++++++
12 files changed, 160 insertions(+), 18 deletions(-)
diff --git a/src/binary_c_derivatives.h b/src/binary_c_derivatives.h
index 6987a8eb2..e801bd212 100644
--- a/src/binary_c_derivatives.h
+++ b/src/binary_c_derivatives.h
@@ -102,9 +102,9 @@
#define DERIVATIVE_STELLAR_NUM_THERMAL_PULSES_SINCE_MCMIN 41
#define DERIVATIVE_STELLAR_H_LAYER_MASS 42
#define DERIVATIVE_STELLAR_ROCHE_RADIUS 43
+#define DERIVATIVE_STELLAR_NUM_NOVAE 44
-
-#define NUMBER_OF_STELLAR_DERIVATIVES 44
+#define NUMBER_OF_STELLAR_DERIVATIVES 45
#define STELLAR_DERIVATIVE_STRINGS \
{ \
@@ -151,7 +151,8 @@
"Number of thermal pulses", \
"Number of thermal pulses since mcmin", \
"nova H layer mass", \
- "Roche radius" \
+ "Roche radius", \
+ "Number of novae" \
}
#define Show_stellar_derivatives(S,SEPARATOR) \
diff --git a/src/binary_c_parameters.h b/src/binary_c_parameters.h
index ef93466ad..2c67ee6d9 100644
--- a/src/binary_c_parameters.h
+++ b/src/binary_c_parameters.h
@@ -1690,7 +1690,7 @@
//#define TPAGBTRACKS
/* log for Window To The Stars 2 */
-#define WTTS_LOG
+//#define WTTS_LOG
/* anti TZ project */
//#define ANTI_TZ
diff --git a/src/binary_c_structures.h b/src/binary_c_structures.h
index fe4f0c4c6..5a4552cf5 100644
--- a/src/binary_c_structures.h
+++ b/src/binary_c_structures.h
@@ -570,6 +570,10 @@ struct preferences_t {
double beta_reverse_nova;
double nova_irradiation_multiplier;
double nova_faml_multiplier;
+
+ double nova_timestep_accelerator_num;
+ double nova_timestep_accelerator_index;
+ double nova_timestep_accelerator_max;
// gamma is the angular momentum factor for mass lost during Roche (-1.0)
double rlof_angmom_gamma; // 0.0
@@ -1331,6 +1335,7 @@ struct star_t {
double num_thermal_pulses;
double num_thermal_pulses_since_mcmin;
+ double num_novae;
double dntp;
double menv_1tp;
double Mc_prev_pulse;
@@ -1741,6 +1746,8 @@ struct model_t {
Boolean split_last_time;
#endif // EVOLUTION_SPLITTING
+
+ double nova_aggression_factor;
};
struct orbit_t {
diff --git a/src/binary_star_functions/apply_stellar_mass_and_angular_momentum_derivatives.c b/src/binary_star_functions/apply_stellar_mass_and_angular_momentum_derivatives.c
index 7470fc2e7..130b5d288 100644
--- a/src/binary_star_functions/apply_stellar_mass_and_angular_momentum_derivatives.c
+++ b/src/binary_star_functions/apply_stellar_mass_and_angular_momentum_derivatives.c
@@ -251,7 +251,17 @@ void apply_stellar_mass_and_angular_momentum_derivatives(struct stardata_t * RES
star->derivative,
dt,
DERIVATIVE_STELLAR_H_LAYER_MASS);
+
}
+
+ /*
+ * Count novae
+ */
+ apply_derivative(stardata,
+ &star->num_novae,
+ star->derivative,
+ dt,
+ DERIVATIVE_STELLAR_NUM_NOVAE);
/*
* Update the stellar type as a result of accretion:
diff --git a/src/binary_star_functions/limit_accretion_rates.c b/src/binary_star_functions/limit_accretion_rates.c
index 74041bfdf..c8ed7d296 100644
--- a/src/binary_star_functions/limit_accretion_rates.c
+++ b/src/binary_star_functions/limit_accretion_rates.c
@@ -336,7 +336,14 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
1.0 / (q*(1.0 + q)) / accretor->mass;
-
+ /*
+ * Calculate the recurrence time, i.e. the time
+ * between successive nova explosions
+ */
+ const double recurrence_time =
+ nova_recurrence_time(accretion_rate,
+ accretor->mass);
+
/*
* Now we have calculated the properties of the
* novae, decide whether we should model them
@@ -372,7 +379,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
* Increase the surface hydrogen layer mass.
*/
accretor->derivative[DERIVATIVE_STELLAR_H_LAYER_MASS] = accretion_rate;
-
+
/*
* Check if we have sufficient mass for explosion
*/
@@ -416,8 +423,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
stardata->model.nova_timeout = stardata->model.time +
NOVA_TIMEOUT_FACTOR * // usually 2
1e-6 * // convert yr to Myr
- nova_recurrence_time(accretion_rate,
- accretor->mass);
+ recurrence_time;
}
else
{
@@ -458,7 +464,15 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
accretor->derivative[DERIVATIVE_STELLAR_MASS_NOVA]);
}
-
+
+
+ /*
+ * Count novae
+ */
+ accretor->derivative[DERIVATIVE_STELLAR_NUM_NOVAE] =
+ 1.0 / recurrence_time;
+
+
/*
donor->derivative[DERIVATIVE_STELLAR_MASS_IRRADIATIVE_LOSS] =
- 1.1e18 * (xi/4.0) * pow(1.0 + donor->mass,2.0/3.0) * pow(donor->mass,2.833)
diff --git a/src/logging/log_every_timestep.c b/src/logging/log_every_timestep.c
index 06628e436..8ed96f655 100644
--- a/src/logging/log_every_timestep.c
+++ b/src/logging/log_every_timestep.c
@@ -2177,7 +2177,7 @@ void log_every_timestep(struct stardata_t * RESTRICT const stardata)
#ifdef LOG_INDIVIDUAL_NOVAE
if(WHITE_DWARF(stardata->star[0].stellar_type))
{
- printf("NOVAE %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %u %u %d %20.12e %20.12e %20.12e %g %g %20.12e %d %d\n",
+ printf("NOVAE %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %u %u %d %20.12e %20.12e %20.12e %g %g %20.12e %d %d %g %g\n",
stardata->model.time,//1
stardata->star[0].mass,
dt,//3
@@ -2197,7 +2197,9 @@ void log_every_timestep(struct stardata_t * RESTRICT const stardata)
stardata->star[1].dmacc,//17
stardata->model.nova_timeout, //18
stardata->star[0].dtlimiter,//19
- stardata->star[1].dtlimiter //20
+ stardata->star[1].dtlimiter, //20
+ stardata->star[0].num_novae, //21
+ stardata->star[1].num_novae //22
);
}
#endif//LOG_INDIVIDUAL_NOVAE
diff --git a/src/perl/scripts2/template.pl b/src/perl/scripts2/template.pl
index 29e8a6073..0194e8883 100755
--- a/src/perl/scripts2/template.pl
+++ b/src/perl/scripts2/template.pl
@@ -123,7 +123,7 @@ if($duplicity == 0)
{
# make a grid of $nstars single binary stars, log-spaced,
# with masses between $mmin and $mmax
- my $nstars = 1000;
+ my $nstars = 10000;
my $mmin = 0.1;
my $mmax = 80.0;
$population->add_grid_variable(
diff --git a/src/perl/scripts2/yield_vs_time.pl b/src/perl/scripts2/yield_vs_time.pl
index 61dd8ff85..11030e5a4 100755
--- a/src/perl/scripts2/yield_vs_time.pl
+++ b/src/perl/scripts2/yield_vs_time.pl
@@ -19,7 +19,7 @@ use Sys::Hostname;
$|=1; # enable this line for auto-flushed output
############################################################
-# Version 2.1 for binary_c 2.1
+# Version 2.1.3 for binary_c 2.1.3
############################################################
my $args = "@ARGV"; # used often
@@ -28,7 +28,7 @@ my $hostname = hostname();
# output directory location
my $voldisk = $ENV{VOLDISK} // $ENV{HOME};
my $outdir = ($args=~/outdir=(\S+)/)[0] //
- $voldisk.'/data/pops/yields_vs_time-SNIa';
+ $voldisk.'/data/pops/yields_vs_time-05112019';
print "Yield vs time : outdir=$outdir : voldisk=$voldisk, outdir=$outdir\n";
use vars qw($dt $logdt $tstart $logtimes $condor_count $differential $distribution $distribution_options $ensemble %ensemble_opts @ensemble_label $exacttimes $force $gnuplot $gridstring @isotopes $idt $nbins $nensemble $nisotopes $nensemble1 $nisotopes1 $norm_output $nsources $sources @source_label $version_string $yieldoutfile );
@@ -142,7 +142,7 @@ sub parser
}
else
{
- print "$header @$linearray\n";
+ #print "$header @$linearray\n";
}
# other output is ignored
}
diff --git a/src/setup/cmd_line_args.h b/src/setup/cmd_line_args.h
index 3ddbd8dea..aa4fefd67 100644
--- a/src/setup/cmd_line_args.h
+++ b/src/setup/cmd_line_args.h
@@ -2994,7 +2994,7 @@ struct cmd_line_arg_t
}, \
\
{ \
- ARG_SECTION_OUTPUT, \
+ ARG_SECTION_BINARY, \
"individual_novae", \
"If individual_novae is True, novae are resolved such that each explosion is performed separtaely.", \
ARG_BOOLEAN, \
@@ -3003,6 +3003,36 @@ struct cmd_line_arg_t
1.0 \
}, \
\
+{ \
+ ARG_SECTION_BINARY, \
+ "nova_timestep_accelerator_num", \
+ "The nova timestep is accelerated if the nova number exceeds nova_timestep_accelerator_num. If zero or negative, acceleration is off. See also nova_timestep_accelerator_index and nova_timestep_accelerator_max.\n", \
+ ARG_FLOAT, \
+ WTTS_USE_DEFAULT, \
+ &(stardata->preferences->nova_timestep_accelerator_num), \
+ 1.0 \
+ }, \
+ \
+{ \
+ ARG_SECTION_BINARY, \
+ "nova_timestep_accelerator_index", \
+ "The index at which the nova timestep is accelerated. A larger value gives longer timesteps. See also nova_timestep_accelerator_num and nova_timestep_accelerator_max.\n", \
+ ARG_FLOAT, \
+ WTTS_USE_DEFAULT, \
+ &(stardata->preferences->nova_timestep_accelerator_index), \
+ 1.0 \
+ }, \
+ \
+{ \
+ ARG_SECTION_BINARY, \
+ "nova_timestep_accelerator_max", \
+ "The nova timestep is accelerated by a factor that is capped at nova_timestep_accelerator_max. This parameter is ignored if it is zero or negative. See also nova_timestep_accelerator_num and nova_timestep_accelerator_index.\n", \
+ ARG_FLOAT, \
+ WTTS_USE_DEFAULT, \
+ &(stardata->preferences->nova_timestep_accelerator_max), \
+ 1.0 \
+ }, \
+ \
{ \
ARG_SECTION_OUTPUT, \
"force_corotation_on_RLOF", \
diff --git a/src/setup/set_default_preferences.c b/src/setup/set_default_preferences.c
index 9daf72989..d4ff10033 100644
--- a/src/setup/set_default_preferences.c
+++ b/src/setup/set_default_preferences.c
@@ -87,6 +87,11 @@ void set_default_preferences(struct preferences_t * RESTRICT const preferences)
preferences->donor_limit_thermal_multiplier=DONOR_LIMIT_THERMAL_MULTIPLIER_DEFAULT;
preferences->donor_limit_dynamical_multiplier=DONOR_LIMIT_DYNAMICAL_MULTIPLIER_DEFAULT;
preferences->donor_limit_dynamical_multiplier=DONOR_LIMIT_ENVELOPE_MULTIPLIER_DEFAULT;
+
+ /* nova timestep accelerators */
+ preferences->nova_timestep_accelerator_num = 100;
+ preferences->nova_timestep_accelerator_max = -1.0;
+ preferences->nova_timestep_accelerator_index = 0.5;
preferences->rlof_angmom_gamma=NONCONSERVATIVE_ANGMOM_GAMMA_DEFAULT;
preferences->nonconservative_angmom_gamma=NONCONSERVATIVE_ANGMOM_GAMMA_DEFAULT;
diff --git a/src/timestep/timestep_limits.c b/src/timestep/timestep_limits.c
index 1848c2537..e394f72ec 100644
--- a/src/timestep/timestep_limits.c
+++ b/src/timestep/timestep_limits.c
@@ -283,14 +283,52 @@ void timestep_limits(Timestep_prototype_args)
stardata->preferences->timestep_multipliers[DT_LIMIT_NOVAE] *
dtnova;
+ /*
+ * Modulate to combine multiple novae
+ *
+ * The factor f is applied only if we have already had more than
+ * nova_timestep_accelerator_num novae.
+ *
+ * If accelerator_num <= 0, f = 1.0.
+ *
+ * If Nnovae < accelerator_num, f = 1.0.
+ *
+ * Then, if acclerator_max > 0.0 we have:
+ *
+ * f = Max(accelerator_max, pow(Nnovae, accelerator_index))
+ *
+ * otherwise
+ *
+ * f = pow(Nnovae, accelerator_index).
+ */
+ const double f =
+
+ stardata->preferences->nova_timestep_accelerator_num < +TINY ?
+ /* num < 0.0 : f = 1.0 */
+ 1.0 :
+ (
+ /* the cap is LARGE_FLOAT (i.e. very large) if acclerator_max <= 0.0 */
+ Max(
+ (stardata->preferences->nova_timestep_accelerator_max < +TINY ?
+ LARGE_FLOAT :
+ stardata->preferences->nova_timestep_accelerator_max),
+
+ /* check we've had sufficient novae */
+ (star->num_novae > stardata->preferences->nova_timestep_accelerator_num ?
+ /* apply the pow(Nnovae, index) */
+ pow(star->num_novae,stardata->preferences->nova_timestep_accelerator_index) :
+ /* or don't */
+ 1.0)
+ )
+ );
+
/*
* limit timestep because of novae
*/
Limit_timestep(*dt,
- dtlim,
+ dtlim*f,
star,
DT_LIMIT_NOVAE);
-
}
if(star->starnum == 0)
diff --git a/tbse b/tbse
index ce154f496..351de6b6b 100755
--- a/tbse
+++ b/tbse
@@ -1181,6 +1181,37 @@ TIMESTEP_SOLVER_FACTOR=1.0
# global timestep modulator for non-remnants (1.0)
TIMESTEP_MODULATOR=1.0
+############################################################
+# Nova timestep acceleration.
+#
+# The idea here
+#
+# The factor f is applied only if we have already had more than
+# nova_timestep_accelerator_num novae.
+#
+# If accelerator_num <= 0, f = 1.0 (i.e. no acceleration).
+#
+# If Nnovae < accelerator_num, f = 1.0 (i.e. no acceleration).
+#
+# Then, if acclerator_max > 0.0 we use,
+#
+# f = Max(accelerator_max, pow(Nnovae, accelerator_index))
+#
+# otherwise
+#
+# f = pow(Nnovae, accelerator_index).
+#
+# Typically, num=100, index=0.5 and max=-1 works well.
+#
+# Note that max=-1 does not mean the timestep will be necessarily
+# large: there are many other limits that can (and do) apply (see
+# below).
+#
+############################################################
+NOVA_TIMESTEP_ACCELERATOR_NUM=100 # 100
+NOVA_TIMESTEP_ACCELERATOR_INDEX=0.5 # 0.5
+NOVA_TIMESTEP_ACCELERATOR_MAX=-1 # -1
+
############################################################
# Individual by-process timestep multipliers
#
@@ -2142,6 +2173,10 @@ $TIMESTEP_MULTIPLIERS \
--overspin_algorithm $OVERSPIN_ALGORITHM \
--gaia_Teff_binwidth $GAIA_TEFF_BINWIDTH \
--gaia_L_binwidth $GAIA_L_BINWIDTH \
+--nova_timestep_accelerator_num $NOVA_TIMESTEP_ACCELERATOR_NUM \
+--nova_timestep_accelerator_index $NOVA_TIMESTEP_ACCELERATOR_INDEX \
+--nova_timestep_accelerator_max $NOVA_TIMESTEP_ACCELERATOR_MAX \
+
";
# valgrind options
--
GitLab