Merge pull request #1361 from brcekadam/ChangeCalculateRadiusOnPhase

Fixes and enhancements to BRCEK core mass prescription

Author: jeffriley

src/BaseBinaryStar.cpp

@@ -2970,7 +2970,13 @@ double BaseBinaryStar::ChooseTimestep(const double p_Multiplier) {
         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
         }

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/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
 
 
 

src/MS_gt_07.h

@@ -51,7 +51,7 @@ class MS_gt_07: virtual public BaseStar, public MainSequence {
             m_InitialMainSequenceCoreMass = MainSequence::CalculateInitialMainSequenceCoreMass(m_MZAMS);
             m_MainSequenceCoreMass        = m_InitialMainSequenceCoreMass;
             m_Luminosity                  = MainSequence::CalculateLuminosityOnPhase(m_Age, m_Mass0, m_LZAMS0);
-            m_Radius                      = MainSequence::CalculateRadiusOnPhase(m_Mass, m_Age, m_RZAMS0);
+            m_Radius                      = MainSequence::CalculateRadiusOnPhase(m_Mass, m_Tau, m_RZAMS0);
             m_Temperature                 = BaseStar::CalculateTemperatureOnPhase_Static(m_Luminosity, m_Radius);
         }
     }