Merge pull request #1371 from jeffriley/sharpened

 Add option "--timestep-multipliers" to enable more granular, phase-dependent, timestep multipliers

Author: jeffriley

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

@@ -343,6 +343,7 @@ Default = 0.0
 
 **--debug-classes** |br|
 Developer-defined debug classes to enable (vector). |br|
+See :doc:`Vector program options <./program-options-vector-options>` for option format. |br|
 Default = `All debug classes enabled (e.g. no filtering)`
 
 **--debug-level** |br|
@@ -694,6 +695,7 @@ Default = HURLEY_ADD |br|
 
 **--log-classes** |br|
 Logging classes to be enabled (vector). |br|
+See :doc:`Vector program options <./program-options-vector-options>` for option format. |br|
 Default = `All debug classes enabled (e.g. no filtering)`
 
 **--logfile-common-envelopes** |br|
@@ -992,11 +994,13 @@ Default = SSE
 
 **--notes** |br|
 Annotation strings (vector). |br|
-Default = ""
+See :doc:`Vector program options <./program-options-vector-options>` for option format. |br|
+Default = "" for each annotation
 
 **--notes-hdrs** |br|
 Annotations header strings (vector). |br|
-Default = `No annotations`
+See :doc:`Vector program options <./program-options-vector-options>` for option format. |br|
+Default = `No annotation headers (no annotations)`
 
 **--number-of-systems [ -n ]** |br|
 The number of systems to simulate. |br|
@@ -1334,10 +1338,29 @@ User-defined timesteps filename. (See :doc:`Timestep files <../timestep-files>`)
 Default = ’’ (None)
 
 **--timestep-multiplier** |br|
-Multiplicative factor for timestep duration.  This multiplier is applied after the timesteps are chosen using other program options 
-such as ``--radial-change-fraction`` and ``--mass-change-fraction``, and will therefore override expected behaviour.  This option is 
-primarily intended for debugging/testing of convergence issues rather than for production runs. |br|
- Default = 1.0
+Multiplicative factor for timestep duration.  |br|
+|br|
+This multiplier is applied after the timesteps are chosen using other program options such as ``--radial-change-fraction`` 
+and ``--mass-change-fraction``, and will therefore override expected behaviour.  |br|
+This option can be used in conjunction with ``--timestep-multipliers``, in which case this multiplier, and the appropriate
+phase-dependent multiplier (specified by ``--timestep-multipliers``) are both applied. |br|
+Default = 1.0 |br| |br|
+This option is primarily intended for debugging/testing of convergence issues rather than for production runs. |br|
+
+**--timestep-multipliers** |br|
+Phase-dependent multiplicative factors for timestep duration. |br|
+See :doc:`Vector program options <./program-options-vector-options>` for option format. |br|
+A multicative factor can be specified for each phase (stellar type), where the ordinal value (zero-based) of the option value
+indicates the stellar type (from ``MS_LTE_07`` to ``CHEMICALLY_HOMOGENEOUS``, see stellar type list at 
+:doc:`../../Developer guide/Headers/typedefs-dot-h`>). |br|
+|br|
+This multiplier is applied after the timesteps are chosen using other program options such as ``--radial-change-fraction`` and 
+``--mass-change-fraction``, and will therefore override expected behaviour. |br|
+This option can be used in conjunction with ``--timestep-multiplier``, in which case that multiplier, and the appropriate
+phase-dependent multiplier (specified by ``--timestep-multipliers``) are both applied. |br|
+Default = 1.0 for each phase (stellar type) |br| |br|
+This option is primarily intended for debugging/testing of convergence issues rather than for production runs. |br|
+
 
 .. _options-props-U:
 

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

@@ -12,6 +12,7 @@ Currently, the only vector program options are:
 - --log-classes
 - --notes-hdrs
 - --notes
+- --timestep-multipliers
 
 
 The notation for vector program options provides for the specification of one or more values. e.g.::
@@ -26,8 +27,8 @@ for program options (to the command-line value, then to the COMPAS default) - le
 `log-class` had been left blank), and specifying an empty string ("") for a value would be ambiguous (as to whether the user wanted the
 option value to default, or just be an empty string).
 
-Option values (in general, but also specifically for vector options) may not begin with the dash character ('-'), because the shell parser
-will parse them as option names before passing them through to COMPAS.
+Option values beyond the first value may not begin with the dash character ('-'), because the shell parser will parse them as option names
+before passing them through to COMPAS (this is a Boost limitation).
 
 COMPAS imposes no limit to the length (number of characters) of an individual option values that are specified as strings, but there may 
 be practical limits imposed by the underlying system.

online-docs/pages/User guide/Running COMPAS/running-grid.rst

@@ -51,7 +51,7 @@ Both option ``--grid-start-line`` and ``--grid-lines-to-process`` are ignored if
 Example
 ~~~~~~~
 
-We will submit a set of COMPAS runs using a grid-file ``grid_demo.txt''
+We will submit a set of COMPAS runs using a grid-file ``grid_demo.txt``
 .. code-block::
 
     COMPAS --grid grid_demo.txt

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

@@ -140,7 +140,7 @@ Bonus Info
 ----------
 
 Dockerfile
-^^^^^^^^^^
+~~~~~~~~~~
 
 The `Dockerfile <https://docs.docker.com/engine/reference/builder/>`__ defines how the docker image is constructed.
 
@@ -159,7 +159,7 @@ The Dockerfile for COMPAS consists of 8 layers:
 Dockerfiles usually end with a `CMD` directive specifying the command to run when the container starts. COMPAS does not have a `CMD` directive because some users will run the executable directly, while others will use `runSubmit.py`.
 
 Makefile.docker
-^^^^^^^^^^^^^^^
+~~~~~~~~~~~~~~~
 
 A separate makefile is required for Docker to:
     1. Separate compiled files from source files.

online-docs/pages/whats-new.rst

@@ -3,10 +3,15 @@ What's new
 
 Following is a brief list of important updates to the COMPAS code.  A complete record of changes can be found in the file ``changelog.h``.
 
+**03.18.00 Apr 13, 2025**
+
+New command line option:
+
+* ``--timestep-multipliers`` to enable more granular, phase-dependent, timestep multipliers
+
 **03.17.00 Mar 22, 2025**
 
-* Added ENVELOPE_STATE_PRESCRIPTION::CONVECTIVE_MASS_FRACTION (default threshold of convective envelope by mass to label envelope 
-convective is 0.1, can be set with --convective-envelope-mass-threshold)
+* Added ENVELOPE_STATE_PRESCRIPTION::CONVECTIVE_MASS_FRACTION (default threshold of convective envelope by mass to label envelope convective is 0.1, can be set with --convective-envelope-mass-threshold)
 * Stable mass transfer now conserves angular momentum after accounting for the rotational angular momentum lost or gained by the stars
 * Imposed Keplerian rotation limit on mass-gaining stars:
 * Response depends on the new --response-to-spin-up option

src/BaseBinaryStar.cpp

@@ -1538,7 +1538,7 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() {
     m_Star2->SetPreCEEValues();                                                                                         // squirrel away pre CEE stellar values for star 2
   	SetPreCEEValues(semiMajorAxisRsol, eccentricity, rRLd1Rsol, rRLd2Rsol);                                             // squirrel away pre CEE binary values
 
-    m_MassTransferTimescale = MASS_TRANSFER_TIMESCALE::CE;
+    m_MassTransferTimescale = MT_TIMESCALE::CE;
     m_MassLossRateInRLOF    = DBL_MAX;
 
 	// double common envelope phase prescription (Brown 1995) to calculate new semi-major axis
@@ -2008,7 +2008,7 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) {
         jLoss = CalculateGammaAngularMomentumLoss();                                                                            // no - re-calculate angular momentum
     }
 
-    m_MassTransferTimescale         = MASS_TRANSFER_TIMESCALE::NONE;                                                            // initial reset
+    m_MassTransferTimescale         = MT_TIMESCALE::NONE;                                                                       // initial reset
     double betaThermal              = 0.0;                                                                                      // fraction of mass accreted if accretion proceeds on thermal timescale
     double maximumAccretionRate     = 0.0;                                                                                      // accretion rate if accretion proceeds on thermal timescale
     double donorMassLossRateThermal = m_Donor->CalculateThermalMassLossRate();
@@ -2037,17 +2037,17 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) {
         double zetaLobe        = CalculateZetaRocheLobe(jLoss, m_FractionAccreted);
         if (utils::Compare(zetaEquilibrium, zetaLobe) > 0  && massDiffDonor > 0.0) {                                            // yes, it's nuclear timescale mass transfer; no need for utils::Compare here
             m_MassLossRateInRLOF    = massDiffDonor / m_Dt;
-            m_MassTransferTimescale = MASS_TRANSFER_TIMESCALE::NUCLEAR;
+            m_MassTransferTimescale = MT_TIMESCALE::NUCLEAR;
             m_ZetaStar              = zetaEquilibrium;
             m_ZetaLobe              = zetaLobe;
         }
     }
-    if (m_MassTransferTimescale != MASS_TRANSFER_TIMESCALE::NUCLEAR) {                                                          // thermal timescale mass transfer (we will check for dynamically unstable / CE mass transfer later)
+    if (m_MassTransferTimescale != MT_TIMESCALE::NUCLEAR) {                                                                     // thermal timescale mass transfer (we will check for dynamically unstable / CE mass transfer later)
         m_ZetaLobe              = CalculateZetaRocheLobe(jLoss, betaThermal);
         m_ZetaStar              = m_Donor->CalculateZetaAdiabatic();
         m_MassLossRateInRLOF    = donorMassLossRateThermal;
         m_FractionAccreted      = betaThermal;
-        m_MassTransferTimescale = MASS_TRANSFER_TIMESCALE::THERMAL;
+        m_MassTransferTimescale = MT_TIMESCALE::THERMAL;
         massDiffDonor           = MassLossToFitInsideRocheLobe(this, m_Donor, m_Accretor, betaThermal, 0.0);                    // use root solver to determine how much mass should be lost from the donor to allow it to fit within the Roche lobe
     }
 
@@ -2090,7 +2090,7 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) {
         bool isEnvelopeRemoved = false;
 
         if (utils::Compare(m_Donor->CoreMass(), 0.0) > 0 && utils::Compare(envMassDonor, 0.0) > 0) {                            // donor has a core and an envelope
-            if (m_MassTransferTimescale == MASS_TRANSFER_TIMESCALE::THERMAL || utils::Compare (massDiffDonor, envMassDonor) >= 0) {
+            if (m_MassTransferTimescale == MT_TIMESCALE::THERMAL || utils::Compare (massDiffDonor, envMassDonor) >= 0) {
                 // remove entire envelope if thermal timescale MT from a giant or if the amount of necessary mass loss exceeds the envelope mass
                 massDiffDonor     = -envMassDonor;
                 isEnvelopeRemoved = true;
@@ -2270,7 +2270,7 @@ void BaseBinaryStar::InitialiseMassTransfer() {
 
 	m_MassTransferTrackerHistory = MT_TRACKING::NO_MASS_TRANSFER;	                                                            // Initiating flag, every timestep, to NO_MASS_TRANSFER. If it undergoes to MT or CEE, it should change.
 
-    m_MassTransferTimescale      = MASS_TRANSFER_TIMESCALE::NONE;
+    m_MassTransferTimescale      = MT_TIMESCALE::NONE;
     m_MassLossRateInRLOF         = 0.0;
 
     m_Star1->InitialiseMassTransfer(m_CEDetails.CEEnow, m_SemiMajorAxis, m_Eccentricity);                                       // initialise mass transfer for star1
@@ -2953,12 +2953,14 @@ void BaseBinaryStar::EmitGravitationalWave(const double p_Dt) {
  *
  * double ChooseTimestep(const double p_Multiplier)
  * 
- * @param   [IN]    p_Multiplier                timestep multiplier
+ * @param   [IN]    p_Factor                    factor applied to timestep (in addition to multipliers)
  * @return                                      new timestep in Myr
  */
-double BaseBinaryStar::ChooseTimestep(const double p_Multiplier) {
+double BaseBinaryStar::ChooseTimestep(const double p_Factor) {
 
-    double dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep());       // dt = smaller of timesteps required by individual stars
+    double dt1 = m_Star1->CalculateTimestep();
+    double dt2 = m_Star2->CalculateTimestep();
+    double dt  = std::min(dt1, dt2);                                                        // dt = smaller of timesteps required by individual stars
 
     if (!IsUnbound()) {                                                                     // check that binary is bound
 
@@ -3007,7 +3009,8 @@ double BaseBinaryStar::ChooseTimestep(const double p_Multiplier) {
         }
     }
 
-    dt *= p_Multiplier;	
+    // apply timestep multipliers
+    dt *= OPTIONS->TimestepMultiplier() * (dt1 < dt2 ? OPTIONS->TimestepMultipliers(static_cast<int>(m_Star1->StellarType())) : OPTIONS->TimestepMultipliers(static_cast<int>(m_Star2->StellarType()))) * p_Factor;
 
     return std::max(std::round(dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, TIDES_MINIMUM_FRACTIONAL_NUCLEAR_TIME * NUCLEAR_MINIMUM_TIMESTEP); // quantised and not less than minimum
 }
@@ -3240,7 +3243,7 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
                 }
 
                 // we want the first timestep to be small - calculate timestep and divide by 1000.0
-                dt = ChooseTimestep(OPTIONS->TimestepMultiplier() / 1000.0);                                                            // calculate timestep - make first step small
+                dt = ChooseTimestep(0.001);                                                                                             // calculate timestep - make first step small
             }
 
             unsigned long int stepNum = 1; 
@@ -3370,7 +3373,7 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
                             dt = timesteps[stepNum];
                         }
                         else {                                                                                                          // no - not using user-provided timesteps
-                            dt = ChooseTimestep(OPTIONS->TimestepMultiplier());
+                            dt = ChooseTimestep();
                         }
 
                         stepNum++;                                                                                                      // increment stepNum

src/BaseBinaryStar.h

@@ -344,7 +344,7 @@ class BaseBinaryStar {
     bool                m_MassTransfer;
     double              m_aMassTransferDiff;
 
-    MASS_TRANSFER_TIMESCALE m_MassTransferTimescale;
+    MT_TIMESCALE        m_MassTransferTimescale;
 
     MT_TRACKING         m_MassTransferTrackerHistory;
 
@@ -418,7 +418,7 @@ class BaseBinaryStar {
     void    CalculateGravitationalRadiation();
     void    EmitGravitationalWave(const double p_Dt);
 
-    double  ChooseTimestep(const double p_Multiplier);
+    double  ChooseTimestep(const double p_Factor = 1.0);
 
     void    CalculateEnergyAndAngularMomentum();
 

src/HG.cpp

@@ -134,7 +134,9 @@ double HG::CalculateLambdaLoveridge(const double p_EnvMass, const bool p_IsMassL
     logBindingEnergy += 33.29866;                                                           // + logBE0
     double bindingEnergy = PPOW(10.0, logBindingEnergy);
 
-    return utils::Compare(bindingEnergy, 0.0) > 0 && utils::Compare(p_EnvMass, 0.0) > 0 ? (G_CGS * m_Mass * MSOL_TO_G * p_EnvMass * MSOL_TO_G) / (m_Radius * RSOL_TO_AU * AU_TO_CM * bindingEnergy) : 1.0;     // default to 1.0 (usual lambda default) if binding energy is not sensible [should never happen] or if envelope mass is not positive [can be zero]
+    return utils::Compare(bindingEnergy, 0.0) > 0 && utils::Compare(p_EnvMass, 0.0) > 0 
+            ? (G_CGS * m_Mass * MSOL_TO_G * p_EnvMass * MSOL_TO_G) / (m_Radius * RSOL_TO_AU * AU_TO_CM * bindingEnergy) 
+            : 1.0;                                                                          // default to 1.0 (usual lambda default) if binding energy is not sensible [should never happen] or if envelope mass is not positive [can be zero]
 }