From 77c30d8ec1a00844e56a23804ac81028426830d9 Mon Sep 17 00:00:00 2001 From: Veome Kapil Date: Wed, 13 Nov 2024 14:10:38 -0500 Subject: [PATCH 1/6] Added MomentOfInertiaPrevAU, and used it to ensure mass and radial changes preserve angular momentum by adjusting Omega accordingly --- src/BaseBinaryStar.cpp | 5 +++++ src/BaseStar.cpp | 2 ++ src/BaseStar.h | 2 ++ src/Star.h | 1 + src/changelog.h | 5 ++++- 5 files changed, 14 insertions(+), 1 deletion(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index a62569050..c07df3c09 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -2468,6 +2468,11 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) { // if m_Omega == 0.0 (should only happen on the first timestep), calculate m_Omega here if (utils::Compare(m_Omega, 0.0) == 0) m_Omega = OrbitalAngularVelocity(); + + // If the stars have changed since the previous timestep, adjust their spins to conserve angular momentum + m_Star1->SetOmega(m_Star1->Omega() * m_Star1->MomentOfInertiaPrevAU() / m_Star1->CalculateMomentOfInertiaAU()); + m_Star2->SetOmega(m_Star2->Omega() * m_Star2->MomentOfInertiaPrevAU() / m_Star2->CalculateMomentOfInertiaAU()); + m_TotalAngularMomentum = CalculateAngularMomentum(); } switch (OPTIONS->TidesPrescription()) { // which tides prescription? diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index 9b4586a7c..b83ed93ac 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -160,6 +160,7 @@ BaseStar::BaseStar(const unsigned long int p_RandomSeed, m_StellarTypePrev = m_StellarType; m_MassPrev = m_MZAMS; m_RadiusPrev = m_RZAMS; + m_MomentOfInertiaPrevAU = CalculateMomentOfInertiaAU(); m_DtPrev = DEFAULT_INITIAL_DOUBLE_VALUE; m_OmegaPrev = m_OmegaZAMS; @@ -4589,6 +4590,7 @@ void BaseStar::UpdateAttributesAndAgeOneTimestepPreamble(const double p_DeltaMas m_StellarTypePrev = m_StellarType; m_MassPrev = m_Mass; m_RadiusPrev = m_Radius; + m_MomentOfInertiaPrevAU = CalculateMomentOfInertiaAU(); } // the GBParams and Timescale calculations need to be done before taking the timestep - since diff --git a/src/BaseStar.h b/src/BaseStar.h index 94bd17095..7c44ae1f9 100644 --- a/src/BaseStar.h +++ b/src/BaseStar.h @@ -144,6 +144,7 @@ class BaseStar { double Mdot() const { return m_Mdot; } double Metallicity() const { return m_Metallicity; } double MinimumCoreMass() const { return m_MinimumCoreMass; } + double MomentOfInertiaPrevAU() const { return m_MomentOfInertiaPrevAU; } double MZAMS() const { return m_MZAMS; } double Omega() const { return m_Omega; } double OmegaCHE() const { return m_OmegaCHE; } @@ -446,6 +447,7 @@ class BaseStar { double m_MassPrev; // Previous mass (Msol) double m_OmegaPrev; // Previous angular frequency (yr^-1) double m_RadiusPrev; // Previous radius (Rsol) + double m_MomentOfInertiaPrevAU; // Previous moment of inertia (Msol AU^2) STELLAR_TYPE m_StellarTypePrev; // Stellar type at previous timestep // Metallicity variables diff --git a/src/Star.h b/src/Star.h index 3c3932c50..94bbce89f 100755 --- a/src/Star.h +++ b/src/Star.h @@ -121,6 +121,7 @@ class Star { double Mass0() const { return m_Star->Mass0(); } double MassPrev() const { return m_Star->MassPrev(); } double Metallicity() const { return m_Star->Metallicity(); } + double MomentOfInertiaPrevAU() const { return m_Star->MomentOfInertiaPrevAU();} double MZAMS() const { return m_Star->MZAMS(); } double Omega() const { return m_Star->Omega(); } double OmegaCHE() const { return m_Star->OmegaCHE(); } diff --git a/src/changelog.h b/src/changelog.h index 6e47847e7..18c005a59 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1375,7 +1375,10 @@ // - Addressed documentation issues in issues #1272 and #1273. // 03.07.05 IM - Nov 07, 2024 - Defect repairs: // - Fix for issue #1270: root finder functions now check if either of the bracket edges provides a sufficiently good solution, which sometimes happens when the initial guess is very close to the truth +// 03.07.06 VK - Nov 13, 2024 - Enhancement: +// - Added m_MomentOfInertiaPrevAU member variable to class BaseStar. +// - Tides behavior updated to adjust initial stellar spins, to account for any mass or radius changes since previous timestep while conserving angular momentum for each star. -const std::string VERSION_STRING = "03.07.05"; +const std::string VERSION_STRING = "03.07.06"; # endif // __changelog_h__ From cc89f2e5bd5d117a5ec695312f36598f36c89c08 Mon Sep 17 00:00:00 2001 From: Veome Kapil Date: Wed, 13 Nov 2024 16:24:44 -0500 Subject: [PATCH 2/6] added condition so that angular momentum is only conserved within the same stellar type --- src/BaseBinaryStar.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index c07df3c09..fcf11a043 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -2469,10 +2469,10 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) { // if m_Omega == 0.0 (should only happen on the first timestep), calculate m_Omega here if (utils::Compare(m_Omega, 0.0) == 0) m_Omega = OrbitalAngularVelocity(); - // If the stars have changed since the previous timestep, adjust their spins to conserve angular momentum - m_Star1->SetOmega(m_Star1->Omega() * m_Star1->MomentOfInertiaPrevAU() / m_Star1->CalculateMomentOfInertiaAU()); - m_Star2->SetOmega(m_Star2->Omega() * m_Star2->MomentOfInertiaPrevAU() / m_Star2->CalculateMomentOfInertiaAU()); - m_TotalAngularMomentum = CalculateAngularMomentum(); + // If the stars have changed since the previous timestep (and not changed stellar type), adjust their spins to conserve angular momentum + if (m_Star1->StellarType() == m_Star1->StellarTypePrev()) m_Star1->SetOmega(m_Star1->Omega() * m_Star1->MomentOfInertiaPrevAU() / m_Star1->CalculateMomentOfInertiaAU()); + if (m_Star2->StellarType() == m_Star2->StellarTypePrev()) m_Star2->SetOmega(m_Star2->Omega() * m_Star2->MomentOfInertiaPrevAU() / m_Star2->CalculateMomentOfInertiaAU()); + CalculateEnergyAndAngularMomentum(); } switch (OPTIONS->TidesPrescription()) { // which tides prescription? From 4a82b9db12c7657363812a431494e2d7962eb2a8 Mon Sep 17 00:00:00 2001 From: Veome Kapil Date: Wed, 20 Nov 2024 11:05:16 -0500 Subject: [PATCH 3/6] Fixed code for retaining core spin after envelope ejection --- src/BaseBinaryStar.cpp | 4 ++++ src/Star.cpp | 5 ----- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index c1cbc3f41..ab2bd8320 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1465,6 +1465,8 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { double periastronRsol = PeriastronRsol(); // periastron, Rsol (before CEE) double rRLd1Rsol = periastronRsol * CalculateRocheLobeRadius_Static(m_Star1->Mass(), m_Star2->Mass()); // Roche-lobe radius at periastron in Rsol at the moment where CEE begins, seen by star1 double rRLd2Rsol = periastronRsol * CalculateRocheLobeRadius_Static(m_Star2->Mass(), m_Star1->Mass()); // Roche-lobe radius at periastron in Rsol at the moment where CEE begins, seen by star2 + double omegaSpin1_pre_CE = m_Star1->Omega(); // star1 spin (before CEE) + double omegaSpin2_pre_CE = m_Star2->Omega(); // star2 spin (before CEE) bool isDonorMS = false; // check for main sequence donor if (OPTIONS->AllowMainSequenceStarToSurviveCommonEnvelope()) { // allow main sequence stars to survive CEE? @@ -1642,11 +1644,13 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { if (envelopeFlag1) { m_Star1->ResolveEnvelopeLossAndSwitch(); // resolve envelope loss for star1 and switch to new stellar type + m_Star1->SetOmega(omegaSpin1_pre_CE); // restore core spin after envelope loss m_MassTransferTrackerHistory = MT_TRACKING::CE_1_TO_2_SURV; } if (envelopeFlag2) { m_Star2->ResolveEnvelopeLossAndSwitch(); // resolve envelope loss for star1 and switch to new stellar type + m_Star2->SetOmega(omegaSpin2_pre_CE); // restore core spin after envelope loss m_MassTransferTrackerHistory = MT_TRACKING::CE_2_TO_1_SURV; } diff --git a/src/Star.cpp b/src/Star.cpp index db2cf1559..1ef0bb660 100644 --- a/src/Star.cpp +++ b/src/Star.cpp @@ -180,16 +180,11 @@ bool Star::RevertState() { * * Resolve the loss of an envelope, including switching the stellar type. * - * Keep the rotation frequency of the core equal to the pre-envelope-loss rotation frequency - * (envelope loss is rapid, so no time for angular momentum transport) - * * void ResolveEnvelopeLossAndSwitch() * */ void Star::ResolveEnvelopeLossAndSwitch() { - double omega = Omega(); (void)SwitchTo(m_Star->ResolveEnvelopeLoss(true)); - SetOmega(omega); } From 8eeefcb1d897f41a773b706aef269a29d745ace3 Mon Sep 17 00:00:00 2001 From: Veome Kapil Date: Wed, 20 Nov 2024 11:07:20 -0500 Subject: [PATCH 4/6] Updated changelog --- src/changelog.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/changelog.h b/src/changelog.h index 99b5c1ada..2aab6ca2a 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1387,7 +1387,9 @@ // - Associated code clean-up // 03.08.01 IM - Nov 19, 2024 - Defect repairs // - Multiple rotation-related fixes to 03.08.01 (units, initialisation, min->max typo, retain omega on envelope loss) +// 03.08.02 VK - Nov 20, 2024 - Defect repairs +// - Fixed behavior for core spin to be retained after envelope loss. -const std::string VERSION_STRING = "03.08.01"; +const std::string VERSION_STRING = "03.08.02"; # endif // __changelog_h__ From 49ceec769186d8c71752c7f6a126e36d04337322 Mon Sep 17 00:00:00 2001 From: Veome Kapil Date: Wed, 20 Nov 2024 11:16:03 -0500 Subject: [PATCH 5/6] Reverted old changes to match dev --- src/BaseStar.cpp | 2 -- src/BaseStar.h | 2 -- src/Star.h | 1 - 3 files changed, 5 deletions(-) diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index 01b06343b..a2d148967 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -158,7 +158,6 @@ BaseStar::BaseStar(const unsigned long int p_RandomSeed, m_StellarTypePrev = m_StellarType; m_MassPrev = m_MZAMS; m_RadiusPrev = m_RZAMS; - m_MomentOfInertiaPrevAU = CalculateMomentOfInertiaAU(); m_DtPrev = DEFAULT_INITIAL_DOUBLE_VALUE; // Lambdas @@ -4591,7 +4590,6 @@ void BaseStar::UpdateAttributesAndAgeOneTimestepPreamble(const double p_DeltaMas m_StellarTypePrev = m_StellarType; m_MassPrev = m_Mass; m_RadiusPrev = m_Radius; - m_MomentOfInertiaPrevAU = CalculateMomentOfInertiaAU(); } // the GBParams and Timescale calculations need to be done before taking the timestep - since diff --git a/src/BaseStar.h b/src/BaseStar.h index d096be158..80eeac1a3 100644 --- a/src/BaseStar.h +++ b/src/BaseStar.h @@ -144,7 +144,6 @@ class BaseStar { double Mdot() const { return m_Mdot; } double Metallicity() const { return m_Metallicity; } double MinimumCoreMass() const { return m_MinimumCoreMass; } - double MomentOfInertiaPrevAU() const { return m_MomentOfInertiaPrevAU; } double MZAMS() const { return m_MZAMS; } double Omega() const { return m_AngularMomentum / CalculateMomentOfInertiaAU(); } double OmegaCHE() const { return m_OmegaCHE; } @@ -446,7 +445,6 @@ class BaseStar { double m_DtPrev; // Previous timestep double m_MassPrev; // Previous mass (Msol) double m_RadiusPrev; // Previous radius (Rsol) - double m_MomentOfInertiaPrevAU; // Previous moment of inertia (Msol AU^2) STELLAR_TYPE m_StellarTypePrev; // Stellar type at previous timestep // Metallicity variables diff --git a/src/Star.h b/src/Star.h index 7c8579a09..65c70cd10 100755 --- a/src/Star.h +++ b/src/Star.h @@ -122,7 +122,6 @@ class Star { double Mass0() const { return m_Star->Mass0(); } double MassPrev() const { return m_Star->MassPrev(); } double Metallicity() const { return m_Star->Metallicity(); } - double MomentOfInertiaPrevAU() const { return m_Star->MomentOfInertiaPrevAU();} double MZAMS() const { return m_Star->MZAMS(); } double Omega() const { return m_Star->Omega(); } double OmegaCHE() const { return m_Star->OmegaCHE(); } From 560e5d65dcd60f085cc9808b4834d34c10b66dd1 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Mon, 25 Nov 2024 10:20:37 +1100 Subject: [PATCH 6/6] Update BaseBinaryStar.cpp Comment update --- src/BaseBinaryStar.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index ab2bd8320..7cd90c864 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1644,13 +1644,13 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { if (envelopeFlag1) { m_Star1->ResolveEnvelopeLossAndSwitch(); // resolve envelope loss for star1 and switch to new stellar type - m_Star1->SetOmega(omegaSpin1_pre_CE); // restore core spin after envelope loss + m_Star1->SetOmega(omegaSpin1_pre_CE); // keep the rotation frequency of the core equal to the pre-envelope-loss rotation frequency m_MassTransferTrackerHistory = MT_TRACKING::CE_1_TO_2_SURV; } if (envelopeFlag2) { m_Star2->ResolveEnvelopeLossAndSwitch(); // resolve envelope loss for star1 and switch to new stellar type - m_Star2->SetOmega(omegaSpin2_pre_CE); // restore core spin after envelope loss + m_Star2->SetOmega(omegaSpin2_pre_CE); // keep the rotation frequency of the core equal to the pre-envelope-loss rotation frequency m_MassTransferTrackerHistory = MT_TRACKING::CE_2_TO_1_SURV; }