A range of fixes: using core radii for CE survival, fixes to #1255, #1258, #1265

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

@@ -847,7 +847,7 @@ Options: { JEANS, ISOTROPIC, CIRCUMBINARY, MACLEOD_LINEAR, ARBITRARY } |br|
 Default = ISOTROPIC
 
 **--mass-transfer-fa** |br|
-Mass Transfer fraction accreted. |br|
+Mass Transfer fraction accreted (beta). |br|
 Used when ``--mass-transfer-accretion-efficiency-prescription = FIXED_FRACTION``. |br|
 Default = 0.5
 
@@ -1122,10 +1122,11 @@ Enable mass loss due to pulsational-pair-instability (PPI). |br|
 Default = TRUE
 
 **--pulsational-pair-instability-prescription** |br|
-Pulsational pair instability prescription. |br|
+Pulsational pair instability prescription (only relevant when using ``--pulsational-pair-instability``). |br|
 Options: { HENDRIKS, COMPAS, STARTRACK, MARCHANT, FARMER } |br|
 ``HENDRIKS`` implements the prescription from Hendriks et al. 2023 |br|
-``COMPAS``, ``STARTRACK`` and ``MARCHANT`` follow Woosley 2017, Belczynski et al. 2016, and Marchant et al. 2018, all as implemented in Stevenson et al. 2019. |br|
+``COMPAS``, ``STARTRACK`` and ``MARCHANT`` follow Woosley 2017, Belczynski et al. 2016, and Marchant et al. 2018, 
+all as implemented in Stevenson et al. 2019. |br|
 ``FARMER`` follows Farmer et al. 2019 |br|
 Default = MARCHANT
 

online-docs/pages/whats-new.rst

@@ -3,6 +3,13 @@ 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.09.00 Nov 25, 2024**
+
+Improved nuclear timescale mass transfer: the nuclear timescale mass transfer rate is now set by the requirement that the star 
+ends the time step just filling its Roche lobe.
+Fixed several significant mass-transfer issues, such as accretors not gaining mass appropriately and failures
+in the root solver for fitting the star into the Roche lobe that were leading to artificial common envelopes and mergers.
+
 **03.08.02 Nov 18, 2024**
 
 Updated implementation of the mass transfer stability critical mass ratio tables from the team of Hongwei Ge.

src/BaseStar.cpp

@@ -2478,6 +2478,8 @@ double BaseStar::CalculateMassLossRateWolfRayetSanderVink2020(const double p_Mu)
  */
 double BaseStar::CalculateMassLossRateWolfRayetTemperatureCorrectionSander2023(const double p_Mdot) const {
 
+    if (p_Mdot <= 0.0) return 0.0;                                  // nothing to adjust
+
     const double teffRef = 141.0E3;                                 // reference effective temperature in Kelvin
     const double teffMin = 100.0E3;                                 // minimum effective temperature in Kelvin to apply correction
 
@@ -2493,6 +2495,8 @@ double BaseStar::CalculateMassLossRateWolfRayetTemperatureCorrectionSander2023(c
         logMdotCorrected = logMdotUncorrected;
     }
 
+    std::cout<<logMdotCorrected<<std::endl;
+
     return PPOW(10.0, logMdotCorrected);
 }
 

src/BaseStar.h

@@ -281,6 +281,9 @@ class BaseStar {
     virtual double          CalculateRadialExtentConvectiveEnvelope() const                                     { return 0.0; }                                                     // Default for stars with no convective envelope
 
     virtual double          CalculateRadiusOnPhaseTau(const double p_Mass, const double p_Tau) const            { return 0.0; }                                                     // Only defined for MS stars
+
+    virtual double          CalculateRemnantRadius() const                                                      { return Radius(); }    // relevant for MS stars, over-written for GB stars
+
 
             void            CalculateSNAnomalies(const double p_Eccentricity);
 
@@ -340,6 +343,8 @@ class BaseStar {
             void            StashSupernovaDetails(const STELLAR_TYPE p_StellarType,
                                                   const SSE_SN_RECORD_TYPE p_RecordType = SSE_SN_RECORD_TYPE::DEFAULT) { LOGGING->StashSSESupernovaDetails(this, p_StellarType, p_RecordType); }
 
+    virtual double          TAMSCoreMass() const                                                                { return 0.0; }                                                                                                             // except MS stars
+
     virtual void            UpdateAfterMerger(double p_Mass, double p_HydrogenMass) { }                                                                                             // Default is NO-OP
     virtual void            UpdateAgeAfterMassLoss() { }                                                                                                                            // Default is NO-OP
 

src/GiantBranch.h

@@ -20,9 +20,8 @@ class GiantBranch: virtual public BaseStar, public MainSequence {
     GiantBranch(){};
 
     GiantBranch(const BaseStar &p_BaseStar) : BaseStar(p_BaseStar), MainSequence(p_BaseStar) {}
-
-    virtual double          CalculateRemnantRadius() const;
-
+
+    double          CalculateRemnantRadius() const;
 
 protected:
 

src/HG.cpp

@@ -2,6 +2,7 @@
 #include "HeMS.h"
 #include "HeWD.h"
 
+
 ///////////////////////////////////////////////////////////////////////////////////////
 //                                                                                   //
 //                     COEFFICIENT AND CONSTANT CALCULATIONS ETC.                    //
@@ -861,6 +862,7 @@ double HG::CalculateRadiusOnPhase(const double p_Mass, const double p_Tau, const
 }
 
 
+
 ///////////////////////////////////////////////////////////////////////////////////////
 //                                                                                   //
 //                                 MASS CALCULATIONS                                 //

src/HG.h

@@ -11,52 +11,49 @@
 #include "GiantBranch.h"
 
 
-// JR: todo: revisit this one day - sometimes HG works better as GiantBranch, sometimes not...
-// Right now it is GiantBranch - figure out which is more appropriate
-
 class BaseStar;
 class GiantBranch;
 
 class HG: virtual public BaseStar, public GiantBranch {
-
+    
 public:
-
+    
     HG() { m_StellarType = STELLAR_TYPE::HERTZSPRUNG_GAP; };
 
     HG(const BaseStar &p_BaseStar, const bool p_Initialise = true) : BaseStar(p_BaseStar), GiantBranch(p_BaseStar) {
-        m_StellarType = STELLAR_TYPE::HERTZSPRUNG_GAP;                                                                                                                          // Set stellar type 
+        m_StellarType = STELLAR_TYPE::HERTZSPRUNG_GAP;                                                                                                                          // Set stellar type
         if (p_Initialise) Initialise();                                                                                                                                         // Initialise if required
     }
-
+    
     HG* Clone(const OBJECT_PERSISTENCE p_Persistence, const bool p_Initialise = true) {
-        HG* clone = new HG(*this, p_Initialise); 
-        clone->SetPersistence(p_Persistence); 
-        return clone; 
+        HG* clone = new HG(*this, p_Initialise);
+        clone->SetPersistence(p_Persistence);
+        return clone;
     }
-
+    
     static HG* Clone(HG& p_Star, const OBJECT_PERSISTENCE p_Persistence, const bool p_Initialise = true) {
-        HG* clone = new HG(p_Star, p_Initialise); 
-        clone->SetPersistence(p_Persistence); 
-        return clone; 
+        HG* clone = new HG(p_Star, p_Initialise);
+        clone->SetPersistence(p_Persistence);
+        return clone;
     }
-
+    
     MT_CASE DetermineMassTransferTypeAsDonor() const { return MT_CASE::B; }                                                                                                     // Always case B
-
-
+    
+    
 protected:
-
+    
     void Initialise() {
-
+        
         m_Tau = 0.0;                                                                                                                                                            // Start of phase
-
+        
         // update stellar properties at start of HG phase (since core definition changes)
         CalculateGBParams();
         CalculateTimescales();
-                                                                                                                                                                                // Initialise timescales
+        // Initialise timescales
         m_Age = m_Timescales[static_cast<int>(TIMESCALE::tMS)];                                                                                                                 // Set age appropriately
-
+        
         // update effective "initial" mass (m_Mass0) so that the core mass is at least equal to the minimum core mass but no more than total mass
-        // (only relevant if RetainCoreMassDuringCaseAMassTransfer()) 
+        // (only relevant if RetainCoreMassDuringCaseAMassTransfer())
         if (utils::Compare(CalculateCoreMassOnPhase(m_Mass0, m_Age), std::min(m_Mass, MinimumCoreMass())) < 0) {
             double desiredCoreMass = std::min(m_Mass, MinimumCoreMass());                                                                                                       // desired core mass
             m_Mass0 = Mass0ToMatchDesiredCoreMass(this, desiredCoreMass);                                                                                                       // use root finder to find new core mass estimate
@@ -69,56 +66,56 @@ class HG: virtual public BaseStar, public GiantBranch {
         }
         EvolveOnPhase(0.0);
     }
-
-
+    
+    
     // member functions - alphabetically
     double          CalculateCOCoreMassAtPhaseEnd() const                           { return 0.0; }                                                                             // McCO(HG) = 0.0
     double          CalculateCOCoreMassOnPhase() const                              { return 0.0; }                                                                             // McCO(HG) = 0.0
-
+    
     double          CalculateCoreMassAt2ndDredgeUp(const DBL_VECTOR &p_GBParams)    { return p_GBParams[static_cast<int>(GBP::McDU)]; }                                         // NO-OP
     double          CalculateCoreMassAtPhaseEnd(const double p_Mass) const;
     double          CalculateCoreMassAtPhaseEnd() const                             { return CalculateCoreMassAtPhaseEnd(m_Mass0); }                                            // Use class member variables
     double          CalculateCoreMassOnPhase(const double p_Mass, const double p_Time) const;
     double          CalculateCoreMassOnPhase() const                                { return CalculateCoreMassOnPhase(m_Mass0, m_Age); }                                        // Use class member variables
     double          CalculateCoreMassOnPhaseIgnoringPreviousCoreMass(const double p_Mass, const double p_Time) const;                                                           //  Ignore previous core mass constraint when computing expected core mass
-
+    
     double          CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const;
     double          CalculateCriticalMassRatioHurleyHjellmingWebbink() const        { return 0.25; }                                                                            // As coded in BSE. Using the inverse owing to how qCrit is defined in COMPAS. See Hurley et al. 2002 sect. 2.6.1 for additional details.
-
+    
     double          CalculateHeCoreMassAtPhaseEnd() const                           { return m_CoreMass; }                                                                      // McHe(HG) = Core Mass
     double          CalculateHeCoreMassOnPhase() const                              { return m_CoreMass; }                                                                      // McHe(HG) = Core Mass
 
-    double          CalculateHeliumAbundanceCoreAtPhaseEnd() const                  { return 1.0 - m_Metallicity; }                                         
-    double          CalculateHeliumAbundanceCoreOnPhase() const                     { return 1.0 - m_Metallicity; }                                                             // Use class member variables                                       
+    double          CalculateHeliumAbundanceCoreAtPhaseEnd() const                  { return 1.0 - m_Metallicity; }
+    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_InitialHeliumAbundance; }                                                        // 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          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          CalculateHydrogenAbundanceSurfaceAtPhaseEnd() const             { return CalculateHydrogenAbundanceSurfaceOnPhase(); }
     double          CalculateHydrogenAbundanceSurfaceOnPhase() const                { return m_InitialHydrogenAbundance; }                                                      // Use class member variables
 
-
-
+    
+    
     double          CalculateLambdaDewi() const;
     double          CalculateLambdaNanjingStarTrack(const double p_Mass, const double p_Metallicity) const;
     double          CalculateLambdaNanjingEnhanced(const int p_MassIndex, const STELLAR_POPULATION p_StellarPop) const;
-
+    
     double          CalculateLuminosityAtPhaseEnd(const double p_Mass) const;
     double          CalculateLuminosityAtPhaseEnd() const                           { return CalculateLuminosityAtPhaseEnd(m_Mass0);}                                           // Use class member variables
     double          CalculateLuminosityOnPhase(const double p_Age, const double p_Mass) const;
     double          CalculateLuminosityOnPhase() const                              { return CalculateLuminosityOnPhase(m_Age, m_Mass0); }                                      // Use class member variables
-
+    
     double          CalculateMassTransferRejuvenationFactor()                       { return 1.0; }
-
+    
     double          CalculateRadiusAtPhaseEnd(const double p_Mass) const;
     double          CalculateRadiusAtPhaseEnd() const                               { return CalculateRadiusAtPhaseEnd(m_Mass); }                                               // Use class member variables
     double          CalculateRadiusOnPhase(const double p_Mass, const double p_Tau, const double p_RZAMS) const;
     double          CalculateRadiusOnPhase() const                                  { return CalculateRadiusOnPhase(m_Mass0, m_Tau, m_RZAMS0); }                                // Use class member variables
-
+    
     double          CalculateRho(const double p_Mass) const;
 
     double          CalculateTauAtPhaseEnd() const                                  { return 1.0; }                                                                             // tau = 1.0 at end of HG
@@ -141,6 +138,7 @@ class HG: virtual public BaseStar, public GiantBranch {
     bool            ShouldEvolveOnPhase() const                                     { return (utils::Compare(m_Age, m_Timescales[static_cast<int>(TIMESCALE::tBGB)]) < 0); }    // Evolve on HG phase if age < Base Giant Branch timescale
     bool            ShouldSkipPhase() const                                         { return false; }                                                                           // Never skip HG phase
 
+    double          TAMSCoreMass() const                                            { return 0.0; }
     void            UpdateAfterMerger(double p_Mass, double p_HydrogenMass) { }                                                                                                 // Nothing to do for stars beyond the Main Sequence for now
     void            UpdateAgeAfterMassLoss();                                                                                                                                   // Per Hurley et al. 2000, section 7.1
 

src/HeHG.cpp

@@ -1,6 +1,7 @@
 #include "HeHG.h"
 #include "HeGB.h"
 #include "COWD.h"
+#include "HeWD.h"
 
 
 /*
@@ -234,6 +235,21 @@ double HeHG::CalculatePerturbationMu() const {
     return std::max(5.0 * ((McMax - m_CoreMass) / McMax), 0.0);         //return non-negative value to avoid round-off issues
 }
 
+/*
+ * Calculate radius of the remnant the star would become if it lost all of its
+ * envelope immediately (i.e. M = Mc, coreMass)
+ *
+ * Hurley et al. 2000, end of section 6
+ *
+ *
+ * double CalculateRemnantRadius()
+ *
+ * @return                                      Radius of remnant core in Rsol
+ */
+double HeHG::CalculateRemnantRadius() const {
+    return HeWD::CalculateRadiusOnPhase_Static(m_CoreMass);
+}
+
 
 ///////////////////////////////////////////////////////////////////////////////////////
 //                                                                                   //

src/HeHG.h

@@ -40,6 +40,7 @@ class HeHG: virtual public BaseStar, public HeMS {
     double  CalculateConvectiveCoreRadius () const                                                          { return std::min(5.0 * CalculateRemnantRadius(), m_Radius); }          // Last paragraph of section 6 of Hurley+ 2000
     static void CalculateGBParams_Static(const double p_Mass0, const double p_Mass, const double p_LogMetallicityXi, const DBL_VECTOR &p_MassCutoffs, const DBL_VECTOR &p_AnCoefficients, const DBL_VECTOR &p_BnCoefficients, DBL_VECTOR &p_GBParams);
 
+    double  CalculateRemnantRadius() const;
 
 protected:
 

src/HeMS.h

@@ -52,6 +52,8 @@ class HeMS: virtual public BaseStar, public TPAGB {
     static DBL_DBL  CalculateRadiusAtPhaseEnd_Static(const double p_Mass, const double p_Luminosity);
     static double   CalculateRadiusAtZAMS_Static(const double p_Mass);
     static double   CalculateRadiusOnPhase_Static(const double p_Mass, const double p_Tau);
+
+           double   CalculateRemnantRadius() const                                                                  { return Radius(); }
 
     MT_CASE         DetermineMassTransferTypeAsDonor() const                                                        { return MT_CASE::OTHER; }                                                      // Not A, B, C, or NONE