Merge pull request #1052 from jeffriley/BSEvsSSE

BSE vs SSE investigation and changes

online-docs/pages/User guide/Program options/program-options-list-defaults.rst

@@ -413,7 +413,7 @@ Default = 5.75
 :ref:`Back to Top <options-props-top>`
 
 **--grid** |br|
-Grid filename. |br|
+Grid filename. (See :doc:`Grid files <../grid-files>`) |br|
 Default = ’’ (None)
 
 **--grid-lines-to-process** |br|
@@ -1161,6 +1161,10 @@ Default = FALSE
 
 :ref:`Back to Top <options-props-top>`
 
+**--timestep-filename** |br|
+User-defined timesteps filename. (See :doc:`Timestep files <../timestep-files>`) |br|
+Default = ’’ (None)
+
 **--timestep-multiplier** |br|
 Multiplicative factor for timestep duration. |br|
 Default = 1.0

online-docs/pages/User guide/timestep-files.rst

@@ -0,0 +1,28 @@
+Timestep files
+==============
+
+A timestep file allows users to specify the timesteps to be taken during Single Star Evolution (SSE) or Binary Star Evolution (BSE).
+
+If a timestep file is provided (via the ``timesteps-filename`` program option), the COMPAS code uses the user-provided timeteps to
+evolve stars/systems instead of calculating a new timestep at each iteration. The COMPAS code will read the first non-blank and
+non-comment line of the timesteps file and set the first timestep to that value, then for each iteration during evolution, each
+subsequent non-blank and non-comment line of the timesteps file is read and the timestep for that iteration set to the value read.
+
+Note that COMPAS will use the timestep exactly as read - the timesteps read from the timesteps file are not quantised by COMPAS, and 
+neither will they be truncated to any limit (e.g. ``ABSOLUTE_MINIMUM_TIMESTEP`` from ``constants.h``) or multiplied by any timestep
+multiplier set by the ``timestep-multiplier`` option.  Furthermore, no check will be made after the timestep is taken to limit the
+radial expansion of the star being evolved.
+
+Timesteps files can have blank lines and comments. Comments begin with a hash/pound character ('#') - the hash/pound character and text 
+beyond it are ignored by COMPAS.
+
+The ``timesteps-filename`` option is valid both on the command line and in a grid file.  If the ``timesteps-filename`` option is specified
+in conjunction with the ``--number-of-systems`` (``-n``) option, or with option ranges and/or sets, the same timeteps file specified will
+be used for each star or system being evolved. On the other hand, if the ``timesteps-filename`` option is specified in a grid file, the
+timesteps file specified on a grid file line will be read and used for the star or system evolved by that grid file line.
+
+Each line of a timesteps file which is not a blank line or comment line must contain a single, non-negative, floating-point number (in
+plain text - COMPAS will read the plain text and convert it to a floating-point number).
+
+COMPAS imposes a hard limit of ``1,000,000`` timesteps in a timesteps file.
+

online-docs/pages/User guide/user-guide.rst

@@ -10,6 +10,7 @@ This section contains the basic user guide for COMPAS.
    ./configuration
    ./Program options/program-options
    ./grid-files
+   ./timestep-files
    ./random-seed
    ./Running COMPAS/running-compas
    ./COMPAS output/output

online-docs/pages/whats-new.rst

@@ -6,6 +6,12 @@ Following is a brief list of important updates to the COMPAS code.  A complete r
 
 **LATEST RELEASE** |br|
 
+**02.42.00 Jan 04, 2023**
+
+* Timesteps are now quantised to an integral multiple of 1e-12Myr.
+* New option provided to allow user-defined timesteps: ``--timesteps-filename`` (See :doc:`Timestep files <../"User guide"/timestep-files>`).
+* Code changes to make SSE and BSE evolution more consistent (See `PR 1052 <https://github.com/TeamCOMPAS/COMPAS/pull/1052>`_).
+
 **02.41.03 Dec 28, 2023**
 
 * The functions ``BaseBinaryStar::CalculateAngularMomentum()``, ``BaseBinaryStar::CalculateTotalEnergy()``, and ``BaseStar::AngularMomentum()`` changed to use moment of inertia instead of gyration radius.

src/BaseBinaryStar.cpp

@@ -1453,7 +1453,7 @@ void BaseBinaryStar::EvaluateSupernovae() {
  * void ResolveCommonEnvelopeEvent()
  */
 void BaseBinaryStar::ResolveCommonEnvelopeEvent() {
-    
+
     double alphaCE = OPTIONS->CommonEnvelopeAlpha();                                                                    // CE efficiency parameter
 
 	double eccentricity      = Eccentricity();								                                            // current eccentricity (before CEE)
@@ -2407,10 +2407,37 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
         SAY("RandomSeed " << m_RandomSeed);
     }
 
+    bool usingProvidedTimesteps = false;                                                                                                    // using user-provided timesteps?
+    DBL_VECTOR timesteps;
+    if (!OPTIONS->TimestepsFileName().empty()) {                                                                                            // have timesteps filename?
+                                                                                                                                            // yes
+        ERROR error;
+        std::tie(error, timesteps) = utils::ReadTimesteps(OPTIONS->TimestepsFileName());                                                    // read timesteps from file
+        if (error != ERROR::NONE) {                                                                                                         // ok?
+            evolutionStatus = EVOLUTION_STATUS::NO_TIMESTEPS;                                                                               // no - set status
+            SHOW_WARN(error);                                                                                                               // show warning
+        }
+        else usingProvidedTimesteps = true;                                                                                                 // have user-provided timesteps
+    }
+
     if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) {                                                                                    // continue evolution
         // 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;                                     // initialise the timestep
-        int    stepNum = 1;                                                                                                                 // initialise step number
+
+        double dt;                                                                                                                          // timestep
+        if (usingProvidedTimesteps) {                                                                                                       // user-provided timesteps?
+            // get new timestep
+            //   - don't quantise
+            //   - don't apply timestep multiplier
+            // (we assume user wants the timesteps in the file)
+            dt = timesteps[0];
+        }
+        else {                                                                                                                              // no - not using user-provided timesteps
+            dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep()) * OPTIONS->TimestepMultiplier() / 1000.0;             // calculate timestep - make first step small
+            dt = std::max(std::round(dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, NUCLEAR_MINIMUM_TIMESTEP);                                  // quantised
+        }
+
+        long unsigned int stepNum = 1;                                                                                                      // initialise step number
+
         while (evolutionStatus == EVOLUTION_STATUS::CONTINUE) {                                                                             // perform binary evolution - iterate over timesteps until told to stop
 
             EvolveOneTimestep(dt);                                                                                                          // evolve the binary system one timestep
@@ -2419,6 +2446,11 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
             if (m_Error != ERROR::NONE) {                                                                                                   // SSE error for either constituent star?
                 evolutionStatus = EVOLUTION_STATUS::SSE_ERROR;                                                                              // yes - stop evolution
             }
+
+            // check for reasons to not continue evolution
+            else if (!OPTIONS->EvolveDoubleWhiteDwarfs() && IsWDandWD()) {                                                                  // double WD and their evolution is not enabled?
+                evolutionStatus = EVOLUTION_STATUS::WD_WD;                                                                                  // yes - do not evolve double WD systems
+            }
             else if (HasOneOf({ STELLAR_TYPE::MASSLESS_REMNANT })) {                                                                        // at least one massless remnant?
                 evolutionStatus = EVOLUTION_STATUS::MASSLESS_REMNANT;                                                                       // yes - stop evolution
             }
@@ -2444,7 +2476,7 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
 
                 (void)PrintRLOFParameters();                                                                                                // print (log) RLOF parameters
 
-                // check for problems
+                // check for problems/reasons to not continue evolution
                 if (StellarMerger()) {                                                                                                      // have stars merged?
                     evolutionStatus = EVOLUTION_STATUS::STELLAR_MERGER;                                                                     // for now, stop evolution
                 }
@@ -2486,37 +2518,73 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
                         }
                     }
 
-                    // check for problems
+                    // check whether to continue evolution
                     if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) {                                                                    // continue evolution?
+
+                        // check for problems
                         if (m_Error != ERROR::NONE) {                                                                                       // error in binary evolution?
                             evolutionStatus = EVOLUTION_STATUS::BINARY_ERROR;                                                               // yes - stop evolution
                         }
-                        else if (!OPTIONS->EvolveDoubleWhiteDwarfs() && IsWDandWD()) {                                                      // double WD and their evolution is not enabled?
-                            evolutionStatus = EVOLUTION_STATUS::WD_WD;                                                                      // yes - do not evolve double WD systems
-                        }
+
+                        // check for reasons to stop evolution
                         else if (IsDCO() && m_Time > (m_DCOFormationTime + m_TimeToCoalescence) && !IsUnbound()) {                          // evolution time exceeds DCO merger time?
                             evolutionStatus = EVOLUTION_STATUS::DCO_MERGER_TIME;                                                            // yes - stop evolution
                         }
                         else if (m_Time > OPTIONS->MaxEvolutionTime()) {                                                                    // evolution time exceeds maximum?
                             evolutionStatus = EVOLUTION_STATUS::TIMES_UP;                                                                   // yes - stop evolution
                         }
+                        else if (!OPTIONS->EvolveDoubleWhiteDwarfs() && IsWDandWD()) {                                                      // double WD and their evolution is not enabled?
+                            evolutionStatus = EVOLUTION_STATUS::WD_WD;                                                                      // yes - do not evolve double WD systems
+                        }
+                        else if (HasOneOf({ STELLAR_TYPE::MASSLESS_REMNANT })) {                                                            // at least one massless remnant?
+                            evolutionStatus = EVOLUTION_STATUS::MASSLESS_REMNANT;                                                           // yes - stop evolution
+                        }
+                        else if (StellarMerger()) {                                                                                         // have stars merged?
+                            evolutionStatus = EVOLUTION_STATUS::STELLAR_MERGER;                                                             // for now, stop evolution
+                        }
+                        else if (HasStarsTouching()) {                                                                                      // binary components touching? (should usually be avoided as MT or CE or merger should happen prior to this)
+                            evolutionStatus = EVOLUTION_STATUS::STARS_TOUCHING;                                                             // yes - stop evolution
+                        }
+                        else if (IsUnbound() && !OPTIONS->EvolveUnboundSystems()) {                                                         // binary is unbound and we don't want unbound systems?
+                            m_Unbound       = true;                                                                                         // yes - set the unbound flag (should already be set)
+                            evolutionStatus = EVOLUTION_STATUS::UNBOUND;                                                                    // stop evolution
+                        }
                     }
                 }
             }
 
             (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::TIMESTEP_COMPLETED);                                                  // print (log) detailed output: this is after all changes made in the timestep
 
             if (stepNum >= OPTIONS->MaxNumberOfTimestepIterations()) evolutionStatus = EVOLUTION_STATUS::STEPS_UP;                          // number of timesteps for evolution exceeds maximum
+            else if (usingProvidedTimesteps && stepNum >= timesteps.size()) {
+                evolutionStatus = EVOLUTION_STATUS::TIMESTEPS_EXHAUSTED;                                                                    // using user-provided timesteps and all consumed
+                SHOW_WARN(ERROR::TIMESTEPS_EXHAUSTED);                                                                                      // show warning
+            }
 
             if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) {                                                                            // continue evolution?
+                if (usingProvidedTimesteps) {                                                                                               // user-provided timesteps
+                    // get new timestep
+                    //   - don't quantise
+                    //   - don't apply timestep multiplier
+                    // (we assume user wants the timesteps in the file)
+                    dt = timesteps[stepNum];
+                }
+                else {                                                                                                                      // not using user-provided timesteps
+                    dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep()) * OPTIONS->TimestepMultiplier();              // calculate new timestep
+                    dt = std::max(std::round(dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, NUCLEAR_MINIMUM_TIMESTEP);                          // quantised
+                }
 
-                dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep()) * OPTIONS->TimestepMultiplier();                  // new timestep
                 if ((m_Star1->IsOneOf({ STELLAR_TYPE::MASSLESS_REMNANT }) || m_Star2->IsOneOf({ STELLAR_TYPE::MASSLESS_REMNANT })) || dt < NUCLEAR_MINIMUM_TIMESTEP) {
                     dt = NUCLEAR_MINIMUM_TIMESTEP;                                                                                          // but not less than minimum
 		        }
                 stepNum++;                                                                                                                  // increment stepNum
             }
         }
+
+        if (usingProvidedTimesteps && timesteps.size() > stepNum) {                                                                         // all user-defined timesteps consumed?
+            evolutionStatus = EVOLUTION_STATUS::TIMESTEPS_NOT_CONSUMED;                                                                     // no - set status
+            SHOW_WARN(ERROR::TIMESTEPS_NOT_CONSUMED);                                                                                       // show warning
+        }
     }
 
     m_EvolutionStatus = evolutionStatus;

src/BaseBinaryStar.h

@@ -495,8 +495,7 @@ class BaseBinaryStar {
         return LOGGING->LogBSESystemParameters(this, p_RecordType);
     }
 
-    bool PrintDetailedOutput(const long int p_Id, 
-                             const BSE_DETAILED_RECORD_TYPE p_RecordType) const {
+    bool PrintDetailedOutput(const long int p_Id, const BSE_DETAILED_RECORD_TYPE p_RecordType) const {
         return OPTIONS->DetailedOutput() ? LOGGING->LogBSEDetailedOutput(this, p_Id, p_RecordType) : true;
     }