Skip to content

Commit

Permalink
Merge pull request #1053 from TeamCOMPAS/BSEvsSSE
Browse files Browse the repository at this point in the history 
BSE vs SSE changes - see #1052
  • Loading branch information
@jeffriley
jeffriley authored Feb 20, 2024
2 parents 8193716 + 334a3cd commit 6d4ca01
Show file tree
Hide file tree
Showing 29 changed files with 607 additions and 230 deletions.
8 changes: 6 additions & 2 deletions online-docs/pages/User guide/Program options/program-options-list-defaults.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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|
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -1324,7 +1328,7 @@ Go to :ref:`the top of this page <options-props-top>` for the full alphabetical
**Administrative**

--mode, --number-of-systems, --evolve-double-white-dwarfs, --evolve-pulsars, --evolve-unbound-systems, --maximum-evolution-time, --maximum-number-timestep-iterations,
--random-seed, --timestep-multiplier
--random-seed, --timestep-multiplier, --timestep-filename

--grid, --grid-start-line, --grid-lines-to-process

Expand Down
63 changes: 63 additions & 0 deletions online-docs/pages/User guide/timestep-files.rst
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,63 @@
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.

If the timesteps provided in the timesteps file are completely consumed (i.e. each timestep in the timesteps file has been taken
during evolution), but evolution is not complete (i.e. no stopping condition for evolution has been met), evolution is stopped and
a warning issued, and the status of the star/system being evolved is set to ``"User-provided timesteps exhausted"``.

If the timesteps provided in the timesteps file are not completely consumed when evolution is complete (i.e. a stopping condition for
evolution has been met), a warning issued, and the status of the star/system evolved is set to ``"User-provided timesteps not consumed"``.

Timesteps files can have blank lines and/or comment lines. 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). An excerpt from a timesteps file is shown below::

# timesteps files can have blank lines and/or comment lines
# comment lines begin with the '#' character

# timesteps file for my experiment...
# description of file here
0.201713464000000
0.201935435000000

0.020216053100000
0.020218312700000
0.020220604600000
0.002430922190000
0.002430922190000

# The following timestep is the NUCLEAR_MINIMUM_TIMESTEP (10^6 * TIMESTEP_QUANTUM = 10^6 * 10^-12 = 10^-6 Myr)
0.000001000000000

0.000171981412000
0.000168541784000
0.000165170948000
0.000161867529000
0.000099674841500
0.000097681344700
0.000095727717800
0.000093813163400


COMPAS imposes a hard limit of ``1,000,000`` timesteps in a timesteps file.

1 change: 1 addition & 0 deletions online-docs/pages/User guide/user-guide.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
6 changes: 6 additions & 0 deletions online-docs/pages/whats-new.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
86 changes: 77 additions & 9 deletions src/BaseBinaryStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
}

unsigned long 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
Expand All @@ -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
}
Expand All @@ -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
}
Expand Down Expand Up @@ -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 (evolutionStatus == EVOLUTION_STATUS::CONTINUE && 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;
Expand Down
3 changes: 1 addition & 2 deletions src/BaseBinaryStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}

Expand Down
Loading

0 comments on commit 6d4ca01

Please sign in to comment.