Resolve sudden mass drop at end of HG

Author: ilyamandel

src/BaseBinaryStar.cpp

@@ -2132,14 +2132,14 @@ void BaseBinaryStar::ResolveMassChanges() {
 
     // update mass of star1 according to mass loss and mass transfer, then update age accordingly
     (void)m_Star1->UpdateAttributes(m_Star1->MassPrev() - m_Star1->Mass() + m_Star1->MassLossDiff() + m_Star1->MassTransferDiff(), 0.0);    // update mass for star1
-    m_Star1->UpdateInitialMass();                                                                       // update initial mass of star1 (MS, HG & HeMS)  JR: todo: fix this kludge one day - mass0 is overloaded, and isn't always "initial mass"
+    m_Star1->UpdateInitialMass();                                                                       // update effective initial mass of star1 (MS, HG & HeMS, HeHG)
     m_Star1->UpdateAgeAfterMassLoss();                                                                  // update age of star1
     m_Star1->ApplyMassTransferRejuvenationFactor();                                                     // apply age rejuvenation factor for star1
     m_Star1->UpdateAttributes(0.0, 0.0, true);
 
     // rinse and repeat for star2
     (void)m_Star2->UpdateAttributes(m_Star2->MassPrev() - m_Star2->Mass() + m_Star2->MassLossDiff() + m_Star2->MassTransferDiff(), 0.0);    // update mass for star2
-    m_Star2->UpdateInitialMass();                                                                       // update initial mass of star 2 (MS, HG & HeMS)  JR: todo: fix this kludge one day - mass0 is overloaded, and isn't always "initial mass"
+    m_Star2->UpdateInitialMass();                                                                       // update effective initial mass of star 2 (MS, HG & HeMS, HeHG)
     m_Star2->UpdateAgeAfterMassLoss();                                                                  // update age of star2
     m_Star2->ApplyMassTransferRejuvenationFactor();                                                     // apply age rejuvenation factor for star2
     m_Star2->UpdateAttributes(0.0, 0.0, true);

src/BaseBinaryStar.h

@@ -543,7 +543,7 @@ class BaseBinaryStar {
             (void)donorCopy->UpdateAttributes(-dM, -dM*donorCopy->Mass0()/donorCopy->Mass());
 
             // Modify donor Mass0 and Age for MS (including HeMS) and HG stars
-            donorCopy->UpdateInitialMass();         // update initial mass (MS, HG & HeMS)  JR: todo: fix this kludge - mass0 is overloaded, and isn't always "initial mass"
+            donorCopy->UpdateInitialMass();         // update initial mass (MS, HG & HeMS)  
             donorCopy->UpdateAgeAfterMassLoss();    // update age (MS, HG & HeMS)
 
             (void)donorCopy->AgeOneTimestep(0.0);   // recalculate radius of star - don't age - just update values

src/BaseStar.cpp

@@ -1055,7 +1055,7 @@ double BaseStar::CalculateLogBindingEnergyLoveridge(bool p_IsMassLoss) const {
         logBindingEnergy += lCoefficients.alpha_mr * utils::intPow(log10(m_Mass), lCoefficients.m) * utils::intPow(log10(m_Radius + deltaR), lCoefficients.r);
     }
 
-    double MZAMS_Mass = (m_MZAMS - m_Mass) / m_MZAMS;
+    double MZAMS_Mass = (m_MZAMS - m_Mass) / m_MZAMS;                                       // Should m_ZAMS really be m_Mass0 (i.e., account for change in effective mass through mass loss in winds, MS mass transfer?)
     logBindingEnergy *= p_IsMassLoss ? 1.0 + (0.25 * MZAMS_Mass * MZAMS_Mass) : 1.0;        // apply mass-loss correction factor (lambda)
 
     constexpr double logBE0 = 33.29866;                                                     // JR: todo: what is this for?  Should it be in constants.h?
@@ -1239,7 +1239,7 @@ double BaseStar::CalculateLuminosityAtZAMS(const double p_MZAMS) {
  *
  * double CalculateLuminosityAtBAGB(double p_Mass)
  *
- * @param   [IN]    p_Mass                      Mass in Msol
+ * @param   [IN]    p_Mass                      (Effective) mass in Msol
  * @return                                      Luminosity at BAGB in Lsol
  */
 double BaseStar::CalculateLuminosityAtBAGB(double p_Mass) const {
@@ -1842,7 +1842,7 @@ void BaseStar::ResolveMassLoss() {
         m_COCoreMass=std::min(m_COCoreMass,m_Mass);                             // Not expected, only a precaution to avoid inconsistencies
         m_CoreMass=std::min(m_CoreMass, m_Mass);
 
-        UpdateInitialMass();                                                    // update initial mass (MS, HG & HeMS)  JR: todo: fix this kludge one day - mass0 is overloaded, and isn't always "initial mass"
+        UpdateInitialMass();                                                    // update effective initial mass (MS, HG & HeMS, HeHG)
         UpdateAgeAfterMassLoss();                                               // update age (MS, HG & HeMS)
         ApplyMassTransferRejuvenationFactor();                                  // apply age rejuvenation factor
     }
@@ -2964,14 +2964,10 @@ void BaseStar::UpdateAttributesAndAgeOneTimestepPreamble(const double p_DeltaMas
  *      attributes are updated.  The p_DeltaMass parameter may be zero, in which case no change is made to the
  *      star's mass before the attributes of the star are calculated.
  *
- *    - if required, the star's initial mass is changed by the amount passed as the p_DeltaMass0 parameter before
- *      other attributes are updated.  The p_DeltaMass parameter may be zero, in which case no change is made to
- *      the star's mass before the attributes of the star are calculated.  This should be used infrequently, and
- *      is really a kludge because the Mass0 attribute in Hurley et al. 2000 was overloaded after the introduction
- *      of mass loss (see section 7.1).  We should really separate the different uses of Mass0 in the code and
- *      use a different variable - initial mass shouldn't change (other than to initial mass upon entering a
- *      stellar phase - it doesn't make a lot of sense for initial mass to change during evolution through the
- *      phase).         JR: todo: action this paragraph.
+ *    - if required, the star's effective initial mass is changed by the amount passed as the p_DeltaMass0 parameter before
+ *      other attributes are updated.  The p_DeltaMass0 parameter may be zero, in which case no change is made to
+ *      the star's mass before the attributes of the star are calculated.  The Mass0 attribute in Hurley et al. 2000 is overloaded by the introduction
+ *      of mass loss (see section 7.1).
  *
  *    - if required, the star is aged by the amount passed as the p_DeltaTime parameter, and the simulation time is
  *      advanced by the same amount, before other attributes are updated.  The p_deltaTime parameter may be zero,

src/EAGB.cpp

@@ -82,7 +82,7 @@ double EAGB::CalculateLambdaNanjing() const {
 	DBL_VECTOR b        = {};                                                       // 0..5 b_coefficients
 
     if (utils::Compare(m_Metallicity, LAMBDA_NANJING_ZLIMIT) > 0) {                 // Z>0.5 Zsun: popI
-        if (utils::Compare(m_MZAMS, 1.5) < 0) {
+        if (utils::Compare(m_MZAMS, 1.5) < 0) {                                     // Should probably use effective mass m_Mass0 instead for Lambda calculations
             maxBG = { 2.5, 1.5 };
             if (utils::Compare(m_Radius, 200.0) > 0) lambdaBG = { 0.05, 0.05 };
             else  {
@@ -408,7 +408,7 @@ double EAGB::CalculateLambdaNanjing() const {
 /*
  * Calculate luminosity on the Early Asymptotic Giant Branch
  *
- * Hurley et al. 2000, eqs 61, 62 & 63
+ * Hurley et al. 2000, eq 37
  *
  *
  * double CalculateLuminosityOnPhase(const double p_CoreMass)

src/FGB.cpp

@@ -165,7 +165,7 @@ STELLAR_TYPE FGB::ResolveEnvelopeLoss(bool p_NoCheck) {
 
     if (p_NoCheck || utils::Compare(m_CoreMass, m_Mass) > 0) {                                      // Envelope loss
 
-        if (utils::Compare(m_Mass, massCutoffs(MHeF)) < 0) {                                        // Star evolves to Helium White Dwarf
+        if (utils::Compare(m_Mass0, massCutoffs(MHeF)) < 0) {                                        // Star evolves to Helium White Dwarf
 
             stellarType  = STELLAR_TYPE::HELIUM_WHITE_DWARF;
 
@@ -200,13 +200,10 @@ STELLAR_TYPE FGB::ResolveEnvelopeLoss(bool p_NoCheck) {
 /*
  * Modify the star due to (possible) helium flash
  *
- * JR: Is this described in Hurley et al. 2000?
- *
  *
  * void ResolveHeliumFlash()
  *
- * Deletermine if Helium Flash occurs, and if so modify the star accordingly.
- * The attributes of the star are updated.
+ * Deletermine if Helium Flash occurs, and if so set m_Mass0 equal to current mass as described in Hurley+ (2000), last paragraph before start of 7.1.1.
  */
 void FGB::ResolveHeliumFlash() {
 #define timescales(x) m_Timescales[static_cast<int>(TIMESCALE::x)]      // for convenience and readability - undefined at end of function

src/GiantBranch.cpp

@@ -530,7 +530,7 @@ double GiantBranch::CalculateLuminosityAtHeIgnition_Static(const double      p_M
 double GiantBranch::CalculateRemnantLuminosity() const {
 #define massCutoffs(x) m_MassCutoffs[static_cast<int>(MASS_CUTOFF::x)]  // for convenience and readability - undefined at end of function
 
-    return (utils::Compare(m_Mass, massCutoffs(MHeF)) < 0)
+    return (utils::Compare(m_Mass0, massCutoffs(MHeF)) < 0)
             ? HeMS::CalculateLuminosityAtZAMS_Static(m_CoreMass)
             : WhiteDwarfs::CalculateLuminosityOnPhase_Static(m_CoreMass, 0.0, m_Metallicity, WD_Baryon_Number.at(STELLAR_TYPE::HELIUM_WHITE_DWARF));
 
@@ -670,7 +670,7 @@ double GiantBranch::CalculateRadiusAtHeIgnition(const double p_Mass) const {
 double GiantBranch::CalculateRemnantRadius() const {
 #define massCutoffs(x) m_MassCutoffs[static_cast<int>(MASS_CUTOFF::x)]  // for convenience and readability - undefined at end of function
 
-    return (utils::Compare(m_Mass, massCutoffs(MHeF)) < 0)
+    return (utils::Compare(m_Mass0, massCutoffs(MHeF)) < 0)
             ? HeMS::CalculateRadiusAtZAMS_Static(m_CoreMass)
             : WhiteDwarfs::CalculateRadiusOnPhase_Static(m_CoreMass);
 

src/HG.cpp

@@ -594,11 +594,8 @@ double HG::CalculateCoreMassOnPhase(const double p_Mass, const double p_Time) co
 /*
  * Calculate rejuvenation factor for stellar age based on mass lost/gained during mass transfer
  *
- * JR: Description?
- *
- * Always returns 1.0 for HG - the rejuvenation factor is unity for convective main sequence stars.
- * The rest of the code is here so sanity checks can be made and an error displayed if a bad prescription
- * was specified in the program options
+ * Always returns 1.0 for HG
+ * The rest of the code is here so sanity checks can be made and an error displayed if a bad prescription was specified in the program options
  *
  *
  * double CalculateMassTransferRejuvenationFactor()
@@ -904,9 +901,23 @@ STELLAR_TYPE HG::EvolveToNextPhase() {
     }
     else {
         stellarType = STELLAR_TYPE::CORE_HELIUM_BURNING;
-    }
-
+    }    
     return stellarType;
 
 #undef massCutoffs
 }
+
+/*
+ * Update effective initial mass
+ *
+ * Per Hurley et al. 2000, section 7.1, the effective initial mass on the HG tracks the stellar mass -- unless it would yield an unphysical decrease in the core mass
+ *
+ *
+ * void UpdateInitialMass()
+ *
+ */
+void HG::UpdateInitialMass() {
+    if (utils::Compare(m_CoreMass,HG::CalculateCoreMassOnPhase(m_Mass, m_Age)) < 0 ) {        //The current mass would yield a core mass larger than the current core mass -- i.e., no unphysical core mass decrease would ensure
+        m_Mass0 = m_Mass;
+    }
+}

src/HG.h

@@ -108,7 +108,7 @@ class HG: virtual public BaseStar, public GiantBranch {
 
     void            UpdateAgeAfterMassLoss();                                                                                                                                   // Per Hurley et al. 2000, section 7.1
 
-    void            UpdateInitialMass()                                             { if (m_Mass0 > m_CoreMass) m_Mass0 = m_Mass; }                                             // Per Hurley et al. 2000, section 7.1
+    void            UpdateInitialMass();                                                   // Per Hurley et al. 2000, section 7.1
 
 };
 

src/MainSequence.cpp

@@ -57,7 +57,7 @@ double MainSequence::CalculateDeltaL(const double p_Mass) const {
 
     double deltaL;
 
-    if (utils::Compare(p_Mass, massCutoffs(MHook)) <= 0) {              // JR: todo: added "=" per Hurley et al. 2000, eq 16
+    if (utils::Compare(p_Mass, massCutoffs(MHook)) <= 0) {              // per Hurley et al. 2000, eq 16
         deltaL = 0.0;                                                   // this really is supposed to be zero
     }
     else if (utils::Compare(p_Mass, a[33]) < 0) {

src/changelog.h

@@ -848,7 +848,11 @@
 // 02.25.10     JR - Nov 19, 2021    - Defect repairs:
 //                                      - clamp timestep returned in BaseStar::CalculateTimestep() to NUCLEAR_MINIMUM_TIMESTEP
 //                                      - change NUCLEAR_MINIMUM_TIMESTEP to 1 year (from 100 years) in constants.h
+// 02.26.00     IM - Nov 30, 2021    - Defect repairs:
+//                                      - only decrease effective initial mass for HG and HeHG stars on mass loss when this decrease would not drive an unphysical decrease in the core mass
+//                                      - change mass comparisons (e.g., mass vs. He flash mass threshold) to compare effective initial mass rather than current mass
+//                                      - minor code and comment cleanup
 
-const std::string VERSION_STRING = "02.25.10";
+const std::string VERSION_STRING = "02.26.00";
 
 # endif // __changelog_h__