Timesteps - Githubissues

Describe the bug The way we calculate timesteps is fairly convoluted and could be streamlined. In particular, many of the functions have similar names and they should be renamed to clarify what each one specifically does. We also currently do not use Ilya's smart step-size reduction, which is in a separate (but related) issue.

Below are all the relevant functions which are called for binary stars when dt is calculated (with some unrelated parts redacted and replaced as my own comments). I'll use these for reference below, and they are pretty scattered throughout the code base, so its good to have them all collected in one place. The commented values formatted as // [n] function as equation numbers to easily see which function calls which other function

Label the issue Please label the 'severity' and 'urgency' of this issue. You can choose:

urgency_low - This issue is not urgent

severity_moderate - This is a moderately severe bug

To Reproduce N/A

Expected behavior This should be streamlined for readability, and Ilya's smart step-sizes should be implemented.

Versioning (please complete the following information):

COMPAS [e.g. v02.11.03]

######################################################

// [0]
 EVOLUTION_STATUS BaseBinaryStar::Evolve() {   

     EVOLUTION_STATUS evolutionStatus = EVOLUTION_STATUS::CONTINUE;

     if (HasStarsTouching()) {   
         // stop evolution
     }

     PrintDetailedOutput(m_Id);                                                                                                              
     // print (log) detailed output, first time

     // evolve the current binary up to the maximum evolution time (and number of steps)
     double dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep()) / 1000.0;                                 

     while (evolutionStatus == EVOLUTION_STATUS::CONTINUE) {                                                                    
         // perform binary evolution - iterate 

         // Evolve each star individually according to SSE
         EvolveOneTimestep(dt);                                                                                                        

         // Check for errors from SSE, reset EVOLUTION_STATUS if so
         if anyErrorsFromSSE {
              evolutionStatus = EVOLUTION_STATUS::SOME_ERROR
         }
         else {                                                                                                                      
              // continue evolution

             // Evolve the full binary 
             EvaluateBinary(dt);                                                                                                                     

             PrintDetailedOutput(m_Id);                                                                                              
             // print (log) detailed output for each timestep 

             // Check for errors from BSE, or too many timesteps, reset EVOLUTION_STATUS if so
             if anyErrorsFromBSE {
                  evolutionStatus = EVOLUTION_STATUS::SOME_ERROR
             }
         }

         // Calculate timestep 
         //                          [1]                          [1]
         dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep());   // new timestep
         //          [6]              [4]
         dt = CalculateTimestep(ChooseTimestep(dt));                       // calculate next timestep
     }

     return evolutionStatus;
 }

// [1]
 double BaseStar::CalculateTimestep() { 

     // Calculate of GB Params and Timescales for the star

     //              [2]
     double dt = ChooseTimestep(m_Age);

     // if we decide to apply an additional timestep reduction factor, do so here
     if (utils::Compare(TIMESTEP_REDUCTION_FACTOR, 1.0) != 0) {    
         if (!IsOneOf({ STELLAR_TYPE::MS_LTE_07, STELLAR_TYPE::MS_GT_07 })) {   
             dt /= TIMESTEP_REDUCTION_FACTOR;         
         }
     }

     //          [3]
     return LimitTimestep(dt);
 }

Note: TIMESTEP_REDUCTION_FACTOR = 1.0, so unless we change that, the middle loop will never run. It's not clear to me why we would want to artificially reduce the timestep, and especially why we would want to only want to do so for non-MS stars.

// [2]
(BaseStar.h)
virtual double ChooseTimestep(const double p_Time)     { return m_Dt; }

Note: This is a virtual function which is overridden by all the stellar type classes. Effectively, each one calculates the SSE timestep using the Hurley 2000 values for \delta t_k and \delta t_e, where \delta t_k corresponds to a fixed number of timesteps per phase, and \delta t_e is the time to the end of the phase (so we don't overshoot).

We also limit by the NUCLEAR_MINIMUM_TIMESTEP, although this is hardcoded in COMPAS, instead of variable on the stellar mass and type (corresponding to \delta t_N in Hurley). In these methods, we do not limit by the other Hurley timestep parameters \delta t_ml, and \delta t_R.

// [3]
 double BaseStar::LimitTimestep(const double p_Dt) {

     double dt = p_Dt;

     // COEN NEIJSSEL 16/11/2016: cap timestep according to max 1% change in mass star.
     if (OPTIONS->UseMassLoss()) {

         double mDot     = CalculateMassLossRate();   
         // First, calculate mass loss rate

         double massLoss = CalculateMassLoss_Static(m_Mass, mDot, dt);     
         // Next, calculate mass loss - limited to (mass * MAXIMUM_MASS_LOSS_FRACTION)

         if (utils::Compare(massLoss, 0.0) > 0) {      // No change if no mass loss
             double dtWind = massLoss / (mDot * 1.0E6);    
             // Calculate timestep to match (possibly limited) mass loss
                    dt     = min(dt, dtWind);        // choose dt

                    m_Mdot = massLoss / (dt * 1.0E6);          
                    // Reset mass loss rate to match (possibly limited) mass loss
         }
     }

     return dt;
 }

Note: Jeff comments here that this may modify m_Mdot, which we probably do not want to do here. This section should just limit the timestep according to the given mass loss rate, which should be checked and capped wherever it is set. This is the Hurley \delta t_ml, though that should be clearer in the name.

Also, 1.0E6 is a magic number, and it's unclear (to me at least) where it comes from. Is this for unit conversion?

// [4]
 double BaseBinaryStar::ChooseTimestep(const double p_Dt) {

     double timestep = 0.0;      // ALEJANDRO - 02/12/2015 - MASSLESS_REMNANT timescale = 0.0

     if (!m_Star1->IsOneOf({ STELLAR_TYPE::MASSLESS_REMNANT }) && !m_Star2->IsOneOf({ STELLAR_TYPE::MASSLESS_REMNANT })) {

         //              [5]
         double dt1 = CalculateDt(p_Dt, m_Star1, m_Star2);
         double dt2 = CalculateDt(p_Dt, m_Star2, m_Star1);

         timestep = std::min(dt1, dt2);
     }

     return std::max(timestep, NUCLEAR_MINIMUM_TIMESTEP);
 }

Note: This function manages the possible dt limit from either star if its near RLOF (see method [5]). However, limiting to NUCLEAR_MINIMUM_TIMESTEP is redundant with the derived versions of [2], ChooseTimestep().

// [5]
 double BaseBinaryStar::CalculateDt(p_Dt, p_Primary, p_Secondary) {

     double newDt = p_Dt;                                                                                                    
     // default is unchanged

     double tmp     = RSOL_TO_AU / (m_SemiMajorAxisPrime * (1.0 - m_Eccentricity));
     double RLRatio = tmp * p_Primary->Radius() / CalculateRocheLobeRadius_Static(p_Primary->Mass(), p_Secondary->Mass());

     if (utils::Compare(RLRatio, MASS_TRANSFER_THRESHOLD) >= 0) {                                                            

         if (p_Primary->DetermineEnvelopeType() == ENVELOPE::RADIATIVE) {
             newDt = m_FastPhaseCaseA ? p_Primary->CalculateThermalTimescale() 
                                      : p_Primary->CalculateDynamicalTimescale();  
         }
         else {
             newDt = p_Primary->CalculateThermalTimescale();
         }
     }
     else if (utils::Compare(RLRatio, ROCHE_LOBE_UPPER_THRESHOLD) >= 0) {
         newDt /= 10.0;
     }
     else if (utils::Compare(RLRatio, ROCHE_LOBE_LOWER_THRESHOLD) >= 0) {
         newDt /= 5.0;
     }

     return newDt;
 }

Note: The tmp variable should be renamed, and some of the other names could be clearer, but this correctly captures the idea that as a star gets close to its RL, the timestep should be steadily reduced. Also, Jeff comments that m_FastPhaseCaseA (binary variable, not stellar) is NEVER set TRUE, but is not sure for single stars. For reference, ROCHE_LOBE_LOWER_THRESHOLD = 0.70, ROCHE_LOBE_UPPER_THRESHOLD = 0.9, and MASS_TRANSFER_THRESHOLD = 0.95, so these are the RL filling fractions that correspond to the timestep reductions above.

// [6]
 double BaseBinaryStar::CalculateTimestep(const double p_Dt) {

     BinaryConstituentStar* star1Copy = new BinaryConstituentStar(*m_Star1);     // copy star1

     //                            [7]
     double dt1 = star1Copy->EvolveOneTimestep(p_Dt);                                    
     // evolve the copy one timestep and get suggested timestep

     delete star1Copy; star1Copy = nullptr;                                      // nuke the copy

     // rinse and repeat for star2
     BinaryConstituentStar* star2Copy = new BinaryConstituentStar(*m_Star2);     // copy star2

     //                           [7]
     double dt2 = star2Copy->EvolveOneTimestep(p_Dt);                            
     // evolve the copy one timestep and get suggested timestep

     delete star2Copy; star2Copy = nullptr;                                      // nuke the copy

     return std::min(dt1, dt2);
 }

Note: It seems to me to be a waste to evolve both star copies to see if the timesteps match, only to find the correct time and re-evolve the actual stars by the correct amount. If we're evolving copies anyway, why don't we just copy the copy back to the original star after the correct dt has been found? I'm not an expert in all this though, so maybe there's a good reason not to do that.

// [7]
 double Star::EvolveOneTimestep(const double p_Dt) {    

     double       dt = p_Dt;    

     STELLAR_TYPE stellarType;    

     bool         takeTimestep = false;    
     int          retryCount   = 0;    

     while (!takeTimestep) {                                                                                     
     // do this until a suitable timestep is found

         SaveState();                                                                                            
         // save the state of the star - in case we want to revert later

         double minTimestep = std::max(m_Star->CalculateDynamicalTimescale(), ABSOLUTE_MINIMUM_TIMESTEP);        

         stellarType = AgeOneTimestep(dt, false);                                                                
         // age the star one time step - modifies attributes, but only of copy star

         // determine if the timestep should be taken    
         // don't take the timestep if we stepped too far    

         takeTimestep = true;                                                                                    

         if (utils::Compare(m_Star->CalculateRadialChange(), MAXIMUM_RADIAL_CHANGE) >= 0) { 
         // too much change?    

             if (utils::Compare(dt, minTimestep) <= 0) { 
             // yes - already at or below minimum timestep?    

                 takeTimestep = true;      
                 // yes - just take the last timestep    
             }    
             else {                                                                                              
             // not at or below dynamical - reduce timestep and tr

                 retryCount++;                                                                                   
                 // increment retry count    

                 if (retryCount > MAX_TIMESTEP_RETRIES) { 
                 // too many retries?    

                     takeTimestep = true;                                                                        
                     // yes - take the last timestep anyway    
                 }    
                 else {                                                                                          
                 // not too many retries - retry with smaller timestep

                     if (RevertState()) {                                                                        
                     // revert to last state ok?

                         dt = dt / 2.0;                                                                          
                         // yes - halve the timestep (limit to minimum)      J

                         takeTimestep = false;                                                                   
                         // previous timestep discared - use new one
                     }
                     else {                                                                                      
                     // revert failed

                         takeTimestep = true;       
                         // take the last timestep anyway
                     }
                 }
             }
         }
     }

     // take the timestep

     (void)SwitchTo(stellarType);                                                                                
     // switch phase if required 

     (void)m_Star->ResolveMassLoss();                                                                            
     // apply wind mass loss if required  

      m_Star->CalculateAllTimescales();                                                                          
      // calculate dynamical, thermal, nuclear and radial expansion timescales

     return dt;                                                                                                  
     // return the timestep actually taken
 }

Note: As far as I can tell, this function essentially just estimates the Hurley \delta t_R parameter (which caps the relative radial difference) iteratively by actually evolving the star copy (through SSE only) and retrying if the step was too big. Could we do this analytically, or does it depend in some complex way on the wind mass loss? If so, will the SSE evolution tell us that? We calculate wind mass loss in the EvolveBinary function, so this might be implemented incorrectly anyway. Why do we need to bother with setting and reverting the state? This also calculates Timescales and Mass loss, but for the copy star. This is unnecessary, and should be bypassed if possible.

###################################################

Apologies for the really long post. I think we could brainstorm a simpler solution to this somewhat convoluted network of function calls, but I need to run, so I'll leave this here for now and write up a proposed fix (in pseudo code) a little later. But if anyone has an answer to any of my questions above, please clarify them here.

Yes, timestep calculation is a bit convoluted. Best if, as you suggest, we brainstorm at the pow-wow.

Also, 1.0E6 is a magic number, and it's unclear (to me at least) where it comes from. Is this for unit conversion?

(From memory) I think all the mass loss rate calculations return mDot in Msol/yr, but dT is in Myr. I don't know why they are in different units of time, but I think I didn't change either of them because I wasn't sure at the time whether that would have side effects elsewhere in the code.

If we're evolving copies anyway, why don't we just copy the copy back to the original star after the correct dt has been found?

Yeah, I think Ilya and I discussed the same thing when I was refactoring. It's possible to do (that's one of the reasons I wrote the Star::Clone(), Star::operator = overload, Star::SaveState() and Star::RevertState() functions), but we have to be sure that when we evolve the copy to find the right dT we're actually evolving all the attributes that need to be evolved - and I think it got messy for binaries (but don't remember the details - I just remember that the binary code has hooks into the timestep calculations). I deferred that because we were going to implement the adaptive timestep. I also didn't really want to change the timestep calculation until we were sure that the refactored code reproduced the legacy results - but now's a good time to look at implementing the adaptive timestep I think.

Why do we need to bother with setting and reverting the state?

Because that's how we find the right timestep - in this implementation. If we stepped too far we go back one step and shorten the timestep - and keep doing that until we find the best timestep. I suspect I'm not understanding your question.

This also calculates Timescales and Mass loss, but for the copy star. This is unnecessary, and should be bypassed if possible.

I don't remember well enough, so this could be wrong (I should go back and re-read the code), but isn't it the copy star we are evolving through multiple timesteps in order to find the best timestep (i.e. the biggest timestep that doesn't overstep something interesting)? If we don't calculate timescales and mass loss, won't that mean it doesn't evolve as it should, so we couldn't trust the timestep calculation? Also, if we don't evolve all the attributes properly, then we couldn't just copy the evolved values to the star as you suggest above, could we?

Yes, timestep calculation is a bit convoluted. Best if, as you suggest, we brainstorm at the pow-wow.

Also, 1.0E6 is a magic number, and it's unclear (to me at least) where it comes from. Is this for unit conversion?

(From memory) I think all the mass loss rate calculations return mDot in Msol/yr, but dT is in Myr. I don't know why they are in different units of time, but I think I didn't change either of them because I wasn't sure at the time whether that would have side effects elsewhere in the code.

If we're evolving copies anyway, why don't we just copy the copy back to the original star after the correct dt has been found?

Yeah, I think Ilya and I discussed the same thing when I was refactoring. It's possible to do (that's one of the reasons I wrote the Star::Clone(), Star::operator = overload, Star::SaveState() and Star::RevertState() functions), but we have to be sure that when we evolve the copy to find the right dT we're actually evolving all the attributes that need to be evolved - and I think it got messy for binaries (but don't remember the details - I just remember that the binary code has hooks into the timestep calculations). I deferred that because we were going to implement the adaptive timestep. I also didn't really want to change the timestep calculation until we were sure that the refactored code reproduced the legacy results - but now's a good time to look at implementing the adaptive timestep I think.

^Thanks for clarifying all these.

Why do we need to bother with setting and reverting the state?

Because that's how we find the right timestep - in this implementation. If we stepped too far we go back one step and shorten the timestep - and keep doing that until we find the best timestep. I suspect I'm not understanding your question.

This was a continuation of my thoughts from the previous sentence, which you answered in the block above.

This also calculates Timescales and Mass loss, but for the copy star. This is unnecessary, and should be bypassed if possible.

I don't remember well enough, so this could be wrong (I should go back and re-read the code), but isn't it the copy star we are evolving through multiple timesteps in order to find the best timestep (i.e. the biggest timestep that doesn't overstep something interesting)? If we don't calculate timescales and mass loss, won't that mean it doesn't evolve as it should, so we couldn't trust the timestep calculation? Also, if we don't evolve all the attributes properly, then we couldn't just copy the evolved values to the star as you suggest above, could we?

In this function, dt is calculated and is the thing that is returned. Since we're working with a copy star, all the other functions that change the attributes of the copy (but leave dt unaffected) don't actually change anything. I think this is maybe because we use this same function elsewhere in the code, where changing those attributes is important, but I could be wrong.

In this function, dt is calculated and is the thing that is returned. Since we're working with a copy star, all the other functions that change the attributes of the copy (but leave dt unaffected) don't actually change anything. I think this is maybe because we use this same function elsewhere in the code, where changing those attributes is important, but I could be wrong.

Oh, yeah, sorry. I thought what you were talking about was inside the loop. Inside the loop we do need to change the attributes, even of the copy - because we need to know if we overstepped. But yes, if we knew that we didn't want the changes outside the loop applied we wouldn't need to do them - but as the code stands at the moment we don't know that. The only reason not to execute those calls is for performance - I don't know that we do it often enough for it to matter in the overall scheme, but yes we could pass a parameter to the function to indicate that all we are interested in is the return value.

I think this is maybe because we use this same function elsewhere in the code, where changing those attributes is important, but I could be wrong.

You're talking about Star::EvolveOneTimestep(), right? Yeah, it's the main function that gets called to advance the evolution of the star being evolved by a single timestep - it's called for all stars, not just copies. And it's called to evolve copies of stars where we want to know what the (copy of the) star looks like after evolving one timestep - so in those cases we want to calculate the timescales and apply the mass loss (we're not just looking for dT).

But we don't actually do anything with the copy star, right? We calculate the timescales and mass loss, but in function [6], the copy is deleted once the dt has been retrieved.

There are other places we evolve a copy star. (sorry for brevity - I'm out, so can't look at the code, so can't point out where)

On Tue, 7 Jul 2020 20:35 Reinhold Willcox, notifications@github.com wrote:

But we don't actually do anything with the copy star, right? We calculate the timescales and mass loss, but in function [6], the copy is deleted once the dt has been retrieved.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/TeamCOMPAS/COMPAS/issues/290#issuecomment-654762091, or unsubscribe https://github.com/notifications/unsubscribe-auth/AGM5XNSD2Q5HTR6HSSXITL3R2L26BANCNFSM4ORJCXQQ .

Look at:

BaseStar::CalculateZetaNuclear() BaseStar::CalculateZetaThermal() BaseBinaryStar::CalculateAdaptiveRocheLobeOverFlow() BaseBinaryStar::CalculateMassTransferFastPhaseCaseA() BaseBinaryStar::CalculateMassTransfer() (may not actually evolve here)

Summary of current plan

My understanding of timesteps is that we have the following constraints.

From Hurley, for SSE only:

dt_k - a fixed fraction of either 1) the total time a star will spend in a given phase or 2) the remaining time on that phase (i.e the size of the timesteps will drop exponentially)
- dt_e - the time remaining to the end of the current phase
- dt_N - the remaining nuclear lifetime
- dt_ml - a timescale which limits the relative mass loss (presumably due to winds)
- dt_R - a timescale which limits the relative radial change

dt_R also implicitly includes an absolute minimum timescale of either the dynamical timescale of the star, or ABSOLUTE_MINIMUM_TIMESTEP, call these:

dt_dyn
dt_MIN

We discussed adding in another SSE parameter for dt_mc for core mass. This will be a bit trickier to handle, because when a star changes type, we might expect the core mass to balloon from 0 to some large-ish number. But that should be the only caveat.

dt_mc

Since they apply to SSE only, we can calculate them separately for each star to get dt_star = max( min(dt_k, dt_e, dt_N, dt_ml, dt_R, dt_mc), dt_dyn, dt_MIN) and then dt_SSE = min(dt_star1, dt_star2). This is just the Hurley definition, with the additional lower limits from dynamical and total minimum timescales. Hurley comments that some of the timescale choices are for aesthetic plotting purposes, and could be increased to improve speed. I think this pertains exclusively to the dt_k parameter, specifically the size of the fixed fraction of the phase time.

In addition, we have the following BSE constraint: If either star has R/RL close to 1 (i.e 0.7, 0.9), we reduce the timestep to better resolve the critical moments leading up to RLOF. If it is overflowing (R/RL >= 0.95), we fix the timestep to either the dynamical or thermal timescale of the donor, as appropriate.

Note 1: We use a fixed Nuclear timescale, which is not star dependent, and which is not the same as the remaining Nuclear lifetime. This may or may not be important - I'm also not sure if it's fundamentally different to the remaining lifetime on the phase, dt_e.

Note 2: The dt_ml is a little funky. It calculates the mass loss, possibly resetting the rate of mass loss if it's too large. In my opinion, the rate of mass loss is physical and should not change (after it's been calculated), and we should instead shorten the timestep as necessary in order to keep the overall relative mass difference below our threshold.

Note 3: It's not clear to me what should happen if, e.g, the dynamical timescale of the star is longer than the mass or radial relative difference timescales - do we stick to the shortest timescale to avoid big relative jumps, or do we go with the longer timescale to avoid over-resolving?

======================================

Proposal

We could simplify a lot of this significantly. The dt_k, dt_e, and dt_N are all straightforward to calculate analytically, but they together depend on the stellar type and age (and are thus deflected to derived classes). For the relative change timescales (dt_R, dt_ml, dt_mc), we currently evolve a copy star repeatedly until the relative change is low enough. But I think we could make this easier if we only evolve the copy star once, and use the endstate to estimate Mdot, Rdot, and MCdot. These, together with the maximum relative difference fractions give us the 3 timescales.

Here is a brief rundown of the code for what I'm suggesting, in words first, and (pseudo-)code after:

For each star, calculate the dt_k, dt_e, dt_N based on its stellar type.
From the minimum of these, get a proposed dt and evolve a copy of the star for this timestep.
Get the average change in parameters (radius, mass, and core mass) divided by the proposed timestep as a proxy for the time derivative, and then calculate the maximal timestep which keeps the relative changes in these parameters below our thresholds. (This saves having to re-evolve a copy star several times).
This value is the final timestep for single stars, but for binaries, we take the minimum of this timestep for the two stars.
If we are in a binary, we may further scale down this timescale depending on how close we are to RLOF, so we check the R*/RL ratio for both and rescale accordingly.

(In BaseBinaryStar::Evolve() )

// Get the minimum timestep that satisfies both stars individually
dtSSE = std::min( m_Star1->CalculateTimestepFromSSE(), 
                  m_Star2->CalculateTimestepFromSSE());

// Reduce the timestep as necessary if either star is close to RLOF
dt = std::min( ReduceTimestepIfCloseToRLOF(dtSSE, m_Star1, m_Star2),
               ReduceTimestepIfCloseToRLOF(dtSSE, m_Star2, m_Star1));

double BaseBinaryStar::ReduceTimestepIfCloseToRLOF(p_Dt, p_Primary, p_Secondary) {

    double newDt = p_Dt;     // default is unchanged

    double minSeparation = (m_SemiMajorAxisPrime * (1.0 - m_Eccentricity)); // AU
    double rocheRadius = minSeparation * (CalculateRocheLobeRadius_Static(p_Primary->Mass(), p_Secondary->Mass())); // AU
    double RLRatio = (p_Primary->Radius()*RSOL_TO_AU) / rocheRadius;

    if (utils::Compare(RLRatio, MASS_TRANSFER_THRESHOLD) >= 0) {                                                            

        if (p_Primary->DetermineEnvelopeType() == ENVELOPE::RADIATIVE) {
            newDt = m_FastPhaseCaseA ? p_Primary->CalculateThermalTimescale() 
                                      : p_Primary->CalculateDynamicalTimescale();  
        }
        else {
            newDt = p_Primary->CalculateThermalTimescale();
        }
    }
    else if (utils::Compare(RLRatio, ROCHE_LOBE_UPPER_THRESHOLD) >= 0) {
        newDt /= 10.0;
    }
    else if (utils::Compare(RLRatio, ROCHE_LOBE_LOWER_THRESHOLD) >= 0) {
        newDt /= 5.0;
    }

    return newDt;
}

double BaseStar::CalculateTimestepFromSSE() { // This replaces CalculateTimestep

    CalculateAllTimescales();
    CalculateGBParams();

    dtPhase = ChooseTimestep(m_Age); // Virtual function, depends on Stellar Type. Should rename to CalculateTimestepFromPhase() // should minimize  dt_k, dt_e, and dt_N

    dtRelative = ReduceTimestepIfRelativeChangesTooLarge(dtPhase) // Evolves copy star by given timestep, reduces if relative changes are too large // minimizes dt_ml, dt_R, dt_mc

    dtDyn = m_DynamicalTimescale // dt should never go below standard dynamical timescale of the star

    dtStar = std::max( dtRelative, dtDyn, ABSOLUTE_MINIMUM_TIMESTEP) 
}

BaseStar::ReduceTimestepIfRelativeChangesTooLarge(p_Dt) {

    BaseStar* starCopy = new BaseStar(*this);     // copy current star

    (void) starCopy->UpdateAttributesAndAgeOneTimestep(0.0, 0.0, p_Dt, p_Switch=true, p_ForceRecalculate=false);   // Really should just evolve the star according to SSE for one timestep

    // Calculate the maximal timescale at which relative differences do not exceed our thresholds

    // Minimize relative radial difference
    double Rdot = (starCopy->Radius()-starCopy->RadiusPrev())/p_Dt;  
    double dtR = std::abs((starCopy->RadiusPrev()/Rdot)) * MAXIMUM_RELATIVE_RADIAL_DIFFERENCE

    // Minimize relative mass loss difference
    double Mdot = (starCopy->Mass()-starCopy->MassPrev())/p_Dt;
    double dtM = std::abs((starCopy->MassPrev()/Rdot)) * MAXIMUM_RELATIVE_MASS_LOSS_DIFFERENCE

    // Minimize relative core mass loss difference - if stellar type changes, don't worry about this
    double dtMc; // default is a value not less than dtM, dtR
    if starCopy->StellarType() == starCopy->StellarTypePrev() {
        double Mcdot = (starCopy->MassCore()-starCopy->MassCorePrev())/p_Dt
        dtMc = std::abs((starCopy->MassCorePrev()/Rdot)) * MAXIMUM_RELATIVE_MASS_LOSS_DIFFERENCE
    }
    else {
        dtMc = dtM;
    }

    double dtRelative = std::min(dtR, dtM, dtMc, p_Dt); // Minimize to greatest constraint

    delete starCopy; starCopy = nullptr;                                      // nuke the copy

    return dtRelative;
}

Note: The call to UpdateAttributesAndAgeOneTimestep might not be used correctly. I think we should revisit this function (and its related versions) to be more clear about what is being done when it gets called. Here, it is intended to just age the copy star by the proposed timestep to see how the radius and mass parameters change over that timestep.

@ilyamandel I'm wondering where your proposed timestep re-evaluations fit into this. I believe all the timesteps, except dtk, are limited by physically motivations - but dtk is (as Hurley comments) more about resolving features in plots, e.g, the hook at the end of the main sequence. As this is only artificial, we may be able to incorporate your timestep increases here without affecting any of the events which need to be more closely resolved.

@reinhold-willcox -- just to let you know, I've already done away with some of the functions discussed above. We can still do better, and I agree with the general suggestion of moving to adaptive time-stepping, or at least allowing for variable time stepping to check for convergence.

Closing because some of this is out of date. We should take another look at Timesteps though, there's definitely still room for improvement.

TeamCOMPAS / COMPAS

Timesteps #290