diff --git a/src/binary_c_macros.h b/src/binary_c_macros.h
index c24b6d2d1c2041cc70a4564a0ec24ca4a7f27183..2eb41d9d4a8e7e8fb0fdad92394501a66ef3004f 100644
--- a/src/binary_c_macros.h
+++ b/src/binary_c_macros.h
@@ -837,18 +837,23 @@
}
+/*
+ * Nova states
+ */
#define NOVA_STATE_NONE 0
-#define NOVA_STATE_TRIGGER 1
-#define NOVA_STATE_POST 2
+#define NOVA_STATE_PRE_EXPLOSION 1
+#define NOVA_STATE_EXPLODING 2
#define NOVA_STATE_QUIESCENT 3
#define NOVA_STATE_STRINGS { \
"None", \
- "Triggered", \
- "Post-nova", \
+ "Pre-explosion", \
+ "Exploding", \
"Quiescent" \
}
-
+/*
+ * TPAGB luminosity algorithm
+ */
#define TPAGB_LUMINOSITY_AVERAGE 0
#define TPAGB_LUMINOSITY_DIPS 1
diff --git a/src/binary_c_parameters.h b/src/binary_c_parameters.h
index 634e1c7b2aff255a0576cc14170c9160b6baf85f..5e7d96474a5d4c9cbdbe2abf331343ceed453237 100644
--- a/src/binary_c_parameters.h
+++ b/src/binary_c_parameters.h
@@ -720,6 +720,12 @@
/* deprecated: instant mixing on the companion */
//#define NOVA_BACK_ACCRETION_INSTANT_MIXING
+/*
+ * after a nova explosion, we continue short timesteps
+ * until NOVA_TIMEOUT_FACTOR * nova_recurrance_time
+ */
+#define NOVA_TIMEOUT_FACTOR 2.0
+
/************************************************************
* Mass transfer by stellar winds
************************************************************/
diff --git a/src/binary_c_structures.h b/src/binary_c_structures.h
index 21057d23ca6a8308bd5b3835439fc5fae0255566..28e553a382964290646a5792ceb6d9d22ea6f035 100644
--- a/src/binary_c_structures.h
+++ b/src/binary_c_structures.h
@@ -1054,10 +1054,10 @@ struct diffstats_t
double orbital_eccentricity;
double orbital_angular_frequency;
double novarate[NUMBER_OF_STARS],v_eq_ratio[NUMBER_OF_STARS];
- int novastate[NUMBER_OF_STARS];
- double time;
+ Time time;
Stellar_type stellar_type[NUMBER_OF_STARS],hybrid_HeCOWD[NUMBER_OF_STARS];
Supernova_type SN_type[NUMBER_OF_STARS];
+ Nova_type nova[NUMBER_OF_STARS];
Boolean sgl,artificial_accretion;
};
@@ -1359,9 +1359,7 @@ struct star_t {
double interpulse_period;
double mc_1tp; /* core mass at first pulse */
double lambda_3dup;
- double time_prev_nova,time_next_nova,dtnova;
-
- double dm_novaH,nova_timeout;
+ double dm_novaH;
double nova_beta,nova_faml;
#ifdef RUNAWAY_STARS
@@ -1503,7 +1501,6 @@ struct star_t {
Stellar_type stellar_type; /* stellar type */
Stellar_type detached_stellar_type; /* stellar type when detached */
Star_number starnum; // set to 0 for primary, 1 for secondary. useful for logging
- int novastate;
#ifdef DISCS
int ndiscs;
@@ -1534,6 +1531,7 @@ struct star_t {
Reject_index reject;
Stellar_type compact_core_type;
+ Nova_state nova;
Boolean hybrid_HeCOWD;
#ifdef WD_KICKS
@@ -1543,7 +1541,7 @@ struct star_t {
#ifdef LOG_SUPERNOVAE
Boolean supernova;
#endif//LOG_SUPERNOAVE
- Boolean nova;
+
#ifdef NUCSYN
Boolean rstar; /* true if star is R type */
Boolean first_dredge_up; /* stores whether we've had first dredge up */
@@ -1632,7 +1630,7 @@ struct model_t {
double time; /* the time, was tphys in the fortran version */
double max_evolution_time; /* Maximum evolution time, in Myr */
double dt,dtm,true_dtm,time_remaining,dtmwas; /* timestepping */
- double dt_zoomfac,RLOF_recommended_timestep;
+ double dt_zoomfac,RLOF_recommended_timestep,nova_timeout;
double prev_log_t;
#ifdef FABIAN_IMF_LOG
double next_fabian_imf_log_time;
diff --git a/src/binary_c_types.h b/src/binary_c_types.h
index a0cdbecac0d4c03de9b9971b55e0e60cfae9a735..6134ca5d347f3128383f88e15dceec57ca60d42f 100644
--- a/src/binary_c_types.h
+++ b/src/binary_c_types.h
@@ -123,6 +123,7 @@ typedef unsigned int Boolean;
#define Spectral_type __int__
#define Long_spectral_type __double__
#define Nova_type __unsigned__ __int__
+#define Nova_state __unsigned__ __int__
#define Reaction_network __unsigned__ __int__
#ifdef CODESTATS
diff --git a/src/binary_star_functions/limit_accretion_rates.c b/src/binary_star_functions/limit_accretion_rates.c
index eadbe7c992748068e8a0a38cd679d3a4617275fe..74041bfdf03bb103fc2d9f5023d545514e2dcf8d 100644
--- a/src/binary_star_functions/limit_accretion_rates.c
+++ b/src/binary_star_functions/limit_accretion_rates.c
@@ -20,8 +20,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
{
struct star_t * accretor = &stardata->star[k];
struct star_t * donor = &stardata->star[Other_star(k)];
- accretor->nova = FALSE;
-
+
Dprint("limit acc rate star %d : is currently %g\n",
k,
Mdot_net(accretor));
@@ -240,7 +239,6 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
{
steady_burn_rate = steady_burn_rate_He;
nova_type = NOVA_TYPE_HELIUM;
- accretor->nova = TRUE;
}
/* restore stellar type */
@@ -337,6 +335,14 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
Pythag2(vrat,1.0) * geometric_accretion_fraction *
1.0 / (q*(1.0 + q)) / accretor->mass;
+
+
+ /*
+ * Now we have calculated the properties of the
+ * novae, decide whether we should model them
+ * individually (with short timesteps), or in an
+ * average way (with long timesteps).
+ */
if(stardata->preferences->individual_novae == TRUE)
{
/*
@@ -382,7 +388,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
new_nova->accretion_rate = accretion_rate;
new_nova->steady_burn_rate = steady_burn_rate;
new_nova->f = f;
- accretor->nova = TRUE;
+ accretor->nova = NOVA_STATE_EXPLODING;
if(Add_new_event(stardata,
BINARY_C_EVENT_NOVA,
&nova_event_handler,
@@ -403,95 +409,34 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
*/
Safe_free(new_nova);
}
- }
-#ifdef __DEPRECATED
- /*
- * The following code was for before we used
- * the events framework for individiual novae.
- */
-
- /*
- * Check if we have sufficient mass for explosion
- */
- if(accretor->dm_novaH > dMcrit)
- {
+
/*
- * Nova explosion :
- *
- * In theory, we should burn dMcrit of material.
- * However, there is likely a bit more than dMcrit
- * in the layer. We burn it all to maintain mass
- * conservation.
+ * If we're exploding now, set the next timeout
*/
- accretor->nova = TRUE;
- if(accretor->novastate == NOVA_STATE_TRIGGER)
- {
- /*
- * A nova has been triggered.
- *
- * The timestep is short.
- *
- * Calculate dm_lost, the mass lost from
- * the white dwarf in the nova explosion,
- * hence the nova mass loss rate.
- */
-
- const double dm_lost = (1.0-f) * accretor->dm_novaH;
- const double dm_loss_rate = - dm_lost / stardata->model.dt;
- accretor->derivative[DERIVATIVE_STELLAR_MASS_NOVA] +=
- dm_loss_rate;
-
- Dprint("remove %g from M=%g (f = %g) at rate %g with dt %g \n",
- dm_lost,
- accretor->mass,
- f,
- accretor->derivative[DERIVATIVE_STELLAR_MASS_NOVA],
- stardata->model.dt
- );
-
- /*
- * Remove all the accreted layer : anything
- * remaining becomes part of the white dwarf.
- */
- accretor->derivative[DERIVATIVE_STELLAR_H_LAYER_MASS] +=
- dm_loss_rate;
-
- /*
- * Trigger complete, change novastate state
- * to NOVA_STATE_POST which enforces another
- * short timestep with no nova.
- */
- accretor->novastate = NOVA_STATE_POST;
- }
- else
- {
- /*
- * trigger a nova on next (short) timestep
- */
- accretor->novastate = NOVA_STATE_TRIGGER;
- }
+ 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);
}
-#endif // __DEPRECATED (was INDIVIDUAL_NOVAE_EVENTS)
-
else
{
/*
* Phase between novae.
*
* If there has been a previous nova explosion, i.e.
- * novastate is not NOVA_STATE_NONE (it is likely
- * NOVA_STATE_POST), then we set NOVA_STATE_QUIESCENT.
+ * accretor->nova is not NOVA_STATE_NONE,
+ * then we set accretor->nova to NOVA_STATE_QUIESCENT.
*
- * If there has never been a nova, stick to
- * the current state which must be NOVA_STATE_NONE.
+ * If there has not yet been a nova, but we expect one
+ * soon, set to NOVA_STATE_PRE_EXPLOSION.
*/
accretor->derivative[DERIVATIVE_STELLAR_MASS_NOVA] = 0.0;
-
- if(accretor->novastate != NOVA_STATE_NONE)
- {
- accretor->nova = TRUE;
- accretor->novastate = NOVA_STATE_QUIESCENT;
- }
+ accretor->nova =
+ (accretor->nova == NOVA_STATE_QUIESCENT ||
+ accretor->nova == NOVA_STATE_EXPLODING) ?
+ NOVA_STATE_QUIESCENT :
+ NOVA_STATE_PRE_EXPLOSION;
}
}
else
@@ -499,11 +444,15 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
/*
* Novae are modelled as average events, which are put in
- * the normal mass loss
+ * the normal mass loss.
+ *
+ * We set nova to be NOVA_STATE_EXPLODING on the assumption
+ * that the timestep is long so there are going to be some
+ * explosions, but it's not really relevant.
*/
accretor->derivative[DERIVATIVE_STELLAR_MASS_NOVA] +=
- (1.0-f) * accretion_rate;
- accretor->nova = TRUE;
+ accretor->nova = NOVA_STATE_EXPLODING;
Dprint("accretion rate %g -> dM/dt nova = %g\n",
accretion_rate,
accretor->derivative[DERIVATIVE_STELLAR_MASS_NOVA]);
@@ -531,11 +480,10 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
donor->derivative[DERIVATIVE_STELLAR_MASS_NOVA] +=
accretor->nova_beta * -accretor->derivative[DERIVATIVE_STELLAR_MASS_NOVA];
}
-
else
{
/*
- * We have steady initial burning.
+ * We have steady burning.
*
* Accretion rate exceeds the steady burning rate
* either burning is steady, or the star forms
@@ -544,9 +492,9 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
* Note: we do not treat the case where
* there is steady accretion of hydrogen, but then
* helium flashes underneath which would leave to
- * helium novae. These should be caught here.
+ * helium novae. These should be caught here (todo!).
*/
- accretor->novastate = NOVA_STATE_NONE;
+ accretor->nova = NOVA_STATE_NONE;
if(accretion_rate < new_envelope_rate)
{
@@ -569,6 +517,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
/*
* Accretion onto a neutron star/black hole, what to do?
*/
+ accretor->nova = NOVA_STATE_NONE;
}
}
else
@@ -576,46 +525,29 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
/*
* No longer a compact object : cannot have novae
*/
- accretor->novastate = NOVA_STATE_NONE;
+ accretor->nova = NOVA_STATE_NONE;
}
} // accretion rate > 0 check
else
{
/*
- * We may have had novae, but are no longer transferring mass.
- * In this case, force the NOVA_STATE_QUIESCENT to prevent
- * the timestep remaining very short.
+ * Accretion rate is zero or negative, i.e. mass
+ * is constant or being lost.
*/
- accretor->novastate =
- accretor->novastate == NOVA_STATE_NONE ?
- NOVA_STATE_NONE : NOVA_STATE_QUIESCENT;
-
- if((accretor->novastate == NOVA_STATE_QUIESCENT ||
- Is_not_zero(accretor->derivative[DERIVATIVE_STELLAR_MASS_NOVA])) &&
- (Is_not_zero(accretor->nova_timeout) && stardata->model.time > accretor->nova_timeout))
- {
- accretor->nova = TRUE;
- }
+
+
+ /*
+ * If we were previously having novae, keep the state
+ * at NOVA_STATE_QUIESCENT until nova_timeout has been reached
+ */
+ accretor->nova =
+ (
+ accretor->nova == NOVA_STATE_QUIESCENT ||
+ accretor->nova == NOVA_STATE_EXPLODING ||
+ stardata->model.time < stardata->model.nova_timeout
+ )
+ ? NOVA_STATE_QUIESCENT : NOVA_STATE_NONE;
}
- }
+ }
}
-/*
- static Pure_function double MM93_k(const double m1,
- const double m2,
- const double xL)
- {
- return sqrt(
- (
- m1/Pow3(xL) +
- m2/Pow3(1.0-xL) +
- -1.0
- )
- *
- (
- m1/Pow3(xL) +
- m2/Pow3(1.0-xL)
- )
- );
- }
-*/
diff --git a/src/evolution/evolution_difflog.c b/src/evolution/evolution_difflog.c
index 3b605207ea7da2920cecf7643b9d03325705a76c..951a6d04c4c8e531f0d9e8a4d0b9b8a6324ffda4 100644
--- a/src/evolution/evolution_difflog.c
+++ b/src/evolution/evolution_difflog.c
@@ -385,7 +385,7 @@ static void set_diffstats(struct stardata_t * stardata,
now->accretion_luminosity[k] = star->accretion_luminosity;
now->luminosity[k] = star->luminosity;
now->novarate[k] = star->derivative[DERIVATIVE_STELLAR_MASS_NOVA];
- now->novastate[k] = star->nova;
+ now->nova[k] = star->nova;
now->v_eq_ratio[k] = star->v_eq_ratio;
now->SN_type[k] = star->SN_type;
now->artificial_accretion =
@@ -471,12 +471,26 @@ static void cf_diffstats(struct stardata_t * RESTRICT const stardata,
}
/* novae */
- if(now->sgl == FALSE)
{
const Boolean novae =
- stardata->preferences->individual_novae == TRUE ?
- Boolean_(now->novastate[k] == TRUE) :
- Boolean_(now->novarate[k] < 0.0);
+ now->sgl == TRUE ? FALSE :
+ (
+ stardata->preferences->individual_novae == TRUE ?
+
+ /*
+ * indiviudal novae algorithm:
+ * Novae are happening when there has been at least
+ * one explosion
+ */
+ Boolean_(now->nova[k] == NOVA_STATE_EXPLODING ||
+ now->nova[k] == NOVA_STATE_QUIESCENT) :
+ /*
+ * time-averaged novae (like BSE):
+ * Novae are happening when the nova mass derivative is
+ * negative (i.e. mass is being lost, on average).
+ */
+ Boolean_(now->novarate[k] < 0.0)
+ );
if((novae == TRUE) &&
lockflags[k][LOCK_NOVAE] == FALSE)
@@ -490,7 +504,6 @@ static void cf_diffstats(struct stardata_t * RESTRICT const stardata,
lockflags[k][LOCK_NOVAE] = FALSE;
flags[END_NOVAE_LABEL] = TRUE;
}
-
}
/* AGB star shrinkage */
@@ -636,7 +649,6 @@ static void cf_diffstats(struct stardata_t * RESTRICT const stardata,
stardata->common.Xsolar,
stardata);
CEMP = Boolean_(CFe > 1.0);
-
}
else
{
diff --git a/src/evolution/evolution_rejection_tests.c b/src/evolution/evolution_rejection_tests.c
index ceef3abe4cc95dc6284a6e00a6b1b08ecc7012a1..b30aa8b429cf028ab8ff2fb2ffa423e84b7810b5 100644
--- a/src/evolution/evolution_rejection_tests.c
+++ b/src/evolution/evolution_rejection_tests.c
@@ -1,6 +1,6 @@
#include "../binary_c.h"
static Boolean smooth_evolution(struct stardata_t * stardata);
-static Boolean smooth_evolution_checker(
+static Boolean MAYBE_UNUSED smooth_evolution_checker(
struct stardata_t * stardata,
const Boolean check_stellar_type_change,
const Boolean check_duplicity_change,
diff --git a/src/evolution/init_star.c b/src/evolution/init_star.c
index 47f185f036372cf455aee29202d3b06717005588..64ffb07d4cad4f9601f4ce18e0a4de2b95e5be1b 100644
--- a/src/evolution/init_star.c
+++ b/src/evolution/init_star.c
@@ -34,7 +34,7 @@ void init_star(struct stardata_t * const stardata,
star->max_mix_depth=1e3;
#endif /* NUCSYN */
- star->stellar_type=MAIN_SEQUENCE;
-
+ star->stellar_type = MAIN_SEQUENCE;
+ star->nova = NOVA_STATE_NONE;
Dprint("Done for this star\n");
}
diff --git a/src/logging/log_every_timestep.c b/src/logging/log_every_timestep.c
index dfee32b1ccedbfb771aa93c40d899b833f3ae5ef..1998ffe72d20628289c4e381c72601702d2f8737 100644
--- a/src/logging/log_every_timestep.c
+++ b/src/logging/log_every_timestep.c
@@ -2188,25 +2188,31 @@ void log_every_timestep(struct stardata_t * RESTRICT const stardata)
stardata->common.orbit.separation);
*/
+#define LOG_INDIVIDUAL_NOVAE
#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 %d %20.12e %20.12e %20.12e %g %g\n",
- stardata->model.time,
+ 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",
+ stardata->model.time,//1
stardata->star[0].mass,
- dt,
+ dt,//3
stardata->star[1].radius / stardata->star[1].roche_radius,
stardata->star[1].mass,//5
stardata->star[1].radius,
Mdot_net(&stardata->star[1]),
stardata->star[0].dm_novaH, //8
stardata->star[0].derivative[DERIVATIVE_STELLAR_MASS_NOVA],
- (stardata->star[0].nova || stardata->star[1].nova) ? 1 : 0, //10
- Mdot_net(&stardata->star[0]), //11 accretor rate (+ve when accreting, -ve when exploding)
- Mdot_net(&stardata->star[1]), //12 donor rate (-ve)
- stardata->common.orbit.separation, //13
- stardata->star[0].dmacc,
- stardata->star[1].dmacc
+ stardata->star[0].nova,//10
+ stardata->star[1].nova,
+ (stardata->star[0].nova || stardata->star[1].nova) ? 1 : 0, //12
+ Mdot_net(&stardata->star[0]), //13 accretor rate (+ve when accreting, -ve when exploding)
+ Mdot_net(&stardata->star[1]), //14 donor rate (-ve)
+ stardata->common.orbit.separation, //15
+ stardata->star[0].dmacc,//16
+ stardata->star[1].dmacc,//17
+ stardata->model.nova_timeout, //18
+ stardata->star[0].dtlimiter,//19
+ stardata->star[1].dtlimiter //20
);
}
#endif//LOG_INDIVIDUAL_NOVAE
@@ -2667,5 +2673,4 @@ void log_every_timestep(struct stardata_t * RESTRICT const stardata)
}
}
#endif//GAIAHRD
-
}
diff --git a/src/novae/nova_event_handler.c b/src/novae/nova_event_handler.c
index a595dfce1aef3c54d935e6a12b9cfbcf6df76cce..2ca4ad3ae911580a716bf170af0bff0688555874 100644
--- a/src/novae/nova_event_handler.c
+++ b/src/novae/nova_event_handler.c
@@ -68,7 +68,6 @@ Event_handler_function nova_event_handler(void * const eventp,
exploder->dm_novaH = 0.0;
exploder->dmacc = 0.0;
exploder->mass -= dm_lost;
- exploder->novastate = NOVA_STATE_POST;
#ifdef NUCSYN
/*
@@ -109,10 +108,24 @@ Event_handler_function nova_event_handler(void * const eventp,
donor->dmacc += dm_gain;
#ifdef NUCSYN
+ /*
+ * Put accreted material into the star's accretion layer
+ */
nucsyn_dilute_shell(
- donor->dmacc,donor->Xacc,
- dm_gain,Xnova
+ donor->dmacc,
+ donor->Xacc,
+ dm_gain,
+ Xnova
);
+
+ /*
+ * Perhaps mix the star's accretion layer into its envelope
+ */
+ if(nucsyn_thermohaline_unstable(stardata,donor))
+ {
+ nucsyn_thermohaline_mix_star(stardata,
+ donor);
+ }
/*
* "negative-yield" accreted material
diff --git a/src/nucsyn/nucsyn_prototypes.h b/src/nucsyn/nucsyn_prototypes.h
index a1a9f82a44916e8a5632b03091dc05167c20081f..9b4cab799f03a4c60ec5959bd226525220467d3d 100644
--- a/src/nucsyn/nucsyn_prototypes.h
+++ b/src/nucsyn/nucsyn_prototypes.h
@@ -268,11 +268,11 @@ double Pure_function Nonnull_some_arguments(1,2) XXsum(const struct stardata_t *
#ifdef NUCSYN_THERMOHALINE_MIXING_TIME
Constant_function double nucsyn_thermohaline_mixing_time(const double mu0,
- const double mu1,
- const double dmacc,
- const double m,
- const double l,
- const double r);
+ const double mu1,
+ const double dmacc,
+ const double m,
+ const double l,
+ const double r);
#endif
void nucsyn_set_nova_abunds(struct stardata_t * const stardata,
@@ -475,11 +475,10 @@ double Constant_function nucsyn_CN_timescale(const double t9,
void nucsyn_thermohaline_mix(struct stardata_t * const stardata,
const Star_number k, // accreting star number
- const Abundance * const Xacc, // accretion abundances
+ Abundance * const Xacc, // accretion abundances
const double dmacc, // mass accreted
const Abundance * const Xenv, // envelope abundances: pre mix
- Abundance * const nXenv, // envelope abundances: post mix (written)
- const double mforce // force mixing down to this co-ordinate
+ Abundance * const nXenv // envelope abundances: post mix (written)
);
void nucsyn_set_sigmav_from_table(struct stardata_t * const stardata,
@@ -634,6 +633,13 @@ void nucsyn_remove_mass_from_surface(
double * RESTRICT ndmacc,
Boolean *RESTRICT allocated);
+
+
+Boolean nucsyn_thermohaline_unstable(struct stardata_t * const stardata,
+ struct star_t * const star);
+
+void nucsyn_thermohaline_mix_star(struct stardata_t * const stardata,
+ struct star_t * const star);
#endif /* NUCSYN */
#endif /* NUCSYN_PROTOYPES */
diff --git a/src/nucsyn/nucsyn_thermohaline_mix.c b/src/nucsyn/nucsyn_thermohaline_mix.c
index 759ae606a4e7d82322b1ae1383b9d8a4946c4c0f..326fb8f6ee2f7d1b3d5df95b8d6e34df686f7f1f 100644
--- a/src/nucsyn/nucsyn_thermohaline_mix.c
+++ b/src/nucsyn/nucsyn_thermohaline_mix.c
@@ -2,13 +2,12 @@
#include "nucsyn.h"
#ifdef NUCSYN
-void nucsyn_thermohaline_mix(struct stardata_t * stardata,
+void nucsyn_thermohaline_mix(struct stardata_t * const stardata,
const Star_number k, // accreting star number
- const Abundance * const Xacc, // accretion abundances
+ Abundance * const Xacc, // accretion abundances
const double dmacc, // mass accreted
const Abundance * const Xenv, // envelope abundances: pre mix
- Abundance * const nXenv, // envelope abundances: post mix (written)
- const double mforce // force mixing down to this co-ordinate
+ Abundance * const nXenv // envelope abundances: post mix (written)
)
{
SETstar(k); // the accreting star
@@ -27,21 +26,36 @@ void nucsyn_thermohaline_mix(struct stardata_t * stardata,
#endif//MIXDEBUG
/* Accretion layer sinks */
- double menv = nucsyn_envelope_mass(&stardata->star[k]);
+ const double menv = nucsyn_envelope_mass(&stardata->star[k]);
Dprint("menv = %g, dmacc = %g\n",menv,dmacc);
/* Fully mix accretion layer with envelope */
- nucsyn_dilute_shell_to(
- /* Accretion layer... */
- dmacc,
- Xacc,
- /* envelope ... */
- menv,
- Xenv,
- /* New envelope */
- nXenv);
+ if(likely(Xenv != nXenv))
+ {
+ nucsyn_dilute_shell_to(
+ /* Accretion layer... */
+ dmacc,
+ Xacc,
+ /* envelope ... */
+ menv,
+ Xenv,
+ /* New envelope */
+ nXenv);
+ }
+ else
+ {
+ nucsyn_dilute_shell(
+ /* Accretion layer... */
+ dmacc,
+ Xacc,
+ /* envelope ... */
+ menv,
+ Xenv
+ );
+ }
+
if(ON_MAIN_SEQUENCE(star->stellar_type))
{
/*
diff --git a/src/nucsyn/nucsyn_thermohaline_mix_star.c b/src/nucsyn/nucsyn_thermohaline_mix_star.c
new file mode 100644
index 0000000000000000000000000000000000000000..42895d6781deb4043531987cbbedb2195341d77b
--- /dev/null
+++ b/src/nucsyn/nucsyn_thermohaline_mix_star.c
@@ -0,0 +1,21 @@
+#include "../binary_c.h"
+
+#ifdef NUCSYN
+
+/*
+ * Wrapper for nucsyn_thermohaline_mix to mix a star
+ * given in the star struct
+ */
+
+void nucsyn_thermohaline_mix_star(struct stardata_t * const stardata,
+ struct star_t * const star)
+{
+ nucsyn_thermohaline_mix(stardata,
+ star->starnum,
+ star->Xacc,
+ star->dmacc,
+ star->Xenv,
+ star->Xenv);
+}
+
+#endif
diff --git a/src/nucsyn/nucsyn_thermohaline_unstable.c b/src/nucsyn/nucsyn_thermohaline_unstable.c
new file mode 100644
index 0000000000000000000000000000000000000000..ff133dca522c7509db001654aa61b598d5140f13
--- /dev/null
+++ b/src/nucsyn/nucsyn_thermohaline_unstable.c
@@ -0,0 +1,31 @@
+#include "../binary_c.h"
+
+#ifdef NUCSYN
+
+/*
+ * Function to return TRUE if the star's current accretion layer
+ * has a higher molecular weight than the material below it,
+ * and hence is thermohaline unstable.
+ *
+ * Note: we assume an ionisation fraction of 100%.
+ */
+
+Boolean nucsyn_thermohaline_unstable(struct stardata_t * const stardata,
+ struct star_t * const star)
+{
+
+ const double mu_acc = nucsyn_molecular_weight(star->Xacc,
+ stardata,
+ 1.0);
+ const double mu_env = nucsyn_molecular_weight(star->Xenv,
+ stardata,
+ 1.0);
+
+ return mu_acc > mu_env ? TRUE : FALSE;
+}
+
+
+#endif
+
+
+
diff --git a/src/nucsyn/nucsyn_update_abundances.c b/src/nucsyn/nucsyn_update_abundances.c
index 717a590adbaedf5c9313119408fcf518edd3db4f..e41c273f5915e42e974f432507f0a1cf1eda956f 100644
--- a/src/nucsyn/nucsyn_update_abundances.c
+++ b/src/nucsyn/nucsyn_update_abundances.c
@@ -1,8 +1,10 @@
#include "../binary_c.h"
/*
- * Function to update the abundance of each star due to mass
- * mass transfer, gain and loss
+ * Subroutine to update the abundance of each star due to mass
+ * mass transfer, gain and loss.
+ *
+ * This subroutine
*/
#ifdef NUCSYN
@@ -500,8 +502,8 @@ void Hot_function nucsyn_update_abundances(struct stardata_t * RESTRICT const st
nXacc[k],
ndmacc[k],
Xenv[k],
- nXenv[k],
- -1.0);
+ nXenv[k]);
+
ndmacc[k] = 0.0;
MDEBUG("\nafter mix C12 : acc %g env %g ",
diff --git a/src/timestep/timestep_limits.c b/src/timestep/timestep_limits.c
index 26c9f38e61c9ba49d5183377794df6d64c523780..1848c2537aa8f517e10d0f4f084d1d7c9eb7f171 100644
--- a/src/timestep/timestep_limits.c
+++ b/src/timestep/timestep_limits.c
@@ -257,126 +257,40 @@ void timestep_limits(Timestep_prototype_args)
}
- if(stardata->preferences->individual_novae == TRUE)
+ if(stardata->preferences->individual_novae == TRUE &&
+ WHITE_DWARF(star->stellar_type) &&
+ (star->nova != NOVA_STATE_NONE ||
+ stardata->model.time < stardata->model.nova_timeout))
{
- if(WHITE_DWARF(star->stellar_type))
- {
- const double accretion_rate = Mdot_gain(star);
- if(Is_zero(star->time_next_nova))
- {
- /* first time setup */
- star->time_next_nova = stardata->model.time;
- star->time_prev_nova = stardata->model.time;
- }
- if(star->nova == TRUE)
- {
- /*
- * When novae are switched on and off because the system
- * becomes occasionally detached it's hard to know
- * dtnova (because mdot is sometimes zero).
- *
- * If dtnova is zero, this is the first time, then assume
- * a mass transfer rate of 1e-9 Msun/yr if nothing else.
- */
-
- double dtlim;
-
- if(Is_zero(star->dtnova))
- {
- star->dtnova = nova_recurrence_time(1e-9,
- star->mass);
- }
-
-#ifdef __DEPRECATED
- if(star->novastate == NOVA_STATE_TRIGGER)
- {
- /*
- * Nova triggered: short timestep
- * required
- */
- dtlim = stardata->preferences->minimum_timestep;
- star->dtnova = nova_recurrence_time(accretion_rate,
- star->mass);
-
- /*
- * Keep the nova state either when accreting or
- * until after nova_timeout: this keeps the timestep
- * short enough to resolve novae for a while, just in case
- * they restart.
- */
- star->nova_timeout = stardata->model.time + 3.0 * star->dtnova;
- }
- else if(star->novastate == NOVA_STATE_POST)
- {
- /*
- * Just after nova: also a short timsetep
- */
- dtlim = stardata->preferences->minimum_timestep;
- }
- else
- {
- /*
- * We have between nova explosions, try
- * to resolve them
- */
- dtlim = stardata->preferences->timestep_multipliers[DT_LIMIT_NOVAE] * star->dtnova;
- }
-#endif
-
- star->dtnova = nova_recurrence_time(accretion_rate,
+ /*
+ * When novae are switched on and off because the system
+ * becomes occasionally detached it's hard to know
+ * dtnova (because mdot is sometimes zero).
+ *
+ * If dtnova is zero, this is the first time, then assume
+ * a mass transfer rate of 1e-9 Msun/yr if nothing else.
+ */
+ const double accretion_rate = Mdot_gain(star);
+ const Time dtnova = 1e-6 * nova_recurrence_time(accretion_rate,
star->mass);
- dtlim = stardata->preferences->timestep_multipliers[DT_LIMIT_NOVAE] * star->dtnova;
-
- dtlim *= 1e-6; // convert to Myr
-
- /*
- * limit timestep because of novae
- */
- Limit_timestep(*dt,
- dtlim,
- star,
- DT_LIMIT_NOVAE);
- }
- else
- {
- /* dtnova is in YEARS, */
- double dtnova = nova_recurrence_time(accretion_rate,
- star->mass);
-
- if(0)
- printf("TIME RES NOVAE : dtnova = %g from M=%g, accretion rate %g (net %g)\n",
- dtnova,
- star->mass,
- accretion_rate,
- Mdot_net(star)
- );
- if(More_or_equal(stardata->model.time,
- star->time_next_nova))
- {
- /*
- printf("NOVANOVANOVA at %g (prev at %g, next at %g, dtnova(M=%g,Mdot=%g)=%g Myr)\n",
- stardata->model.time,
- star->time_prev_nova,
- star->time_prev_nova + dtnova*1e-6,
- star->mass,
- accretion_rate,
- dtnova*1e-6);
- */
-
- star->time_prev_nova = stardata->model.time;
- }
-
- /* limit timestep to dtnova * 0.1 */
- Limit_timestep(*dt,0.1*dtnova*1e-6,star,DT_LIMIT_NOVAE);
- }
+ /*
+ * Hence the timestep limit in Myr
+ */
+ const Time dtlim =
+ stardata->preferences->timestep_multipliers[DT_LIMIT_NOVAE] *
+ dtnova;
+ /*
+ * limit timestep because of novae
+ */
+ Limit_timestep(*dt,
+ dtlim,
+ star,
+ DT_LIMIT_NOVAE);
- /* update next nova time */
- star->time_next_nova = star->time_prev_nova + star->dtnova*1e-6;
- }
}
if(star->starnum == 0)
diff --git a/tbse b/tbse
index 2ef5d1c97aaca5edf29ec5efb715c29abe78b4e6..b9d1a5d9e950e04a6c3baa07c8b266259564c128 100755
--- a/tbse
+++ b/tbse
@@ -827,7 +827,7 @@ NOVA_RETENTION_METHOD=NOVA_RETENTION_ALGORITHM_CONSTANT
NOVA_RETENTION_FRACTION=1e-3
# Explode individual novae?
-INDIVIDUAL_NOVAE=True
+INDIVIDUAL_NOVAE=False
# Reverse accretion beta (Shara+ 1986) :
# accrete beta * the ejected mass back onto the companion.