From 1cbe8ba5c8c92f831eaca6a2a9130ecb9e8aef04 Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Fri, 18 Oct 2019 15:43:20 +0200
Subject: [PATCH] more clean up of nova yielding: the reaccretion now works,
and when a nova goes off we set the surface layer and nova layer mass to
zero, on the assumption that any material left after the nova is mixed with
white dwarf material.
one could also argue that it would be hydrogen rich, but well,
who knows.... yet.
---
src/logging/log_every_timestep.c | 6 +-
src/novae/nova_event_handler.c | 66 +++++++++--
src/nucsyn/nucsyn_add_mass_to_surface.c | 30 +++++
src/nucsyn/nucsyn_prototypes.h | 17 +++
src/nucsyn/nucsyn_remove_mass_from_surface.c | 37 ++++++
src/nucsyn/nucsyn_update_abundances.c | 115 +++++--------------
src/nucsyn/nucsyn_update_abundances.h | 15 ---
7 files changed, 174 insertions(+), 112 deletions(-)
create mode 100644 src/nucsyn/nucsyn_add_mass_to_surface.c
create mode 100644 src/nucsyn/nucsyn_remove_mass_from_surface.c
diff --git a/src/logging/log_every_timestep.c b/src/logging/log_every_timestep.c
index e783f1f6d..b91c49d03 100644
--- a/src/logging/log_every_timestep.c
+++ b/src/logging/log_every_timestep.c
@@ -2191,7 +2191,7 @@ void log_every_timestep(struct stardata_t * RESTRICT const stardata)
#ifdef 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\n",
+ 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,
stardata->star[0].mass,
dt,
@@ -2204,7 +2204,9 @@ void log_every_timestep(struct stardata_t * RESTRICT const stardata)
(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->common.orbit.separation, //13
+ stardata->star[0].dmacc,
+ stardata->star[1].dmacc
);
}
#endif//INDIVIDUAL_NOVAE
diff --git a/src/novae/nova_event_handler.c b/src/novae/nova_event_handler.c
index 2b266494b..c59fda166 100644
--- a/src/novae/nova_event_handler.c
+++ b/src/novae/nova_event_handler.c
@@ -4,6 +4,11 @@
/*
* Event handling function for novae
+ *
+ * When you choose individual_novae to be on, nova
+ * derivatives (in star->derivative and stardata->model.derivative)
+ * are zero. Instead, novae are treated as instanteous events
+ * here, much like a supernova would be.
*/
Event_handler_function nova_event_handler(void * const eventp,
@@ -33,27 +38,50 @@ Event_handler_function nova_event_handler(void * const eventp,
const double dm_lost = (1.0 - nova_data->f) * exploder->dm_novaH;
const double dm_loss_rate = - dm_lost / stardata->model.dt;
+ /*
+ * Re-accretion of a fraction
+ */
+ const double dm_gain = exploder->nova_beta * dm_lost;
+ const double dm_gain_rate = +dm_gain / stardata->model.dt;
+
+
/*
* Simulate a nova in the derivative
* (although this is only really useful for logging)
*/
exploder->derivative[DERIVATIVE_STELLAR_MASS_NOVA] = dm_loss_rate;
+ exploder->derivative[DERIVATIVE_STELLAR_MASS_NOVA] = dm_gain_rate;
/*
- * No hydrogen is left on the surface of the white dwarf
+ * Remove material from the white dwarf:
+ * no hydrogen is left on the surface
+ *
+ * Here I just remove dm_lost from dm_novaH... but you
+ * could argue this material is mixed in and becomes part
+ * of the white dwarf after the nova.
*/
exploder->dm_novaH = 0.0;
+ exploder->dmacc = 0.0;
exploder->mass -= dm_lost;
exploder->novastate = NOVA_STATE_POST;
+
+ /*
+ * Accrete some on the companion
+ */
+ donor->mass += dm_gain;
/*
- * yield nova material
+ * calculate nova yield abundances
*/
double Xnova[ISOTOPE_ARRAY_SIZE];
nucsyn_set_nova_abunds(stardata,
exploder,
donor->Xenv,//THIS IS WRONG
Xnova);
+
+ /*
+ * yield material in nova explosion
+ */
nucsyn_calc_yields(stardata,
exploder,
dm_lost,
@@ -64,11 +92,35 @@ Event_handler_function nova_event_handler(void * const eventp,
YIELD_NOT_FINAL,
NUCSYN_SOURCE_NOVAE);
- /*
- * Re-accretion of a fraction
- */
- const double dm_gain = exploder->nova_beta * dm_lost;
- donor->mass += dm_gain;
+ /* should we call nucsyn remove dm here? */
+
+ if(Is_not_zero(dm_gain))
+ {
+
+ /*
+ * Put material on companion star: we put it
+ * in the accretion layer, let it get mixed
+ * on the next timestep
+ */
+ nucsyn_dilute_shell(
+ donor->dmacc,donor->Xacc,
+ dm_gain,Xnova
+ );
+ donor->dmacc += dm_gain;
+
+ /*
+ * "negative-yield" accreted material
+ */
+ nucsyn_calc_yields(stardata,
+ donor,
+ 0.0,
+ NULL,
+ dm_gain,
+ Xnova,
+ exploder->starnum,
+ YIELD_NOT_FINAL,
+ NUCSYN_SOURCE_NOVAE);
+ }
/*
* update stellar and orbital angular momenta
diff --git a/src/nucsyn/nucsyn_add_mass_to_surface.c b/src/nucsyn/nucsyn_add_mass_to_surface.c
new file mode 100644
index 000000000..546450f46
--- /dev/null
+++ b/src/nucsyn/nucsyn_add_mass_to_surface.c
@@ -0,0 +1,30 @@
+#include "../binary_c.h"
+
+#ifdef NUCSYN
+void nucsyn_add_mass_to_surface(const Star_number k,
+ const double dm,
+ const Abundance * RESTRICT const X,
+ double * RESTRICT const ndmacc,
+ const Abundance * const * const RESTRICT Xacc,
+ Abundance ** RESTRICT newXacc)
+{
+ /*
+ * Add material to the surface of a star by putting it
+ * in its accretion layer
+ */
+ if(ndmacc[k]>VERY_TINY)
+ {
+ /* Existing accretion layer : have to mix */
+ nucsyn_dilute_shell(
+ ndmacc[k],newXacc[k],
+ dm,X
+ );
+ }
+ else
+ {
+ /* new accretion layer */
+ Copy_abundances(X,newXacc[k]);
+ }
+ ndmacc[k] += dm;
+}
+#endif//NUCSYN
diff --git a/src/nucsyn/nucsyn_prototypes.h b/src/nucsyn/nucsyn_prototypes.h
index a1c6b888f..a1a9f82a4 100644
--- a/src/nucsyn/nucsyn_prototypes.h
+++ b/src/nucsyn/nucsyn_prototypes.h
@@ -617,6 +617,23 @@ Abundance nucsyn_elemental_abundance_by_number(const char * RESTRICT const eleme
void nucsyn_make_icache(struct store_t * const store);
void nucsyn_free_icache(Isotope ** icache);
+void nucsyn_add_mass_to_surface(const Star_number k,
+ const double dm,
+ const Abundance * RESTRICT const X,
+ double * RESTRICT const ndmacc,
+ const Abundance * const * RESTRICT const Xacc,
+ Abundance ** RESTRICT nXacc);
+
+void nucsyn_remove_mass_from_surface(
+#if defined DEBUG && DEBUG > 0
+ struct stardata_t * stardata,
+#endif
+ struct star_t * RESTRICT star,
+ const double dm,
+ Abundance ** RESTRICT X,
+ double * RESTRICT ndmacc,
+ Boolean *RESTRICT allocated);
+
#endif /* NUCSYN */
#endif /* NUCSYN_PROTOYPES */
diff --git a/src/nucsyn/nucsyn_remove_mass_from_surface.c b/src/nucsyn/nucsyn_remove_mass_from_surface.c
new file mode 100644
index 000000000..1d61c9312
--- /dev/null
+++ b/src/nucsyn/nucsyn_remove_mass_from_surface.c
@@ -0,0 +1,37 @@
+#include "../binary_c.h"
+
+void nucsyn_remove_mass_from_surface(
+#if defined DEBUG && DEBUG > 0
+ struct stardata_t * stardata,
+#endif
+ struct star_t * RESTRICT star,
+ const double dm,
+ Abundance ** RESTRICT Xp,
+ double * RESTRICT ndmacc,
+ Boolean * RESTRICT allocated)
+{
+ /*
+ * Remove material dm (>0.0) from the star, and
+ * set X, the abundance of the material removed,
+ * if X is non-NULL.
+ *
+ * Allocated can be NULL but only if X is NULL.
+ */
+
+ ndmacc[star->starnum] = Max(0.0, ndmacc[star->starnum] - dm);
+ if(likely(Xp!=NULL))
+ {
+ nucsyn_remove_dm_from_surface(
+#if defined DEBUG && DEBUG > 0
+ stardata,
+#endif
+ star,
+ dm,
+ Xp,
+ allocated);
+ }
+ else if(allocated != NULL)
+ {
+ *allocated = FALSE;
+ }
+}
diff --git a/src/nucsyn/nucsyn_update_abundances.c b/src/nucsyn/nucsyn_update_abundances.c
index a85968635..717a590ad 100644
--- a/src/nucsyn/nucsyn_update_abundances.c
+++ b/src/nucsyn/nucsyn_update_abundances.c
@@ -227,12 +227,12 @@ void Hot_function nucsyn_update_abundances(struct stardata_t * RESTRICT const st
/*
* Wind accretion, abundance calculated above
*/
- add_material(k,
- dm_wind_gain,
- Xwind_gain,
- ndmacc,
- Xacc,
- nXacc);
+ nucsyn_add_mass_to_surface(k,
+ dm_wind_gain,
+ Xwind_gain,
+ ndmacc,
+ Xacc,
+ nXacc);
/*
* RLOF accretion and non-conservative loss:
@@ -249,12 +249,12 @@ void Hot_function nucsyn_update_abundances(struct stardata_t * RESTRICT const st
&allocated_XRLOF_gain);
/* only the RLOF gain is accreted */
- add_material(k,
- dm_RLOF_gain,
- XRLOF_gain,
- ndmacc,
- Xacc,
- nXacc);
+ nucsyn_add_mass_to_surface(k,
+ dm_RLOF_gain,
+ XRLOF_gain,
+ ndmacc,
+ Xacc,
+ nXacc);
#ifdef DISCS
{
@@ -265,12 +265,12 @@ void Hot_function nucsyn_update_abundances(struct stardata_t * RESTRICT const st
if(stardata->common.ndiscs>=1 &&
stardata->common.discs[0].type==DISC_CIRCUMBINARY)
{
- add_material(k,
- dm_cbdisc_gain,
- stardata->common.discs[0].X,
- ndmacc,
- Xacc,
- nXacc);
+ nucsyn_add_mass_to_surface(k,
+ dm_cbdisc_gain,
+ stardata->common.discs[0].X,
+ ndmacc,
+ Xacc,
+ nXacc);
}
}
#endif // DISCS
@@ -278,7 +278,7 @@ void Hot_function nucsyn_update_abundances(struct stardata_t * RESTRICT const st
/*
* Wind loss
*/
- remove_material(
+ nucsyn_remove_mass_from_surface(
#if defined DEBUG && DEBUG > 0
stardata,
#endif
@@ -292,7 +292,7 @@ void Hot_function nucsyn_update_abundances(struct stardata_t * RESTRICT const st
/*
* RLOF and excretion disc loss
*/
- remove_material(
+ nucsyn_remove_mass_from_surface(
#if defined DEBUG && DEBUG > 0
stardata,
#endif
@@ -357,7 +357,7 @@ void Hot_function nucsyn_update_abundances(struct stardata_t * RESTRICT const st
/*
* Nova loss from the white dwarf
*/
- remove_material(
+ nucsyn_remove_mass_from_surface(
#if defined DEBUG && DEBUG > 0
stardata,
#endif
@@ -376,12 +376,12 @@ void Hot_function nucsyn_update_abundances(struct stardata_t * RESTRICT const st
/*
* Nova "reaccretion"
*/
- add_material(k,
- dm_novae,
- Xnovae,
- ndmacc,
- Xacc,
- nXacc);
+ nucsyn_add_mass_to_surface(k,
+ dm_novae,
+ Xnovae,
+ ndmacc,
+ Xacc,
+ nXacc);
Add_yield(0.0,NULL,
dm_novae,Xnovae,
NUCSYN_SOURCE_NOVAE);
@@ -627,67 +627,6 @@ static int id_wind_source(const struct star_t * const star,
}
#endif //NUCSYN_ID_SOURCES
-static void add_material(const Star_number k,
- const double dm,
- const Abundance * RESTRICT const X,
- double * RESTRICT const ndmacc,
- const Abundance * const * const RESTRICT Xacc,
- Abundance ** RESTRICT newXacc)
-{
- /*
- * Add material to the surface of a star by putting it
- * in its accretion layer
- */
- if(ndmacc[k]>VERY_TINY)
- {
- /* Existing accretion layer : have to mix */
- nucsyn_dilute_shell(
- ndmacc[k],newXacc[k],
- dm,X
- );
- }
- else
- {
- /* new accretion layer */
- Copy_abundances(X,newXacc[k]);
- }
- ndmacc[k] += dm;
-}
-
-static void remove_material(
-#if defined DEBUG && DEBUG > 0
- struct stardata_t * stardata,
-#endif
- struct star_t * RESTRICT star,
- const double dm,
- Abundance ** RESTRICT Xp,
- double * RESTRICT ndmacc,
- Boolean * RESTRICT allocated)
-{
- /*
- * Remove material dm (>0.0) from the star, and
- * set X, the abundance of the material removed,
- * if X is non-NULL.
- *
- * Allocated can be NULL but only if X is NULL.
- */
- ndmacc[star->starnum] = Max(0.0, ndmacc[star->starnum] - dm);
- if(likely(Xp!=NULL))
- {
- nucsyn_remove_dm_from_surface(
-#if defined DEBUG && DEBUG > 0
- stardata,
-#endif
- star,
- dm,
- Xp,
- allocated);
- }
- else if(allocated != NULL)
- {
- *allocated = FALSE;
- }
-}
#endif /* NUCSYN */
diff --git a/src/nucsyn/nucsyn_update_abundances.h b/src/nucsyn/nucsyn_update_abundances.h
index 5fa8a4db5..1fe48b669 100644
--- a/src/nucsyn/nucsyn_update_abundances.h
+++ b/src/nucsyn/nucsyn_update_abundances.h
@@ -103,21 +103,6 @@ double coratio(double c,double o)
#endif
-static void add_material(const Star_number k,
- const double dm,
- const Abundance * RESTRICT const X,
- double * RESTRICT const ndmacc,
- const Abundance * const * RESTRICT const Xacc,
- Abundance ** RESTRICT nXacc);
-static void remove_material(
-#if defined DEBUG && DEBUG > 0
- struct stardata_t * stardata,
-#endif
- struct star_t * RESTRICT star,
- const double dm,
- Abundance ** RESTRICT X,
- double * RESTRICT ndmacc,
- Boolean *RESTRICT allocated);
#ifdef NUCSYN_ID_SOURCES
--
GitLab