Merge pull request #1374 from jeffriley/sharpened

Requested changes

Author: jeffriley

compas_python_utils/preprocessing/compasConfigDefault.yaml

@@ -1,5 +1,5 @@
 ##~!!~## COMPAS option values
-##~!!~## File Created Tue Apr  1 17:29:43 2025 by COMPAS v03.17.00
+##~!!~## File Created Mon Apr 14 17:15:32 2025 by COMPAS v03.18.00
 ##~!!~##
 ##~!!~## The default COMPAS YAML file (``compasConfigDefault.yaml``), as distributed, has
 ##~!!~## all COMPAS option entries commented so that the COMPAS default value for the
@@ -88,6 +88,7 @@ numericalChoices:
 #    --number-of-systems: 10                                               # Default: 10                                                                                                                                                                               # number of systems per batch
 #    --radial-change-fraction: 0.100000                                    # Default: 0.100000                                                                                                                                                                         # approximate desired fractional changes in stellar radius per timestep
 #    --timestep-multiplier: 1.000000                                       # Default: 1.000000                                                                                                                                                                         # optional multiplier relative to default time step duration
+#    --timestep-multipliers: '{ }'                                         # Default: ''                                                                                                                                                                               # optional phase-dependent multipliers relative to default time step duration
 
     ### STELLAR PROPERTIES
 #    --cool-wind-mass-loss-multiplier: 1.000000                            # Default: 1.000000

online-docs/pages/User guide/COMPAS output/standard-logfiles-record-specification-stellar.rst

@@ -1789,22 +1789,6 @@ same header string.`
    * - Header Strings:
      - Recycled_NS, Recycled_NS(1), Recycled_NS(2), Recycled_NS(SN), Recycled_NS(CP)
 
-.. flat-table::
-   :widths: 25 75 1 1
-   :header-rows: 0
-   :class: aligned-text
-
-   * - :cspan:`2` **RLOF_ONTO_NS**
-     -
-   * - Data type:
-     - DOUBLE
-   * - COMPAS variable:
-     - `derived from` BaseStar::m_SupernovaDetails.events.past
-   * - Description:
-     - Flag to indicate whether the star transferred mass to a neutron star at any time prior to the current timestep.
-   * - Header Strings:
-     - RLOF->NS, RLOF->NS(1), RLOF->NS(2), RLOF->NS(SN), RLOF->NS(CP)
-
 .. flat-table::
    :widths: 25 75 1 1
    :header-rows: 0

online-docs/pages/whats-new.rst

@@ -3,29 +3,36 @@ 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**
+**03.18.00 Apr 14, 2025**
 
 New command line option:
 
 * ``--timestep-multipliers`` to enable more granular, phase-dependent, timestep multipliers
 
+**03.17.03 Apr 14, 2025**
+
+* Neutron stars are now labelled as ``RecycledNS`` when undergoing mass transfer through common envelope (when ``--neutron-star-accretion-in-ce`` is not set to ``ZERO``). 
+* Removed output option ``RLOF_ONTO_NS`` as it can be retrieved from existing RLOF output info. 
+
 **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
-* default (TRANSFER_TO_ORBIT) allows the star to accrete, but excess angular momentum is deposited in the orbit
-* KEPLERIAN_LIMIT forces mass  transfer to become non-conservative once star (approximately) reaches super-critical rotation
-* while the NO_LIMIT variation allows arbitrary super-critical accretion, to match legacy choices
+* Imposed Keplerian rotation limit on mass-gaining stars: response depends on the new ``--response-to-spin-up`` option, with possible values:
+   * ``TRANSFER_TO_ORBIT`` (default) allows the star to accrete, but excess angular momentum is deposited in the orbit
+   * ``KEPLERIAN_LIMIT`` forces mass transfer to become non-conservative once star (approximately) reaches super-critical rotation
+   * ``NO_LIMIT`` allows arbitrary super-critical accretion, to match legacy choices
 
 **03.16.02 Mar 19, 2025**
 
-New output options for supernova:
+New output options for supernova, which allow for full characterization of the binary orientation post-SN:
 
-* ORBITAL_ANGULAR_MOMENTUM_VECTOR_X, ORBITAL_ANGULAR_MOMENTUM_VECTOR_Y, ORBITAL_ANGULAR_MOMENTUM_VECTOR_Z,
-* SYSTEMIC_VELOCITY_X, SYSTEMIC_VELOCITY_Y, SYSTEMIC_VELOCITY_Z,
-* These allow for full characterization of the binary orientation post-SN
+* ORBITAL_ANGULAR_MOMENTUM_VECTOR_X
+* ORBITAL_ANGULAR_MOMENTUM_VECTOR_Y
+* ORBITAL_ANGULAR_MOMENTUM_VECTOR_Z
+* SYSTEMIC_VELOCITY_X
+* SYSTEMIC_VELOCITY_Y
+* SYSTEMIC_VELOCITY_Z
 
 **03.15.00 Mar 5, 2025**
 

src/BaseBinaryStar.cpp

@@ -2146,12 +2146,14 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) {
         }
     }
 
-	// Check for recycled pulsars. Not considering CEE as a way of recycling NSs.
-	if (!m_CEDetails.CEEnow && m_Accretor->IsOneOf({ STELLAR_TYPE::NEUTRON_STAR })) {                                           // accretor is a neutron star
-        m_Donor->SetRLOFOntoNS();                                                                                               // donor donated mass to a neutron star
+    // Check for recycled pulsars. Not considering CEE as a way of recycling NSs.
+    if (!m_CEDetails.CEEnow && m_Accretor->IsOneOf({ STELLAR_TYPE::NEUTRON_STAR })) {                                           // accretor is a neutron star, system is not in CE 
         m_Accretor->SetRecycledNS();                                                                                            // accretor is (was) a recycled NS
-	}
-    
+    }
+    else if (m_CEDetails.CEEnow && m_Accretor->IsOneOf({ STELLAR_TYPE::NEUTRON_STAR })
+             && OPTIONS->NeutronStarAccretionInCE() != NS_ACCRETION_IN_CE::ZERO) {                                              // accretor is a neutron star, system is in CE
+        m_Accretor->SetRecycledNS();                                                                                            // accretor is (was) a recycled NS
+    }
 }
 
 
@@ -2958,8 +2960,9 @@ void BaseBinaryStar::EmitGravitationalWave(const double p_Dt) {
  */
 double BaseBinaryStar::ChooseTimestep(const double p_Factor) {
 
-    double dt1 = m_Star1->CalculateTimestep();
-    double dt2 = m_Star2->CalculateTimestep();
+    double dt1 = m_Star1->CalculateTimestep() * OPTIONS->TimestepMultipliers(static_cast<int>(m_Star1->StellarType()));
+    double dt2 = m_Star2->CalculateTimestep() * OPTIONS->TimestepMultipliers(static_cast<int>(m_Star2->StellarType()));
+
     double dt  = std::min(dt1, dt2);                                                        // dt = smaller of timesteps required by individual stars
 
     if (!IsUnbound()) {                                                                     // check that binary is bound
@@ -2970,7 +2973,13 @@ double BaseBinaryStar::ChooseTimestep(const double p_Factor) {
         if ((utils::Compare(radiusToRL1 * (1.0 + 2.0 * OPTIONS->RadialChangeFraction()), 1.0) >= 0 && utils::Compare(radiusToRL1 * (1.0 + 0.5 * OPTIONS->RadialChangeFraction()), 1.0) <= 0) ||
             (utils::Compare(radiusToRL2 * (1.0 + 2.0 * OPTIONS->RadialChangeFraction()), 1.0) >= 0 && utils::Compare(radiusToRL2 * (1.0 + 0.5 * OPTIONS->RadialChangeFraction()), 1.0) <= 0))
             dt /= 2.0;
-
+        
+        // limit time step for stars losing mass on nuclear timescale
+        if (utils::Compare(radiusToRL1 * (1.0 + 0.5 * OPTIONS->RadialChangeFraction()), 1.0) > 0)
+            dt = std::min(dt, 0.5 * OPTIONS->RadialChangeFraction() * m_Star1->CalculateRadialExpansionTimescaleDuringMassTransfer());
+        if (utils::Compare(radiusToRL2 * (1.0 + 0.5 * OPTIONS->RadialChangeFraction()), 1.0) > 0)
+            dt = std::min(dt, 0.5 * OPTIONS->RadialChangeFraction() * m_Star2->CalculateRadialExpansionTimescaleDuringMassTransfer());
+        
         if (OPTIONS->EmitGravitationalRadiation()) {                                        // emitting GWs?
             dt = std::min(dt, -1.0E-2 * m_SemiMajorAxis / m_DaDtGW);                        // yes - reduce timestep if necessary to ensure that the orbital separation does not change by more than ~1% per timestep due to GW emission
         }
@@ -3009,8 +3018,7 @@ double BaseBinaryStar::ChooseTimestep(const double p_Factor) {
         }
     }
 
-    // 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;
+    dt *= OPTIONS->TimestepMultiplier() * 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
 }

src/BaseStar.cpp

@@ -2569,33 +2569,6 @@ double BaseStar::CalculateMassLossRate() {
 }
 
 
-/*
- * Calculate the nuclear mass loss rate as the mass divided by the radial expansion timescale
- * We do not use CalculateRadialExpansionTimescale(), however, since in the process of mass transfer the previous radius
- * is determined by binary evolution, not nuclear timescale evolution
- *
- *
- * double CalculateNuclearMassLossRate()
- *
- * @return                                      Nuclear mass loss rate
- */
-double BaseStar::CalculateNuclearMassLossRate() {
-    
-    // We create and age it slightly to determine how the radius will change.
-    // To be sure the clone does not participate in logging, we set its persistence to EPHEMERAL.
-    BaseStar *clone = Clone(OBJECT_PERSISTENCE::EPHEMERAL, false);                              // do not re-initialise the clone
-
-    double timestep = std::max(1000.0 * NUCLEAR_MINIMUM_TIMESTEP, m_Age / 1.0E6);
-    clone->UpdateAttributesAndAgeOneTimestep(0.0, 0.0, timestep, true, false);
-    double radiusAfterAging = clone->Radius();
-    delete clone; clone = nullptr;                                                              // return the memory allocated for the clone
-
-    double radialExpansionTimescale = timestep * m_Radius / fabs(m_Radius - radiusAfterAging);
-
-    return m_Mass / radialExpansionTimescale;
-}
-
-
 /*
  * Calculate values for mDot and mass assuming mass loss is applied
  *
@@ -3684,6 +3657,31 @@ double BaseStar::CalculateRadialExpansionTimescale_Static(const STELLAR_TYPE p_S
 }
 
 
+/*
+ * Calculate the radial expansion timescale in the mass transfer regime
+ * We do not use CalculateRadialExpansionTimescale(), since in the process of mass transfer the previous radius
+ * is determined by binary evolution, not nuclear timescale evolution
+ *
+ *
+ * double CalculateRadialExpansionTimescaleDuringMassTransfer()
+ *
+ * @return                                      Radial expansion timescale
+ */
+double BaseStar::CalculateRadialExpansionTimescaleDuringMassTransfer() {
+    
+    // We create and age it slightly to determine how the radius will change.
+    // To be sure the clone does not participate in logging, we set its persistence to EPHEMERAL.
+    BaseStar *clone = Clone(OBJECT_PERSISTENCE::EPHEMERAL, false);                              // do not re-initialise the clone
+
+    double timestep = std::max(1000.0 * NUCLEAR_MINIMUM_TIMESTEP, m_Age / 1.0E6);
+    clone->UpdateAttributesAndAgeOneTimestep(0.0, 0.0, timestep, true, false);
+    double radiusAfterAging = clone->Radius();
+    delete clone; clone = nullptr;                                                              // return the memory allocated for the clone
+
+    return timestep * m_Radius / fabs(m_Radius - radiusAfterAging);
+}
+
+
 /*
  * Calculate mass change timescale
  *

src/BaseStar.h

@@ -260,7 +260,7 @@ class BaseStar {
     virtual double          CalculateMomentOfInertia() const                                                    { return (0.1 * (m_Mass) * m_Radius * m_Radius); }                  // Defaults to MS. k2 = 0.1 as defined in Hurley et al. 2000, after eq 109
     virtual double          CalculateMomentOfInertiaAU() const                                                  { return CalculateMomentOfInertia() * RSOL_TO_AU * RSOL_TO_AU; }
 
-            double          CalculateNuclearMassLossRate();
+            double          CalculateNuclearMassLossRate()                                                      { return m_Mass / CalculateRadialExpansionTimescaleDuringMassTransfer(); }
 
             double          CalculateOmegaCHE(const double p_MZAMS, const double p_Metallicity) const;
 
@@ -270,6 +270,8 @@ class BaseStar {
 
             double          CalculateRadialExpansionTimescale() const                                           { return CalculateRadialExpansionTimescale_Static(m_StellarType, m_StellarTypePrev, m_Radius, m_RadiusPrev, m_DtPrev); } // Use class member variables
 
+            double          CalculateRadialExpansionTimescaleDuringMassTransfer();
+    
     virtual double          CalculateRadialExtentConvectiveEnvelope() const                                     { return 0.0; }                                                     // Default for stars with no convective envelope
 
     virtual double          CalculateRadiusOnMassChange(double p_dM)                                            { return Radius(); } // NO-OP

src/BinaryConstituentStar.cpp

@@ -79,7 +79,6 @@ COMPAS_VARIABLE BinaryConstituentStar::StellarPropertyValue(const T_ANY_PROPERTY
         case ANY_STAR_PROPERTY::RADIAL_EXPANSION_TIMESCALE_POST_COMMON_ENVELOPE:    value = RadialExpansionTimescalePostCEE();              break;
         case ANY_STAR_PROPERTY::RADIAL_EXPANSION_TIMESCALE_PRE_COMMON_ENVELOPE:     value = RadialExpansionTimescalePreCEE();               break;
         case ANY_STAR_PROPERTY::RECYCLED_NEUTRON_STAR:                              value = ExperiencedRecycledNS();                        break;
-        case ANY_STAR_PROPERTY::RLOF_ONTO_NS:                                       value = ExperiencedRLOFOntoNS();                        break;
         case ANY_STAR_PROPERTY::TEMPERATURE_POST_COMMON_ENVELOPE:                   value = TemperaturePostCEE() * TSOL;                    break;
         case ANY_STAR_PROPERTY::TEMPERATURE_PRE_COMMON_ENVELOPE:                    value = TemperaturePreCEE() * TSOL;                     break;
         case ANY_STAR_PROPERTY::THERMAL_TIMESCALE_POST_COMMON_ENVELOPE:             value = ThermalTimescalePostCEE();                      break;

src/BinaryConstituentStar.h

@@ -50,8 +50,7 @@ class BinaryConstituentStar: virtual public Star {
         m_CEDetails.postCEE.radialExpansionTimescale = DEFAULT_INITIAL_DOUBLE_VALUE;
 
         m_Flags.recycledNS                           = false;
-        m_Flags.rlofOntoNS                           = false;
-
+        
         m_MassLossDiff                               = DEFAULT_INITIAL_DOUBLE_VALUE;
         m_MassTransferDiff                           = DEFAULT_INITIAL_DOUBLE_VALUE;
 
@@ -151,7 +150,6 @@ class BinaryConstituentStar: virtual public Star {
 
     bool            ExperiencedRecycledNS() const                                       { return m_Flags.recycledNS; }
     bool            ExperiencedRLOF() const                                             { return m_RLOFDetails.experiencedRLOF; }
-    bool            ExperiencedRLOFOntoNS() const                                       { return m_Flags.rlofOntoNS; }
 
     double          HeCoreMassAtCEE() const                                             { return m_CEDetails.HeCoreMass; }
 
@@ -196,8 +194,6 @@ class BinaryConstituentStar: virtual public Star {
     void            ClearRecycledNS()                                                   { m_Flags.recycledNS = false; }
     void            SetRecycledNS()                                                     { m_Flags.recycledNS = true; }
 
-    void            ClearRLOFOntoNS()                                                   { m_Flags.rlofOntoNS = false; }
-    void            SetRLOFOntoNS()                                                     { m_Flags.rlofOntoNS = true; }
 
     void            CalculateCommonEnvelopeValues();
 
@@ -248,7 +244,6 @@ class BinaryConstituentStar: virtual public Star {
 
     struct FLAGS {                                                  // Miscellaneous flags
         bool recycledNS;                                            // Indicate whether the accretor was a recycled neutron star
-        bool rlofOntoNS;                                            // Indicates whether the donor donated mass to neutron star through RLOF
     }                       m_Flags;
 
     double                  m_MassLossDiff;

src/HG.h

@@ -90,14 +90,14 @@ class HG: virtual public BaseStar, public GiantBranch {
     double          CalculateHeliumAbundanceCoreOnPhase() const                     { return 1.0 - m_Metallicity; }                                                             // Use class member variables
 
     double          CalculateHeliumAbundanceSurfaceAtPhaseEnd() const               { return CalculateHeliumAbundanceSurfaceOnPhase(); }
-    double          CalculateHeliumAbundanceSurfaceOnPhase() const                  { return m_InitialHeliumAbundance; }                                                        // Use class member variables
+    double          CalculateHeliumAbundanceSurfaceOnPhase() const                  { return m_HeliumAbundanceSurface; }                                                        // Use class member variables
 
     double          CalculateHydrogenAbundanceCoreAtPhaseEnd() const                { return CalculateHydrogenAbundanceCoreOnPhase(); }
     double          CalculateHydrogenAbundanceCoreOnPhase(const double p_Tau) const;
     double          CalculateHydrogenAbundanceCoreOnPhase() const                   { return 0.0; }                                                                             // Star has exhausted hydrogen in its core
 
     double          CalculateHydrogenAbundanceSurfaceAtPhaseEnd() const             { return CalculateHydrogenAbundanceSurfaceOnPhase(); }
-    double          CalculateHydrogenAbundanceSurfaceOnPhase() const                { return m_InitialHydrogenAbundance; }                                                      // Use class member variables
+    double          CalculateHydrogenAbundanceSurfaceOnPhase() const                { return m_HydrogenAbundanceSurface; }                                                      // Use class member variables