Skip to content

Commit

Permalink
Merge pull request #1102 from TeamCOMPAS/issue-978
Browse files Browse the repository at this point in the history 
Fixes for issue 978 and 1075.
  • Loading branch information
@jeffriley
jeffriley committed May 3, 2024
2 parents c9650e5 + 2e31e48 commit 27221c4
Show file tree
Hide file tree
Showing 25 changed files with 4,598 additions and 538 deletions.
443 changes: 219 additions & 224 deletions src/BaseStar.cpp
View file

Large diffs are not rendered by default.

12 changes: 6 additions & 6 deletions src/BaseStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -251,7 +251,7 @@ class BaseStar {

double CalculateRadialExpansionTimescale() const { return CalculateRadialExpansionTimescale_Static(m_StellarType, m_StellarTypePrev, m_Radius, m_RadiusPrev, m_DtPrev); } // Use class member variables

virtual double CalculateRadialExtentConvectiveEnvelope() const { return 0.0; } // default for stars with no convective envelope
virtual double CalculateRadialExtentConvectiveEnvelope() const { return 0.0; } // default for stars with no convective envelope

void CalculateSNAnomalies(const double p_Eccentricity);

Expand Down Expand Up @@ -376,10 +376,10 @@ class BaseStar {
double m_TZAMS0; // Effective ZAMS Temperature

// Current timestep variables
double m_Age; // Current effective age (changes with mass loss/gain)(myrs)
double m_Age; // Current effective age (changes with mass loss/gain) (Myr)
double m_COCoreMass; // Current CO core mass (Msol)
double m_CoreMass; // Current core mass (Msol)
double m_Dt; // Current timestep (myrs)
double m_Dt; // Size of current timestep (Myr)
bool m_EnvelopeJustExpelledByPulsations; // Flag to know if the convective envelope has just been expelled by pulsations
double m_HeCoreMass; // Current He core mass (Msol)
bool m_LBVphaseFlag; // Flag to know if the star satisfied the conditions, at any point in its evolution, to be considered a Luminous Blue Variable (LBV)
Expand All @@ -395,7 +395,7 @@ class BaseStar {
double m_Radius; // Current radius (Rsol)
double m_Tau; // Relative time
double m_Temperature; // Current temperature (Tsol)
double m_Time; // Current physical time the star has been evolved(myrs)
double m_Time; // Current physical time the star has been evolved (Myr)

// Previous timestep variables
double m_DtPrev; // Previous timestep
Expand Down Expand Up @@ -526,7 +526,7 @@ class BaseStar {

static double CalculateMassLoss_Static(const double p_Mass, const double p_Mdot, const double p_Dt);

double CalculateMassLossRate();
virtual double CalculateMassLossRate();
virtual double CalculateMassLossRateHurley();
double CalculateMassLossRateKudritzkiReimers() const;
double CalculateMassLossRateLBV(const LBV_PRESCRIPTION p_LBV_prescription);
Expand Down Expand Up @@ -557,7 +557,7 @@ class BaseStar {
double CalculateMassLossRateWolfRayetSanderVink2020(const double p_Mu) const;
double CalculateMassLossRateWolfRayetTemperatureCorrectionSander2023(const double p_Mdot) const;
double CalculateMassLossRateHeliumStarVink2017() const;
double CalculateMassLossRateWolfRayetShenar2019() const;
virtual double CalculateMassLossRateWolfRayetShenar2019() const;

virtual double CalculateMassTransferRejuvenationFactor() const;

Expand Down
48 changes: 30 additions & 18 deletions src/CHeB.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ void CHeB::CalculateTimescales(const double p_Mass, DBL_VECTOR &p_Timescales) {
// particularly eq 58 and beyond). COMPAS does not allow for a 0-length blue loop - some of
// the equations used (e.g. to calculate Radius) result in nan or inf. As a temporary workaround
// until we work out how to skip the blue loop (when it is 0-length) we will set the length of a
// 0-length blue loop to the absolute minimum timestep (currently 100 seconds).
// 0-length blue loop to the absolute minimum timestep (currently ~100 seconds).
//
// Note that this works around a long-standing problem, which was worked around in legacy COMPAS
// by the following code in calculateBluePhaseFBL() in star.cpp:
Expand Down Expand Up @@ -931,17 +931,24 @@ double CHeB::CalculateLuminosityOnPhase(const double p_Mass, const double p_Tau)
double tx = timescales(tauX_BL); // 0 for LM and HM stars, non-zero for IM stars
double Lx = CalculateLuminosityAtBluePhaseStart(p_Mass);

if (utils::Compare(p_Tau, tx) >= 0) {
if (utils::Compare(p_Tau, tx) >= 0) { // on the blue loop
double Rx = CalculateRadiusAtBluePhaseStart(p_Mass);
double RmHe = CalculateMinimumRadiusOnPhase_Static(p_Mass, m_CoreMass, m_Alpha1, massCutoffs(MHeF), massCutoffs(MFGB), m_MinimumLuminosityOnPhase, m_BnCoefficients);
double LBAGB = CalculateLuminosityAtBAGB(p_Mass);

// the following check for high mass stars was added to match the Hurley sse code
// - see Hurley sse `hrdiag.f` line 297
if (p_Mass > HIGH_MASS_THRESHOLD) {
Rx = RmHe;
}

double epsilon = std::min(2.5, std::max(0.4, (RmHe / Rx)));
double lambda = (utils::Compare(p_Tau, tx) == 0) ? 0.0 : PPOW((p_Tau - tx) / (1.0 - tx), epsilon); // JR: tx can be 1.0 here - if so, lambda = 0.0
double lambda = (utils::Compare(p_Tau, tx) == 0) ? 0.0 : PPOW((p_Tau - tx) / (1.0 - tx), epsilon); // tx can be 1.0 here - if so, lambda = 0.0
lCHeB = Lx * PPOW(LBAGB / Lx, lambda);
}
else {
else { // before the blue loop
double LHeI = GiantBranch::CalculateLuminosityAtHeIgnition_Static(p_Mass, m_Alpha1, massCutoffs(MHeF), m_BnCoefficients); // pow() is slow - use multiplication
double tmp = (tx - p_Tau) / tx; // JR: tx cannot be 0.0 here, so safe (tx > tau, tau = [0, 1])
double tmp = (tx - p_Tau) / tx; // tx cannot be 0.0 here, so safe (tx > tau, tau = [0, 1])
double lambdaPrime = tmp * tmp * tmp;
lCHeB = Lx * PPOW((LHeI / Lx), lambdaPrime);
}
Expand Down Expand Up @@ -1016,17 +1023,17 @@ double CHeB::CalculateMinimumRadiusOnPhase_Static(const double p_Mass,
const DBL_VECTOR &p_BnCoefficients) {
#define b p_BnCoefficients // for convenience and readability - undefined at end of function

double minR = 0.0; // Minimum radius
double minR = 0.0; // Minimum radius

if (utils::Compare(p_MHeF, p_Mass) < 0) {
double m_b28 = PPOW(p_Mass, b[28]); // pow() is slow - do it once only
double m_b28 = PPOW(p_Mass, b[28]); // pow() is slow - do it once only
minR = ((b[24] * p_Mass) + (PPOW((b[25] * p_Mass), b[26]) * m_b28)) / (b[27] + m_b28);
}
else {
double LZAHB_MHeF = GiantBranch::CalculateLuminosityOnZAHB_Static(p_MHeF, p_CoreMass, p_Alpha1, p_MHeF, p_MFGB, p_MinimumLuminosityOnPhase, p_BnCoefficients);
double LZAHB = GiantBranch::CalculateLuminosityOnZAHB_Static(p_Mass, p_CoreMass, p_Alpha1, p_MHeF, p_MFGB, p_MinimumLuminosityOnPhase, p_BnCoefficients);
double mu = p_Mass / p_MHeF;
double MHeF_b28 = PPOW(p_MHeF, b[28]); // pow() is slow - do it once only
double MHeF_b28 = PPOW(p_MHeF, b[28]); // pow() is slow - do it once only
double top = ((b[24] * p_MHeF) + (PPOW((b[25] * p_MHeF), b[26]) * MHeF_b28)) / (b[27] + MHeF_b28);
double bottom = GiantBranch::CalculateRadiusOnPhase_Static(p_MHeF, LZAHB_MHeF, p_BnCoefficients);

Expand Down Expand Up @@ -1253,7 +1260,7 @@ double CHeB::CalculateTauOnPhase() const {
* double CalculateLifetimeOnPhase(const double p_Mass)
*
* @param [IN] p_Mass Mass in Msol
* @return Lifetime of Core Helium Burning in MYRs (tHe)
* @return Lifetime of Core Helium Burning in Myr (tHe)
*
* JR: changed this to use m_Timescales[TS::tBGB] instead of parameter
*/
Expand Down Expand Up @@ -1305,12 +1312,12 @@ double CHeB::CalculateBluePhaseFBL(const double p_Mass) {

// Might be that we are supposed to use min(RmHe, Rx=RHeI)
double RHeI = CalculateRadiusAtHeIgnition(p_Mass);
double LHeI = GiantBranch::CalculateLuminosityAtHeIgnition_Static(p_Mass, m_Alpha1, massCutoffs(MHeF), m_BnCoefficients);
double LHeI = GiantBranch::CalculateLuminosityAtHeIgnition_Static(p_Mass, m_Alpha1, massCutoffs(MHeF), b);

top = std::min(top, RHeI);

// Calculate RAGB(LHeI(M)) for M > MFGB > MHeF
double bottom = EAGB::CalculateRadiusOnPhase_Static(p_Mass, LHeI, massCutoffs(MHeF), m_BnCoefficients);
double bottom = EAGB::CalculateRadiusOnPhase_Static(p_Mass, LHeI, massCutoffs(MHeF), b);
double brackets = 1.0 - (top / bottom);

return PPOW(p_Mass, b[48]) * PPOW(brackets, b[49]);
Expand All @@ -1331,7 +1338,7 @@ double CHeB::CalculateBluePhaseFBL(const double p_Mass) {
* @param [IN] p_Mass Mass in Msol
* @return Relative lifetime of blue phase of Core Helium Burning, clamped to [0, 1]
*/
double CHeB::CalculateLifetimeOnBluePhase(const double p_Mass) {
double CHeB::CalculateLifetimeOnBluePhase(const double p_Mass) {
#define b m_BnCoefficients // for convenience and readability - undefined at end of function
#define massCutoffs(x) m_MassCutoffs[static_cast<int>(MASS_CUTOFF::x)] // for convenience and readability - undefined at end of function

Expand All @@ -1341,9 +1348,9 @@ double CHeB::CalculateLifetimeOnBluePhase(const double p_Mass) {
tbl = 1.0;
}
else if (utils::Compare(p_Mass, massCutoffs(MFGB)) <= 0) {
double m_MFGB = p_Mass / massCutoffs(MFGB);
double firstTerm = (b[45] * PPOW(m_MFGB, 0.414));
tbl = firstTerm + ((1.0 - firstTerm) * PPOW((log10(m_MFGB) / log10(massCutoffs(MHeF) / massCutoffs(MFGB))), b[46]));
double mass_MFGB = p_Mass / massCutoffs(MFGB);
double firstTerm = (b[45] * PPOW(mass_MFGB, 0.414));
tbl = firstTerm + ((1.0 - firstTerm) * PPOW((log10(mass_MFGB) / log10(massCutoffs(MHeF) / massCutoffs(MFGB))), b[46]));
}
else {
double fblM = CalculateBluePhaseFBL(p_Mass);
Expand All @@ -1368,11 +1375,16 @@ double CHeB::CalculateLifetimeOnBluePhase(const double p_Mass) {
* @return true if evolution should continue on phase, false otherwise
*/
bool CHeB::ShouldEvolveOnPhase() const {
bool afterHeIgnition = (m_Age >= m_Timescales[static_cast<int>(TIMESCALE::tHeI)]);
bool beforeEndOfHeBurning = (m_Age < (m_Timescales[static_cast<int>(TIMESCALE::tHeI)] + m_Timescales[static_cast<int>(TIMESCALE::tHe)]));
bool coreIsNotTooMassive = (m_HeCoreMass < m_Mass);
#define timescales(x) m_Timescales[static_cast<int>(TIMESCALE::x)] // for convenience and readability - undefined at end of function

bool afterHeIgnition = m_Age >= timescales(tHeI);
bool beforeEndOfHeBurning = m_Age < (timescales(tHeI) + timescales(tHe));
bool coreIsNotTooMassive = utils::Compare(m_HeCoreMass, m_Mass) < 0;

// Evolve on CHeB phase if age after He Ign and while He Burning and He core mass does not exceed total mass (could happen due to mass loss)
return (afterHeIgnition && beforeEndOfHeBurning && coreIsNotTooMassive && !ShouldEnvelopeBeExpelledByPulsations());

#undef timescales
}


Expand Down
14 changes: 7 additions & 7 deletions src/EAGB.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -383,7 +383,7 @@ double EAGB::CalculateLambdaNanjingStarTrack(const double p_Mass, const double p
DBL_VECTOR b = {}; // 0..5 b_coefficients

if (utils::Compare(p_Metallicity, LAMBDA_NANJING_ZLIMIT) > 0) { // Z>0.5 Zsun: popI
if (utils::Compare(p_Mass, 1.5) < 0) { // Should probably use effective mass m_Mass0 instead for Lambda calculations
if (utils::Compare(p_Mass, 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 {
Expand Down Expand Up @@ -773,9 +773,8 @@ double EAGB::CalculateRadiusOnPhase_Static(const double p_Mass,
#define b p_BnCoefficients // for convenience and readability - undefined at end of function

// sanity check for mass and luminosity - just return 0.0 if mass or luminosity <= 0
// doing this will save some compute cycles - but it is not strictly what the equation in Hurley says
// (Hurley et al. 2000, eq 74 and immediately following) - there if mass is 0.0 but luminosity is
// non-zero we get a non-zero value for radius (shouldn't happen, but we need to code for all possibilities)
// doing this will save some compute cycles

if (utils::Compare(p_Mass, 0.0) <= 0 || utils::Compare(p_Luminosity, 0.0) <= 0) return 0.0;

// calculate radius constant A (Hurley et al. 2000, eq 74)
Expand Down Expand Up @@ -855,8 +854,8 @@ double EAGB::CalculateRemnantRadius() const {
* @return Core mass on the Early Asymptotic Giant Branch in Msol
*/
double EAGB::CalculateCOCoreMassOnPhase(const double p_Time) const {
#define timescales(x) m_Timescales[static_cast<int>(TIMESCALE::x)] // for convenience and readability - undefined at end of function
#define gbParams(x) m_GBParams[static_cast<int>(GBP::x)] // for convenience and readability - undefined at end of function
#define timescales(x) m_Timescales[static_cast<int>(TIMESCALE::x)] // for convenience and readability - undefined at end of function
#define gbParams(x) m_GBParams[static_cast<int>(GBP::x)] // for convenience and readability - undefined at end of function

return utils::Compare(p_Time, timescales(tMx_FAGB)) <= 0
? PPOW((gbParams(p) - 1.0) * gbParams(AHe) * gbParams(D) * (timescales(tinf1_FAGB) - p_Time), 1.0 / (1.0 - gbParams(p)))
Expand Down Expand Up @@ -1027,6 +1026,7 @@ STELLAR_TYPE EAGB::ResolveEnvelopeLoss(bool p_NoCheck) {
timescales(tinf2_HeGB) = timescales(tx_HeGB) + ((1.0 / (q1 * gbParams(AHe) * gbParams(B))) * PPOW((gbParams(B) / gbParams(Lx)), q1_q));

m_Age = HeGB::CalculateAgeOnPhase_Static(m_Mass, m_COCoreMass, timescales(tHeMS), m_GBParams);

HeHG::CalculateGBParams_Static(m_Mass0, m_Mass, LogMetallicityXi(), m_MassCutoffs, m_AnCoefficients, m_BnCoefficients, m_GBParams); // IM: order of type change and parameter updates to be revisited (e.g., why not just CalculateGBParams(m_Mass0, m_GBParams)?) JR: static function has no access to class variables
m_Luminosity = HeGB::CalculateLuminosityOnPhase_Static(m_COCoreMass, gbParams(B), gbParams(D));

Expand Down Expand Up @@ -1107,7 +1107,7 @@ bool EAGB::IsSupernova() const {
double McCOBAGB = CalculateCOCoreMassOnPhase(timescales(tHeI) + timescales(tHe));
double McSN = std::max(gbParams(McSN), 1.05 * McCOBAGB); // hack from Hurley fortran code, doesn't seem to be in the paper JR: do we know why?

return (utils::Compare(McSN, m_COCoreMass) < 0); // core is heavy enough to go Supernova
return (utils::Compare(McSN, m_COCoreMass) <= 0); // core is heavy enough to go Supernova

#undef gbParams
#undef timescales
Expand Down
2 changes: 1 addition & 1 deletion src/FGB.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ double FGB::CalculateLuminosityOnPhase(const double p_Time) const {
* double CalculateCoreMassOnPhase(const double p_Mass, const double p_Time)
*
* @param [IN] p_Mass Mass in Msol
* @param [IN] p_Time Time after ZAMS in MYRS (tBGB <= time <= tHeI)
* @param [IN] p_Time Time after ZAMS in Myr (tBGB <= time <= tHeI)
* @return Core mass on the First Giant Branch in Msol
*/
double FGB::CalculateCoreMassOnPhase(const double p_Mass, const double p_Time) const {
Expand Down

0 comments on commit 27221c4

Please sign in to comment.