From de2bf7a544422e0222675c518af5a0e61679f37d Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Sat, 30 Nov 2024 14:53:35 +1100 Subject: [PATCH 1/7] Bug Fixes --- src/BaseBinaryStar.cpp | 12 +++++++++++- src/BaseStar.cpp | 4 +--- src/HeHG.cpp | 8 +++++++- 3 files changed, 19 insertions(+), 5 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 951daddd0..3531540e0 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -2086,6 +2086,8 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { // no m_MassTransferTrackerHistory = m_Donor == m_Star1 ? MT_TRACKING::STABLE_1_TO_2_SURV : MT_TRACKING::STABLE_2_TO_1_SURV; // record what happened - for later printing double massGainAccretor = -massDiffDonor * m_FractionAccreted; // set accretor mass gain to mass loss * conservativeness + + double omegaDonor_pre_MT = m_Donor->Omega(); // used if full donor envelope is removed m_Donor->SetMassTransferDiffAndResolveWDShellChange(massDiffDonor); // set new mass of donor m_Accretor->SetMassTransferDiffAndResolveWDShellChange(massGainAccretor); // set new mass of accretor @@ -2095,7 +2097,12 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { m_aMassTransferDiff = aFinal - aInitial; // set change in orbit (semi-major axis) STELLAR_TYPE stellarTypeDonor = m_Donor->StellarType(); // donor stellar type before resolving envelope loss - if (isEnvelopeRemoved) m_Donor->ResolveEnvelopeLossAndSwitch(); // if this was an envelope stripping episode, resolve envelope loss + + if (isEnvelopeRemoved) { // if this was an envelope stripping episode, resolve envelope loss + m_Star1->ResolveEnvelopeLossAndSwitch(); // resolve envelope loss for the donor and switch to new stellar type + m_Star1->SetOmega(omegaDonor_pre_MT); // keep the rotation frequency of the core equal to the pre-envelope-loss rotation frequency + } + if (m_Donor->StellarType() != stellarTypeDonor) { // stellar type change? (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::STELLAR_TYPE_CHANGE_DURING_MT); // yes - print (log) detailed output } @@ -2430,6 +2437,8 @@ void BaseBinaryStar::ResolveMassChanges() { STELLAR_TYPE stellarType1 = m_Star1->StellarTypePrev(); // star 1 stellar type before updating attributes STELLAR_TYPE stellarType2 = m_Star2->StellarTypePrev(); // star 2 stellar type before updating attributes + //std::cout<<"Time"<MassPrev()<<"cur"<Mass()<<"change"<MassLossDiff() + m_Star1->MassTransferDiff()<<"loss"<MassLossDiff()<<"transfer"<MassTransferDiff()<MassPrev(), m_Star1->Mass()) == 0) { // mass already updated? // no - resolve mass changes double massChange = m_Star1->MassLossDiff() + m_Star1->MassTransferDiff(); // mass change due to winds and mass transfer + if (utils::Compare(massChange, 0.0) != 0) { // winds/mass transfer changes mass? // yes - calculate new angular momentum; assume accretor is adding angular momentum from a circular orbit at the stellar radius double angularMomentumChange = (utils::Compare(massChange, 0.0) > 0) ? diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index 500fb76d3..fa0273cef 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -2494,8 +2494,6 @@ double BaseStar::CalculateMassLossRateWolfRayetTemperatureCorrectionSander2023(c else { logMdotCorrected = logMdotUncorrected; } - - std::cout<(GBP::x)] // for convenience and readability - undefined at end of function - + GiantBranch::CalculateGBParams_Static(p_Mass, p_LogMetallicityXi, p_MassCutoffs, p_AnCoefficients, p_BnCoefficients, p_GBParams); // calculate common values (actually, all) // recalculate HeHG specific values @@ -125,6 +125,12 @@ void HeHG::CalculateGBParams_Static(const double p_Mass0, gbParams(Mx) = GiantBranch::CalculateCoreMass_Luminosity_Mx_Static(p_GBParams); // depends on B, D, p & q - recalculate if any of those are changed gbParams(Lx) = GiantBranch::CalculateCoreMass_Luminosity_Lx_Static(p_GBParams); // depends on B, D, p, q & Mx - recalculate if any of those are changed + // need to reset p and q to HeHG specific versions -- while they are set in GiantBranch::CalculateGBParams_Static(), that uses + // the GiantBranch versions of CalculateCoreMass_Luminosity_p_Static and CalculateCoreMass_Luminosity_q_Static + // should return to this and understand desired behavior... + gbParams(p) = CalculateCoreMass_Luminosity_p_Static(p_Mass, p_MassCutoffs); + gbParams(q) = CalculateCoreMass_Luminosity_q_Static(p_Mass, p_MassCutoffs); + gbParams(McBAGB) = p_Mass0; gbParams(McBGB) = GiantBranch::CalculateCoreMassAtBGB_Static(p_Mass, p_MassCutoffs, p_AnCoefficients, p_GBParams); From d8c19cf4dcf6cb9a3fcacbf850b615511eed15e2 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Sat, 30 Nov 2024 14:57:40 +1100 Subject: [PATCH 2/7] Draft PR --- src/changelog.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/changelog.h b/src/changelog.h index 7aedd9410..2fb4cd86d 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1411,6 +1411,9 @@ // - Minor fixes, including to #1255, #1258 // 03.10.00 JR - Nov 29, 2024 - Enhancement: // - added functionality to allow stellar mergers (for BSE) to be logged to switchlog file (see documentation for details) -const std::string VERSION_STRING = "03.10.00"; +// 03.10.01 IM - Nov 30, 2024 - Defect repair: +// - corrected treatment of rotation to retain pre-mass-loss spin frequency, not angular momentum, on complete envelope removal during stable mass transfer +// - fixed issue with updating helium giants that manifested as supernovae with nan core mass (see #1245) +const std::string VERSION_STRING = "03.10.01"; # endif // __changelog_h__ From 6d755f974b2f7aa46593cd62faf866e473934bfb Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Sat, 30 Nov 2024 16:09:07 +1100 Subject: [PATCH 3/7] added check for exceeding Chandrasekhar mass when computing white dwarf radius (resolves issue #1264) --- src/WhiteDwarfs.cpp | 2 ++ src/changelog.h | 1 + 2 files changed, 3 insertions(+) diff --git a/src/WhiteDwarfs.cpp b/src/WhiteDwarfs.cpp index 3e49e3b07..2444958bb 100644 --- a/src/WhiteDwarfs.cpp +++ b/src/WhiteDwarfs.cpp @@ -155,6 +155,8 @@ double WhiteDwarfs::CalculateRadiusOnPhase_Static(const double p_Mass) { // sanity check for mass - just return 0.0 if mass <= 0 if (utils::Compare(p_Mass, 0.0) <= 0) return 0.0; + if (utils::Compare(p_Mass, MCH) >= 0) return NEUTRON_STAR_RADIUS; // only likely to come up if asking for the core or remnant radius of a giant star + double MCH_Mass_one_third = std::cbrt(MCH / p_Mass); double MCH_Mass_two_thirds = MCH_Mass_one_third * MCH_Mass_one_third; diff --git a/src/changelog.h b/src/changelog.h index 2fb4cd86d..8a0f539c3 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1414,6 +1414,7 @@ // 03.10.01 IM - Nov 30, 2024 - Defect repair: // - corrected treatment of rotation to retain pre-mass-loss spin frequency, not angular momentum, on complete envelope removal during stable mass transfer // - fixed issue with updating helium giants that manifested as supernovae with nan core mass (see #1245) +// - added check for exceeding Chandrasekhar mass when computing white dwarf radius (resolves issue #1264) const std::string VERSION_STRING = "03.10.01"; # endif // __changelog_h__ From fd798aa02b2554daab16c3cc013cab88ca98f6fc Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Sat, 30 Nov 2024 16:35:09 +1100 Subject: [PATCH 4/7] added check to only compute McBGB for stars with mass above MHeF, following text above Eq. 44 in Hurley+, 2000 (resolves issue #1256) --- src/GiantBranch.cpp | 4 +++- src/changelog.h | 1 + 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/GiantBranch.cpp b/src/GiantBranch.cpp index 1602a0b31..fb48a2e0b 100644 --- a/src/GiantBranch.cpp +++ b/src/GiantBranch.cpp @@ -770,10 +770,12 @@ double GiantBranch::CalculateCoreMassAtBGB(const double p_Mass, const DBL_VECTOR #define gbParams(x) p_GBParams[static_cast(GBP::x)] // for convenience and readability - undefined at end of function #define massCutoffs(x) m_MassCutoffs[static_cast(MASS_CUTOFF::x)] // for convenience and readability - undefined at end of function + if (utils::Compare(p_Mass, massCutoffs(MHeF)) <= 0) return 0.0; // No McBGB for stars with mass below the helium flash threshold, see text above Eq. (44) of Hurley+ (2000) + double luminosity = GiantBranch::CalculateLuminosityAtPhaseBase_Static(massCutoffs(MHeF), m_AnCoefficients); double Mc_MHeF = BaseStar::CalculateCoreMassGivenLuminosity_Static(luminosity, p_GBParams); double c = (Mc_MHeF * Mc_MHeF * Mc_MHeF * Mc_MHeF) - (MC_L_C1 * PPOW(massCutoffs(MHeF), MC_L_C2)); // pow() is slow - use multiplication - + return std::min((0.95 * gbParams(McBAGB)), std::sqrt(std::sqrt(c + (MC_L_C1 * PPOW(p_Mass, MC_L_C2))))); // sqrt is much faster than PPOW() #undef massCutoffs diff --git a/src/changelog.h b/src/changelog.h index 8a0f539c3..f6e40f14f 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1415,6 +1415,7 @@ // - corrected treatment of rotation to retain pre-mass-loss spin frequency, not angular momentum, on complete envelope removal during stable mass transfer // - fixed issue with updating helium giants that manifested as supernovae with nan core mass (see #1245) // - added check for exceeding Chandrasekhar mass when computing white dwarf radius (resolves issue #1264) +// - added check to only compute McBGB for stars with mass above MHeF, following text above Eq. 44 in Hurley+, 2000 (resolves issue #1256) const std::string VERSION_STRING = "03.10.01"; # endif // __changelog_h__ From fc8dbffb66fa41e404bca8f737d2f7607c90ba5a Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Sat, 30 Nov 2024 16:38:10 +1100 Subject: [PATCH 5/7] Typo fix --- src/BaseBinaryStar.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index cf215befc..1d0fde17a 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -2099,8 +2099,8 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { STELLAR_TYPE stellarTypeDonor = m_Donor->StellarType(); // donor stellar type before resolving envelope loss if (isEnvelopeRemoved) { // if this was an envelope stripping episode, resolve envelope loss - m_Star1->ResolveEnvelopeLossAndSwitch(); // resolve envelope loss for the donor and switch to new stellar type - m_Star1->SetOmega(omegaDonor_pre_MT); // keep the rotation frequency of the core equal to the pre-envelope-loss rotation frequency + m_Donor->ResolveEnvelopeLossAndSwitch(); // resolve envelope loss for the donor and switch to new stellar type + m_Donor->SetOmega(omegaDonor_pre_MT); // keep the rotation frequency of the core equal to the pre-envelope-loss rotation frequency } if (m_Donor->StellarType() != stellarTypeDonor) { // stellar type change? From 1be607daed6dd2e6cffbeb85df22d54872dfae18 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Sat, 30 Nov 2024 16:38:58 +1100 Subject: [PATCH 6/7] Cleaning --- src/BaseBinaryStar.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 1d0fde17a..624fa9f93 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -2437,8 +2437,6 @@ void BaseBinaryStar::ResolveMassChanges() { STELLAR_TYPE stellarType1 = m_Star1->StellarTypePrev(); // star 1 stellar type before updating attributes STELLAR_TYPE stellarType2 = m_Star2->StellarTypePrev(); // star 2 stellar type before updating attributes - //std::cout<<"Time"<MassPrev()<<"cur"<Mass()<<"change"<MassLossDiff() + m_Star1->MassTransferDiff()<<"loss"<MassLossDiff()<<"transfer"<MassTransferDiff()< Date: Sat, 30 Nov 2024 16:41:21 +1100 Subject: [PATCH 7/7] Comment cleanup --- src/HeHG.cpp | 2 +- src/WhiteDwarfs.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/HeHG.cpp b/src/HeHG.cpp index b6d434c03..59ecb63f8 100755 --- a/src/HeHG.cpp +++ b/src/HeHG.cpp @@ -127,7 +127,7 @@ void HeHG::CalculateGBParams_Static(const double p_Mass0, // need to reset p and q to HeHG specific versions -- while they are set in GiantBranch::CalculateGBParams_Static(), that uses // the GiantBranch versions of CalculateCoreMass_Luminosity_p_Static and CalculateCoreMass_Luminosity_q_Static - // should return to this and understand desired behavior... + // should return to this and understand desired behavior - *ILYA* gbParams(p) = CalculateCoreMass_Luminosity_p_Static(p_Mass, p_MassCutoffs); gbParams(q) = CalculateCoreMass_Luminosity_q_Static(p_Mass, p_MassCutoffs); diff --git a/src/WhiteDwarfs.cpp b/src/WhiteDwarfs.cpp index 2444958bb..494647678 100644 --- a/src/WhiteDwarfs.cpp +++ b/src/WhiteDwarfs.cpp @@ -155,7 +155,7 @@ double WhiteDwarfs::CalculateRadiusOnPhase_Static(const double p_Mass) { // sanity check for mass - just return 0.0 if mass <= 0 if (utils::Compare(p_Mass, 0.0) <= 0) return 0.0; - if (utils::Compare(p_Mass, MCH) >= 0) return NEUTRON_STAR_RADIUS; // only likely to come up if asking for the core or remnant radius of a giant star + if (utils::Compare(p_Mass, MCH) >= 0) return NEUTRON_STAR_RADIUS; // only expected to come up if asking for the core or remnant radius of a giant star double MCH_Mass_one_third = std::cbrt(MCH / p_Mass); double MCH_Mass_two_thirds = MCH_Mass_one_third * MCH_Mass_one_third;