diff --git a/src/binary_c_derivatives.h b/src/binary_c_derivatives.h
index 6987a8eb2aa4f75af487a6d5fc652321f60ef1b5..e801bd2129507bcce9e79700fe816cd50ed081ec 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 ef93466ad75e5a9b7e3bb1ccbc633722f5349162..2c67ee6d95d4343a25aa90a74b47ff204f802ffc 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 fe4f0c4c63f032a6b9031a8c0d8cd57623e36511..5a4552cf5b3060ddc95431393d1a2e196bf84f04 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 7470fc2e7246af685b16c50761aaa55e8fb91944..130b5d2885ce544c7bbd8a1660df43a88f49b6fa 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 74041bfdf03bb103fc2d9f5023d545514e2dcf8d..c8ed7d2965326f49a38d0e12d6b32048266ce64d 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 06628e4367203fe7589d056c0d23e7a1916ac393..8ed96f655708eac5ec7745019370e913fe5ee4d1 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 29e8a607334ac371e917850ce10fe96d542218e6..0194e8883b55b5e7b9e8802b7876ef7a3ec16f80 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 61dd8ff85f21b5cf74190a9f8e307ab0890afcf3..11030e5a46ef72558ecbeb8d3836207402727ff0 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 3ddbd8dea1329e242af2fecc2875dfc32a706087..aa4fefd67fbffc000cbe48be45290b8f0d810c66 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 9daf72989def8cb7a7a35d6f49673c750b572eb3..d4ff10033ee7153284986cc2998eaf1465459124 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 1848c2537aa8f517e10d0f4f084d1d7c9eb7f171..e394f72ec2d47a64d64aec77d90f9cc7717fb839 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 ce154f496c4f1391f4ab06f69b2af16925c9fd53..351de6b6b48db9f10ce1e9c8ddd26d6803d1f7ef 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