From 0a4bacdc6e5803828dfb68ac9c9b4d9867726871 Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Fri, 18 Oct 2019 11:55:29 +0200
Subject: [PATCH] add event-based changes because of novae:

cleaned up the star->nova flag
added a new function to take into account angular momentum changes because of novae, this is the same function for the old and new codes
updated the difflog and the timestepping
---
 src/binary_c_parameters.h                     |  2 +-
 src/binary_c_structures.h                     |  4 +-
 ...ar_momentum_and_eccentricity_derivatives.c | 31 ++-----
 .../limit_accretion_rates.c                   | 29 ++++---
 src/evolution/evolution_difflog.c             |  6 +-
 src/logging/log_every_timestep.c              |  8 +-
 src/novae/nova_angular_momentum_changes.c     | 51 +++++++++++
 src/novae/nova_event_handler.c                | 86 +++++++++++++++----
 src/novae/novae_prototypes.h                  |  8 ++
 src/timestep/setup_fixed_timesteps.c          |  5 +-
 src/timestep/timestep_limits.c                | 42 ++++++---
 .../timestep_set_default_multipliers.c        |  2 +-
 12 files changed, 194 insertions(+), 80 deletions(-)
 create mode 100644 src/novae/nova_angular_momentum_changes.c

diff --git a/src/binary_c_parameters.h b/src/binary_c_parameters.h
index 253a2b109..ca7d451f6 100644
--- a/src/binary_c_parameters.h
+++ b/src/binary_c_parameters.h
@@ -718,7 +718,7 @@
  * NOVAE_TIMESTEP_FACTOR of the recurrance time
  * (time between nova explosions).
  */
-//#define TIME_RESOLVE_NOVAE
+#define TIME_RESOLVE_NOVAE
 //#define NOVAE_TIMESTEP_FACTOR 0.1
 
 /*
diff --git a/src/binary_c_structures.h b/src/binary_c_structures.h
index 93289a383..cf032a6f1 100644
--- a/src/binary_c_structures.h
+++ b/src/binary_c_structures.h
@@ -1368,7 +1368,7 @@ struct star_t {
 #endif// TIME_RESOLVE_NOVAE
     
 #ifdef INDIVIDUAL_NOVAE
-    double dm_novaH;
+    double dm_novaH,nova_timeout;
 #endif //INDIVIDUAL_NOVAE
     
     double nova_beta,nova_faml;
@@ -1554,7 +1554,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 */
diff --git a/src/binary_star_functions/angular_momentum_and_eccentricity_derivatives.c b/src/binary_star_functions/angular_momentum_and_eccentricity_derivatives.c
index 075770d68..908fa54af 100644
--- a/src/binary_star_functions/angular_momentum_and_eccentricity_derivatives.c
+++ b/src/binary_star_functions/angular_momentum_and_eccentricity_derivatives.c
@@ -201,32 +201,11 @@ void angular_momentum_and_eccentricity_derivatives(
             SETstars(k);
             if(Is_not_zero(star->derivative[DERIVATIVE_STELLAR_MASS_NOVA]))
             {
-                double q = companion->mass / star->mass;
-
-                /*
-                 * Assume novae lose mass spherically, 
-                 * and the companion gains mass in a Keplerian disc.
-                 */
-                double k = star->derivative[DERIVATIVE_STELLAR_MASS_NOVA] < 0.0 
-                    ? (2.0/3.0) : 1.0;
-                
-                star->derivative[DERIVATIVE_STELLAR_ANGMOM_NOVA] =
-                    star->derivative[DERIVATIVE_STELLAR_MASS_NOVA] *
-                    k * Pow2(star->radius) * star->omega; 
-                
-                stardata->model.derivative[DERIVATIVE_ORBIT_ANGMOM_NOVA] +=
-                    q / (1.0 + q) *
-                    star->derivative[DERIVATIVE_STELLAR_MASS_NOVA] / star->mass
-                    * stardata->common.orbit.angular_momentum;
-                
-                if(star->derivative[DERIVATIVE_STELLAR_MASS_NOVA] < 0.0)
-                {
-                    /* this is the noving WD, only let this contribute to FAML */
-                    stardata->model.derivative[DERIVATIVE_ORBIT_ANGMOM_NOVA] +=
-                        star->nova_faml *
-                        star->derivative[DERIVATIVE_STELLAR_MASS_NOVA] *
-                        stardata->common.orbit.angular_momentum;
-                }
+                nova_angular_momentum_changes(stardata,
+                                              star,
+                                              star->derivative[DERIVATIVE_STELLAR_MASS_NOVA],
+                                              &star->derivative[DERIVATIVE_STELLAR_ANGMOM_NOVA],
+                                              &stardata->model.derivative[DERIVATIVE_ORBIT_ANGMOM_NOVA]);
             }
         }
 
diff --git a/src/binary_star_functions/limit_accretion_rates.c b/src/binary_star_functions/limit_accretion_rates.c
index 322c2311f..a7c1261cd 100644
--- a/src/binary_star_functions/limit_accretion_rates.c
+++ b/src/binary_star_functions/limit_accretion_rates.c
@@ -20,6 +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,
@@ -50,7 +51,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
           mdot *= YEAR_LENGTH_IN_SECONDS / M_SUN;
           accretor->derivative[DERIVATIVE_STELLAR_MASS_RLOF_GAIN] += mdot;
         */
-
+        
         /*
          * Limit mass accretion rate.
          * Only applies if the accretion rate > 0
@@ -167,7 +168,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
              * (e.g. novae)
              */
             double accretion_rate = Mdot_net(accretor);
-
+                        
             if(accretion_rate > 0.0 &&
                COMPACT_OBJECT(accretor->stellar_type) &&
                !COMPACT_OBJECT(donor->stellar_type))
@@ -239,6 +240,7 @@ 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 */
@@ -280,6 +282,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
                          * Angular momentum is dealth with in 
                          * angular_momentum_and_eccentricity_derivatives()
                          */
+                        
                         const double f = nova_retention_fraction(stardata,
                                                                  accretor,
                                                                  donor,
@@ -287,6 +290,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
                                                                  steady_burn_rate);
                         
                         Dprint("nova_retention_factor f = %g\n",f);
+
                         
                         /*
                          * Accretion-induced irradiation 
@@ -364,13 +368,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
                              */
                             accretor->derivative[DERIVATIVE_STELLAR_H_LAYER_MASS] = accretion_rate;
                             
-                            printf("dm NOVA %g\n",
-                                   accretor->dm_novaH);
-
 #ifdef INDIVIDUAL_NOVAE_EVENTS
-                            printf("cf DM %g vs dMCRIT %g\n",
-                                   accretor->dm_novaH,
-                                   dMcrit);
                             /*
                              * Check if we have sufficient mass for explosion
                              */
@@ -386,7 +384,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;
-                                printf("novae event\n");
+                                accretor->nova = TRUE;
                                 if(Add_new_event(stardata,
                                                  BINARY_C_EVENT_NOVA,
                                                  &nova_event_handler,
@@ -423,17 +421,19 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
                                  * in the layer. We burn it all to maintain mass
                                  * conservation.
                                  */
+                                accretor->nova = TRUE;
                                 if(accretor->novastate == NOVA_STATE_TRIGGER)
                                 {
                                     /* 
                                      * A nova has been triggered.
                                      *
-                                     * The timestep is (at most) 1 year.
+                                     * 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] +=
@@ -487,6 +487,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
 
                                 if(accretor->novastate != NOVA_STATE_NONE)
                                 {
+                                    accretor->nova = TRUE;
                                     accretor->novastate = NOVA_STATE_QUIESCENT;
                                 }
                             }
@@ -501,6 +502,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
                              */
                             accretor->derivative[DERIVATIVE_STELLAR_MASS_NOVA] += 
                                 - (1.0-f) * accretion_rate;
+                            accretor->nova = TRUE;
                             Dprint("accretion rate %g -> dM/dt nova = %g\n",
                                    accretion_rate,
                                     accretor->derivative[DERIVATIVE_STELLAR_MASS_NOVA]);
@@ -590,6 +592,13 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
             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;
+            }
 #endif
         }
     } 
diff --git a/src/evolution/evolution_difflog.c b/src/evolution/evolution_difflog.c
index 9274b6781..a62434d91 100644
--- a/src/evolution/evolution_difflog.c
+++ b/src/evolution/evolution_difflog.c
@@ -386,7 +386,7 @@ static void set_diffstats(struct stardata_t * stardata,
         now->luminosity[k] = star->luminosity;
         now->novarate[k] = star->derivative[DERIVATIVE_STELLAR_MASS_NOVA];
 #ifdef INDIVIDUAL_NOVAE
-        now->novastate[k] = star->novastate;
+        now->novastate[k] = star->nova;
 #endif
         now->v_eq_ratio[k] = star->v_eq_ratio;
         now->SN_type[k] = star->SN_type;
@@ -480,13 +480,13 @@ static void cf_diffstats(struct stardata_t * RESTRICT const stardata,
 #else
             Boolean novae = Boolean_(now->novarate[k] < 0.0);
 #endif
-            if(novae == TRUE &&
+            if((novae == TRUE) &&
                lockflags[k][LOCK_NOVAE] == FALSE)
             {
                 lockflags[k][LOCK_NOVAE] = TRUE;
                 flags[BEG_NOVAE_LABEL] = TRUE;
             }
-            else if(novae == FALSE &&
+            else if((novae == FALSE) &&
                     lockflags[k][LOCK_NOVAE] == TRUE)
             {
                 lockflags[k][LOCK_NOVAE] = FALSE;
diff --git a/src/logging/log_every_timestep.c b/src/logging/log_every_timestep.c
index 3a3427202..e783f1f6d 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\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\n",
                stardata->model.time,
                stardata->star[0].mass,
                dt,
@@ -2200,7 +2200,11 @@ void log_every_timestep(struct stardata_t * RESTRICT const stardata)
                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].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
             );
     }
 #endif//INDIVIDUAL_NOVAE
diff --git a/src/novae/nova_angular_momentum_changes.c b/src/novae/nova_angular_momentum_changes.c
new file mode 100644
index 000000000..c669ff91e
--- /dev/null
+++ b/src/novae/nova_angular_momentum_changes.c
@@ -0,0 +1,51 @@
+#include "../binary_c.h"
+
+void nova_angular_momentum_changes(
+    struct stardata_t * const stardata,
+    struct star_t * const star,
+    const double dm,
+    double * const dJstar,
+    double * const dJorb
+    )
+{
+    /*
+     *  Angular momentum changes that happen because of novae.
+     *
+     * dm is either a mass lost (e.g. in a nova event) or the 
+     * time deriative. If it is a derivative, then similarly
+     * dJstar and dJorb, the changes of the stellar and orbital
+     * angular momenta, should point to the appropriate derivatives.
+     *
+     * Usually dm is negative, because the exploding star loses mass.
+     * 
+     * If dm is positive, this function should deal with nova re-accretion.
+     */
+
+    struct star_t * const companion =
+        Other_star_struct(star);
+
+    const double q = companion->mass / star->mass;
+
+    /*
+     * Assume novae lose mass spherically, 
+     * and the companion gains mass in a Keplerian disc.
+     */
+    const double k = dm < 0.0 ? (2.0/3.0) : 1.0;
+                
+    *dJstar =
+        dm * k * Pow2(star->radius) * star->omega; 
+                
+    *dJorb =
+        dm / star->mass *
+        q / (1.0 + q) *
+        stardata->common.orbit.angular_momentum;
+    
+    if(dm < 0.0)
+    {
+        /* this is the noving WD, only let this contribute to FAML */
+        *dJorb +=
+            star->nova_faml *
+            dm *
+            stardata->common.orbit.angular_momentum;
+    }
+}
diff --git a/src/novae/nova_event_handler.c b/src/novae/nova_event_handler.c
index e1e1c8a9a..2b266494b 100644
--- a/src/novae/nova_event_handler.c
+++ b/src/novae/nova_event_handler.c
@@ -16,50 +16,100 @@ Event_handler_function nova_event_handler(void * const eventp,
         (struct binary_c_nova_event_t *) data;
 
     struct star_t * const donor = nova_data->donor;
-    struct star_t * const accretor = nova_data->accretor;
-    const double f = nova_data->f;
+    struct star_t * const exploder = nova_data->accretor;
     
-    printf("caught nova event at %20.12g : accretor = %d with H-layer mass %g, donor = %d : f = %g\n",
+    printf("caught nova event at %20.12g : exploder = %d with H-layer mass %g, donor = %d : f = %g\n",
            stardata->model.time,
-           accretor->starnum,
-           accretor->dm_novaH,
+           exploder->starnum,
+           exploder->dm_novaH,
            donor->starnum,
-           f
+           nova_data->f
         );
-    
-    const double dm_lost = (1.0-f) * accretor->dm_novaH; 
+
+    /*
+     * Calculate mass lost: note that dm_lost is positive,
+     * while dm_loss_rate is negative.
+     */
+    const double dm_lost = (1.0 - nova_data->f) * exploder->dm_novaH; 
     const double dm_loss_rate = - dm_lost / stardata->model.dt; 
-    accretor->derivative[DERIVATIVE_STELLAR_MASS_NOVA] +=
-        dm_loss_rate;
 
+    /*
+     * Simulate a nova in the derivative 
+     * (although this is only really useful for logging) 
+     */
+    exploder->derivative[DERIVATIVE_STELLAR_MASS_NOVA] = dm_loss_rate;
+    
     /*
      * No hydrogen is left on the surface of the white dwarf
      */
-    accretor->dm_novaH = 0.0;
-    accretor->mass -= dm_lost;
-    accretor->novastate = NOVA_STATE_POST;
-
+    exploder->dm_novaH = 0.0;
+    exploder->mass -= dm_lost;
+    exploder->novastate = NOVA_STATE_POST;
+    
     /*
      * yield nova material
      */
     double Xnova[ISOTOPE_ARRAY_SIZE];
     nucsyn_set_nova_abunds(stardata,
-                           accretor,
+                           exploder,
                            donor->Xenv,//THIS IS WRONG
                            Xnova);
     nucsyn_calc_yields(stardata,
-                       accretor,
+                       exploder,
                        dm_lost,
                        Xnova,
                        0.0,
                        NULL,
-                       accretor->starnum,
+                       exploder->starnum,
                        YIELD_NOT_FINAL,
                        NUCSYN_SOURCE_NOVAE);
 
     /*
-     * Todo : re-accretion of a fraction
+     * Re-accretion of a fraction
+     */
+    const double dm_gain = exploder->nova_beta * dm_lost; 
+    donor->mass += dm_gain;
+    
+    /*
+     * update stellar and orbital angular momenta
+     * caused by mass loss from the exploder
+     */
+    double dJstar,dJorb;
+    nova_angular_momentum_changes(
+        stardata,
+        exploder,
+        -dm_lost,
+        &dJstar,
+        &dJorb
+        );
+    exploder->angular_momentum += dJstar;
+    stardata->common.orbit.angular_momentum += dJorb;
+    fprintf(stderr,"NOVA exploder %g %g\n",dJstar,dJorb);
+
+    /*
+     * update stellar and orbital angular momenta
+     * caused by mass gained onto the donor
+     */
+    nova_angular_momentum_changes(
+        stardata,
+        donor,
+        dm_gain,
+        &dJstar,
+        &dJorb
+        );
+    donor->angular_momentum += dJstar;
+    stardata->common.orbit.angular_momentum += dJorb;
+    fprintf(stderr,"NOVA gainer (reaccretion) %g %g\n",dJstar,dJorb);
+    
+    /*
+     * update orbit
      */
+    fprintf(stderr,"sep was %20.12e .. ",stardata->common.orbit.separation);
+    update_orbital_variables(stardata,
+                             &stardata->common.orbit,
+                             &stardata->star[0],
+                             &stardata->star[1]);
+    fprintf(stderr,"now %20.12e\n",stardata->common.orbit.separation);
     
     return NULL;
 }
diff --git a/src/novae/novae_prototypes.h b/src/novae/novae_prototypes.h
index 6c4142de1..d4f14ac32 100644
--- a/src/novae/novae_prototypes.h
+++ b/src/novae/novae_prototypes.h
@@ -29,5 +29,13 @@ Event_handler_function nova_erase_event_handler(void * eventp,
                                                 void * data);
 #endif // INDIVIDUAL_NOVAE
 
+void nova_angular_momentum_changes(
+    struct stardata_t * const stardata,
+    struct star_t * const exploder,
+    const double dm,
+    double * const dJstar,
+    double * const dJorb
+    );
+
 #endif // NOVAE_PROTOTYPES_H
 
diff --git a/src/timestep/setup_fixed_timesteps.c b/src/timestep/setup_fixed_timesteps.c
index 102ce2280..78147ccde 100644
--- a/src/timestep/setup_fixed_timesteps.c
+++ b/src/timestep/setup_fixed_timesteps.c
@@ -13,7 +13,7 @@ void setup_fixed_timesteps(struct stardata_t * RESTRICT const stardata)
      * use this by default, although allow it 
      * to be turned off on the command line.
      */
-    if(stardata->preferences->use_fixed_timestep[FIXED_TIMESTEP_YIELDS])
+    if(stardata->preferences->use_fixed_timestep[FIXED_TIMESTEP_YIELDS] == TRUE)
     {
         t->enabled = TRUE;
         t->begin =
@@ -47,7 +47,7 @@ void setup_fixed_timesteps(struct stardata_t * RESTRICT const stardata)
      * The test timestep outputs every 1Myr. It is disabled by default.
      */
     t = &stardata->model.fixed_timesteps[FIXED_TIMESTEP_TEST];
-    if(stardata->preferences->use_fixed_timestep[FIXED_TIMESTEP_TEST])
+    if(stardata->preferences->use_fixed_timestep[FIXED_TIMESTEP_TEST] == TRUE)
     {
         t->enabled = FALSE;
         t->begin = 0.0;
@@ -70,5 +70,4 @@ void setup_fixed_timesteps(struct stardata_t * RESTRICT const stardata)
      * for it in timestep.h and increment the NUMBER_OF_FIXED_TIMESTEPS
      * appropriately.
      */
-    
 }
diff --git a/src/timestep/timestep_limits.c b/src/timestep/timestep_limits.c
index 8606f7951..244fd698a 100644
--- a/src/timestep/timestep_limits.c
+++ b/src/timestep/timestep_limits.c
@@ -269,7 +269,8 @@ void timestep_limits(Timestep_prototype_args)
                 }
 
 #ifdef INDIVIDUAL_NOVAE
-                if(stardata->preferences->individual_novae == TRUE)
+                if(star->nova == TRUE &&
+                   stardata->preferences->individual_novae == TRUE)
                 {
                     /*
                      * When novae are switched on and off because the system
@@ -284,35 +285,49 @@ void timestep_limits(Timestep_prototype_args)
             
                     if(Is_zero(star->dtnova))
                     {
-                        star->dtnova = nova_recurrence_time(stardata,
-                                                            1e-9,
+                        star->dtnova = nova_recurrence_time(1e-9,
                                                             star->mass);
                     }
 
                     if(star->novastate == NOVA_STATE_TRIGGER)
                     {
-                        dtlim = 1.0;
                         /* 
-                         * We know we're in a state of mass transfer
-                         * so we can update the recurrance time.
+                         * Nova triggered: short timestep 
+                         * required
                          */
-                        star->dtnova = nova_recurrence_time(stardata,
-                                                            accretion_rate,
+                        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)
                     {
-                        dtlim = 1.0;
+                        /*
+                         * Just after nova: also a short timsetep
+                         */
+                        dtlim = stardata->preferences->minimum_timestep;
                     }
                     else
                     {
-                        dtlim = NOVAE_TIMESTEP_FACTOR 
-                            * star->dtnova;
+                        /*
+                         * We have between nova explosions, try
+                         * to resolve them
+                         */
+                        dtlim = stardata->preferences->timestep_multipliers[DT_LIMIT_NOVAE] * star->dtnova;
                     }
 
                     dtlim *= 1e-6; // convert to Myr
                     
-                    /* limit timestep because of novae */
+                    /* 
+                     * limit timestep because of novae
+                     */
                     Limit_timestep(*dt,
                                    dtlim,
                                    star,
@@ -322,8 +337,7 @@ void timestep_limits(Timestep_prototype_args)
 #endif // INDIVIDUAL_NOVAE
                 {
                     /* dtnova is in YEARS,  */
-                    double dtnova = nova_recurrence_time(stardata,
-                                                         accretion_rate,
+                    double dtnova = nova_recurrence_time(accretion_rate,
                                                          star->mass);
             
                     if(0)
diff --git a/src/timestep/timestep_set_default_multipliers.c b/src/timestep/timestep_set_default_multipliers.c
index 3a179f28a..8bc5d2197 100644
--- a/src/timestep/timestep_set_default_multipliers.c
+++ b/src/timestep/timestep_set_default_multipliers.c
@@ -57,7 +57,7 @@ void timestep_set_default_multipliers(struct preferences_t * preferences)
     preferences->timestep_multipliers[DT_LIMIT_YVT] = 1.0; /* 36 */
     preferences->timestep_multipliers[DT_LIMIT_MINIMUM_TIMESTEP] = 1.0; /* 37 */
     preferences->timestep_multipliers[DT_LIMIT_MAXIMUM_TIMESTEP] = 1.0; /* 38 */
-    preferences->timestep_multipliers[DT_LIMIT_NOVAE] = 1.0; /* 39 */
+    preferences->timestep_multipliers[DT_LIMIT_NOVAE] = 1e-2; /* 39 */
     preferences->timestep_multipliers[DT_LIMIT_ARTIFICIAL_ACCRETION] = 1.0; /* 40 */
     preferences->timestep_multipliers[DT_LIMIT_SN] = 1.0; /* 41 */
     preferences->timestep_multipliers[DT_LIMIT_MASS_GAIN] = 1.0; /* 42 */
-- 
GitLab