From 26ffc9119400f4f91bba54935dd5c6cb35d90264 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Mon, 25 Nov 2024 20:49:24 +1100 Subject: [PATCH 1/2] Enhancements to nuclear timescale MT, several fixes, addresses issue #1285 --- src/BaseBinaryStar.cpp | 89 ++++++++++++++++++++++++------------------ src/BaseBinaryStar.h | 23 +++++++---- src/ErrorCatalog.h | 2 + src/WhiteDwarfs.h | 4 +- src/changelog.h | 6 ++- 5 files changed, 75 insertions(+), 49 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 5ed28e662..81059e914 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1979,33 +1979,46 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { jLoss = CalculateGammaAngularMomentumLoss(); // no - re-calculate angular momentum } + m_MassTransferTimescale = MASS_TRANSFER_TIMESCALE::NONE; // initial reset double betaThermal = 0.0; // fraction of mass accreted if accretion proceeds on thermal timescale - double betaNuclear = 0.0; // fraction of mass accreted if accretion proceeds on nuclear timescale + double maximumAccretionRate = 0.0; // accretion rate if accretion proceeds on thermal timescale double donorMassLossRateThermal = m_Donor->CalculateThermalMassLossRate(); - double donorMassLossRateNuclear = m_Donor->CalculateNuclearMassLossRate(); - - std::tie(std::ignore, betaThermal) = m_Accretor->CalculateMassAcceptanceRate(donorMassLossRateThermal, - m_Accretor->CalculateThermalMassAcceptanceRate(accretorRLradius), - donorIsHeRich); - std::tie(std::ignore, betaNuclear) = m_Accretor->CalculateMassAcceptanceRate(donorMassLossRateNuclear, + + std::tie(maximumAccretionRate, betaThermal) = m_Accretor->CalculateMassAcceptanceRate(donorMassLossRateThermal, m_Accretor->CalculateThermalMassAcceptanceRate(accretorRLradius), donorIsHeRich); - - m_ZetaStar = m_Donor->CalculateZetaAdiabatic(); - double zetaEquilibrium = m_Donor->CalculateZetaEquilibrium(); - - m_ZetaLobe = CalculateZetaRocheLobe(jLoss, betaNuclear); // try nuclear timescale mass transfer first - - if (m_Donor->IsOneOf(ALL_MAIN_SEQUENCE) && utils::Compare(zetaEquilibrium, m_ZetaLobe) > 0) { - m_MassLossRateInRLOF = donorMassLossRateNuclear; - m_FractionAccreted = betaNuclear; - m_MassTransferTimescale = MASS_TRANSFER_TIMESCALE::NUCLEAR; + double massDiffDonor = 0.0; + + // can the mass transfer happen on a nuclear timescale? only considering this for MS donors + if (m_Donor->IsOneOf(ALL_MAIN_SEQUENCE)) { + // technically, we do not know how much mass the accretor should gain until we do the calculation, which impacts the RL size, so we will check whether a nuclear timescale MT was feasible later + double maximumAccretedMass = maximumAccretionRate * m_Dt; + if (OPTIONS->MassTransferAccretionEfficiencyPrescription() == MT_ACCRETION_EFFICIENCY_PRESCRIPTION::THERMALLY_LIMITED) { + massDiffDonor = MassLossToFitInsideRocheLobe(this, m_Donor, m_Accretor, -1.0, maximumAccretedMass); // use root solver to determine how much mass should be lost from the donor to allow it to fit within the Roche lobe, fixed accretion amount + m_FractionAccreted = std::min(maximumAccretedMass, massDiffDonor) / massDiffDonor; + } + else { + m_FractionAccreted = maximumAccretionRate / donorMassLossRateThermal; // relevant for MT_ACCRETION_EFFICIENCY_PRESCRIPTION::FIXED_FRACTION + massDiffDonor = MassLossToFitInsideRocheLobe(this, m_Donor, m_Accretor, m_FractionAccreted, 0.0); // use root solver to determine how much mass should be lost from the donor to allow it to fit within the Roche lobe, fixed beta + } + + // check that the star really would have consistently fit into the Roche lobe + double zetaEquilibrium = m_Donor->CalculateZetaEquilibrium(); // note: value is only meaningful for MS donor + double zetaLobe = CalculateZetaRocheLobe(jLoss, m_FractionAccreted); + if (utils::Compare(zetaEquilibrium, zetaLobe) > 0 && massDiffDonor > 0.0) { // yes, it's nuclear timescale mass transfer; no need for utils::Compare here + m_MassLossRateInRLOF = massDiffDonor / m_Dt; + m_MassTransferTimescale = MASS_TRANSFER_TIMESCALE::NUCLEAR; + m_ZetaStar = zetaEquilibrium; + m_ZetaLobe = zetaLobe; + } } - else { - m_ZetaLobe = CalculateZetaRocheLobe(jLoss, betaThermal); + if (m_MassTransferTimescale != MASS_TRANSFER_TIMESCALE::NUCLEAR) { // thermal timescale mass transfer (we will check for dynamically unstable / CE mass transfer later) + m_ZetaLobe = CalculateZetaRocheLobe(jLoss, betaThermal); + m_ZetaStar = m_Donor->CalculateZetaAdiabatic(); m_MassLossRateInRLOF = donorMassLossRateThermal; m_FractionAccreted = betaThermal; m_MassTransferTimescale = MASS_TRANSFER_TIMESCALE::THERMAL; + massDiffDonor = MassLossToFitInsideRocheLobe(this, m_Donor, m_Accretor, betaThermal, 0.0); // use root solver to determine how much mass should be lost from the donor to allow it to fit within the Roche lobe } double aInitial = m_SemiMajorAxis; // semi-major axis in default units, AU, current timestep @@ -2043,9 +2056,6 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { } else { // stable MT - 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 massDiffDonor; double envMassDonor = m_Donor->Mass() - m_Donor->CoreMass(); bool isEnvelopeRemoved = false; @@ -2054,29 +2064,27 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { isEnvelopeRemoved = true; } else { // donor has no envelope - massDiffDonor = MassLossToFitInsideRocheLobe(this, m_Donor, m_Accretor, m_FractionAccreted); // use root solver to determine how much mass should be lost from the donor to allow it to fit within the Roche lobe - if (massDiffDonor <= 0.0) { // no root found - // if donor cannot lose mass to fit inside Roche lobe, the only viable action is to enter CE phase + // if donor cannot lose mass to fit inside Roche lobe, the only viable action is to enter CE phase -- but this probably should not happen... + SHOW_WARN(ERROR::UNEXPECTED_ROOT_FINDER_FAILURE); m_CEDetails.CEEnow = true; // flag CE } else { // have required mass loss - if (utils::Compare(m_MassLossRateInRLOF,donorMassLossRateNuclear) == 0) // if transferring mass on nuclear timescale, limit mass loss amount to rate * timestep (thermal timescale MT always happens in one timestep) - massDiffDonor = std::min(massDiffDonor, m_MassLossRateInRLOF * m_Dt); massDiffDonor = -massDiffDonor; // set mass difference m_Donor->UpdateMinimumCoreMass(); // reset the minimum core mass following case A } - } if (!m_CEDetails.CEEnow) { // CE flagged? // 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 m_Donor->SetMassTransferDiffAndResolveWDShellChange(massDiffDonor); // set new mass of donor m_Accretor->SetMassTransferDiffAndResolveWDShellChange(massGainAccretor); // set new mass of accretor aFinal = CalculateMassTransferOrbit(m_Donor->Mass(), massDiffDonor, *m_Accretor, m_FractionAccreted); // calculate new orbit + m_aMassTransferDiff = aFinal - aInitial; // set change in orbit (semi-major axis) STELLAR_TYPE stellarTypeDonor = m_Donor->StellarType(); // donor stellar type before resolving envelope loss @@ -2109,15 +2117,16 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { * Uses boost::math::tools::bracket_and_solve_root() * * - * double MassLossToFitInsideRocheLobe(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted) + * double MassLossToFitInsideRocheLobe(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted, double p_MaximumAccretedMass) * * @param [IN] p_Binary (Pointer to) The binary star under examination * @param [IN] p_Donor (Pointer to) The star donating mass * @param [IN] p_Accretor (Pointer to) The star accreting mass - * @param [IN] p_FractionAccreted The fraction of the donated mass accreted by the accretor + * @param [IN] p_FractionAccreted The fraction of the donated mass accreted by the accretor (for thermal timescale accretion) + * @param [IN] p_MaximumAccretedMass The total amount of mass that can be accreted (for nuclear timescale accretion, p_FractionAccreted should be negative for this to be used) * @return Root found: will be -1.0 if no acceptable real root found */ -double BaseBinaryStar::MassLossToFitInsideRocheLobe(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted) { +double BaseBinaryStar::MassLossToFitInsideRocheLobe(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted, double p_MaximumAccretedMass) { const boost::uintmax_t maxit = ADAPTIVE_RLOF_MAX_ITERATIONS; // Limit to maximum iterations. boost::uintmax_t it = maxit; // Initially our chosen max iterations, but updated with actual. @@ -2136,19 +2145,21 @@ double BaseBinaryStar::MassLossToFitInsideRocheLobe(BaseBinaryStar *p_Binary, Bi double factorFrac = ADAPTIVE_RLOF_SEARCH_FACTOR_FRAC; // search step size factor fractional part double factor = 1.0 + factorFrac; // factor to determine search step size (size = guess * factor) - + std::pair root(-1.0, -1.0); // initialise root - default return std::size_t tries = 0; // number of tries bool done = false; // finished (found root or exceed maximum tries)? ERROR error = ERROR::NONE; - RadiusEqualsRocheLobeFunctor func = RadiusEqualsRocheLobeFunctor(p_Binary, p_Donor, p_Accretor, p_FractionAccreted, &error); // no need to check error here + RadiusEqualsRocheLobeFunctor func = RadiusEqualsRocheLobeFunctor(p_Binary, p_Donor, p_Accretor, p_FractionAccreted, p_MaximumAccretedMass, &error); while (!done) { // while no error and acceptable root found - double semiMajorAxis = p_Binary->CalculateMassTransferOrbit(p_Donor->Mass(), -guess , *p_Accretor, p_FractionAccreted); - double RLRadius = semiMajorAxis * (1.0 - p_Binary->Eccentricity()) * CalculateRocheLobeRadius_Static(p_Donor->Mass() - guess, p_Accretor->Mass() + (p_Binary->FractionAccreted() * guess)) * AU_TO_RSOL; - double radiusAfterMassLoss = p_Donor->CalculateRadiusOnPhaseTau(p_Donor->Mass()-guess, p_Donor->Tau()); - bool isRising = radiusAfterMassLoss > RLRadius ? true : false; // guess for direction of search - + bool isRising = true; //guess for direction of search + // while the change in the functor at guess may be more appropriate -- something like + // isRising = (RLRadiusGuess-radiusAfterMassLoss) > (RLRadius - radius)? true : false; + // or isRising = func((const double)guess) >= func((const double)guess * factor) ? false : true; + // -- this choice is more robust given that we will be taking smaller steps anyway (following factor reduction) + // if a bigger step does not find a solution + // run the root finder // regardless of any exceptions or errors, display any problems as a warning, then @@ -2439,7 +2450,7 @@ void BaseBinaryStar::ResolveMassChanges() { // (a sign that ResolveEnvelopeLossAndSwitch() has been called after the full envelope was stripped) // no need to resolve mass changes if the mass has already been updated // calculate mass change due to winds and mass transfer - if (utils::Compare(m_Star1->MassPrev(), m_Star1->Mass()) == 0) { // mass already updated? + if (utils::Compare(m_Star2->MassPrev(), m_Star2->Mass()) == 0) { // mass already updated? // no - resolve mass changes double massChange = m_Star2->MassLossDiff() + m_Star2->MassTransferDiff(); // mass change due to winds and mass transfer if (utils::Compare(massChange, 0.0) != 0) { // winds/mass transfer changes mass? diff --git a/src/BaseBinaryStar.h b/src/BaseBinaryStar.h index 9ff04d830..194dd5ed0 100644 --- a/src/BaseBinaryStar.h +++ b/src/BaseBinaryStar.h @@ -539,12 +539,13 @@ class BaseBinaryStar { * * * Constructor: initialise the class - * template RadiusEqualsRocheLobeFunctor(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted) + * template RadiusEqualsRocheLobeFunctor(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted, double p_MaximumAccretedMass, ERROR *p_Error) * * @param [IN] p_Binary (Pointer to) The binary star under examination * @param [IN] p_Donor (Pointer to) The star donating mass * @param [IN] p_Accretor (Pointer to) The star accreting mass - * @param [IN] p_FractionAccreted The fraction of the donated mass accreted by the accretor + * @param [IN] p_FractionAccreted The fraction of the donated mass accreted by the accretor (for thermal timescale accretion) + * @param [IN] p_MaximumAccretedMass The total amount of mass that can be accreted (for nuclear timescale accretion, p_FractionAccreted should be negative for this to be used) * @param [IN] p_Error (Address of variable to record) Error encountered in functor * * Function: calculate radius difference after mass loss @@ -555,12 +556,13 @@ class BaseBinaryStar { */ template struct RadiusEqualsRocheLobeFunctor { - RadiusEqualsRocheLobeFunctor(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted, ERROR *p_Error) { + RadiusEqualsRocheLobeFunctor(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted, double p_MaximumAccretedMass, ERROR *p_Error) { m_Binary = p_Binary; m_Donor = p_Donor; m_Accretor = p_Accretor; m_Error = p_Error; m_FractionAccreted = p_FractionAccreted; + m_MaximumAccretedMass = p_MaximumAccretedMass; } T operator()(double const& p_dM) { @@ -572,11 +574,17 @@ class BaseBinaryStar { double donorMass = m_Donor->Mass(); double accretorMass = m_Accretor->Mass(); - double semiMajorAxis = m_Binary->CalculateMassTransferOrbit(m_Donor->Mass(), -p_dM , *m_Accretor, m_FractionAccreted); - double RLRadius = semiMajorAxis * (1.0 - m_Binary->Eccentricity()) * CalculateRocheLobeRadius_Static(donorMass - p_dM, accretorMass + (m_Binary->FractionAccreted() * p_dM)) * AU_TO_RSOL; + // beta is the actual accretion efficiency; if p_FractionAccreted is negative (placeholder + // for nuclear timescale accretion efficiency, for which the total accretion mass over the + // duration of the timestep is known), then the ratio of the maximum allowed accreted + // mass / donated mass is used + double beta = (utils::Compare(m_FractionAccreted, 0.0) >=0 ) ? m_FractionAccreted : std::min(m_MaximumAccretedMass/p_dM, 1.0); + + double semiMajorAxis = m_Binary->CalculateMassTransferOrbit(donorMass, -p_dM , *m_Accretor, beta); + double RLRadius = semiMajorAxis * (1.0 - m_Binary->Eccentricity()) * CalculateRocheLobeRadius_Static(donorMass - p_dM, accretorMass + (beta * p_dM)) * AU_TO_RSOL; double radiusAfterMassLoss = m_Donor->CalculateRadiusOnPhaseTau(donorMass-p_dM, m_Donor->Tau()); - + return (RLRadius - radiusAfterMassLoss); } private: @@ -585,10 +593,11 @@ class BaseBinaryStar { BinaryConstituentStar *m_Accretor; ERROR *m_Error; double m_FractionAccreted; + double m_MaximumAccretedMass; }; - double MassLossToFitInsideRocheLobe(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted); + double MassLossToFitInsideRocheLobe(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted, double p_MaximumAccretedMass); double OmegaAfterSynchronisation(const double p_M1, const double p_M2, const double p_I1, const double p_I2, const double p_Ltot, const double p_Guess); diff --git a/src/ErrorCatalog.h b/src/ErrorCatalog.h index 87f47e905..fcf66e6e4 100644 --- a/src/ErrorCatalog.h +++ b/src/ErrorCatalog.h @@ -142,6 +142,7 @@ enum class ERROR: int { UNEXPECTED_PROGRAM_OPTION, // unexpected program option UNEXPECTED_PROPERTY, // unexpected property UNEXPECTED_PROPERTY_TYPE, // unexpected property type + UNEXPECTED_ROOT_FINDER_FAILURE, // unexpected root finder failure UNEXPECTED_SN_EVENT, // unexpected supernova event in this context UNEXPECTED_STELLAR_PROPERTY, // unexpected stellar property UNEXPECTED_STELLAR_PROPERTY_TYPE, // unexpected stellar property type @@ -314,6 +315,7 @@ const COMPASUnorderedMap> ERROR_CATA { ERROR::UNEXPECTED_PROGRAM_OPTION, { ERROR_SCOPE::ALWAYS, "Unexpected program option" }}, { ERROR::UNEXPECTED_PROPERTY, { ERROR_SCOPE::ALWAYS, "Unexpected property" }}, { ERROR::UNEXPECTED_PROPERTY_TYPE, { ERROR_SCOPE::ALWAYS, "Unexpected property type" }}, + { ERROR::UNEXPECTED_ROOT_FINDER_FAILURE, { ERROR_SCOPE::ALWAYS, "Unexpected root finder failure" }}, { ERROR::UNEXPECTED_SN_EVENT, { ERROR_SCOPE::ALWAYS, "Unexpected supernova event in this context" }}, { ERROR::UNEXPECTED_STELLAR_PROPERTY, { ERROR_SCOPE::ALWAYS, "Unexpected stellar property" }}, { ERROR::UNEXPECTED_STELLAR_PROPERTY_TYPE, { ERROR_SCOPE::ALWAYS, "Unexpected stellar property type" }}, diff --git a/src/WhiteDwarfs.h b/src/WhiteDwarfs.h index f4eecaf41..ad4331986 100644 --- a/src/WhiteDwarfs.h +++ b/src/WhiteDwarfs.h @@ -74,10 +74,10 @@ class WhiteDwarfs: virtual public BaseStar, public Remnants { double Calculatel0Ritter() const { return (m_Metallicity > 0.01) ? 1995262.3 : 31622.8; } // Luminosity constant which depends on metallicity in Ritter 1999, eq 10 virtual DBL_DBL CalculateMassAcceptanceRate(const double p_DonorMassRate, - const bool p_IsHeRich) { return std::make_tuple(0.0, 0.0); } // Should never be called JR: is this true? Not implemented in ONeWD clas? + const bool p_IsHeRich) { return std::make_tuple(0.0, 0.0); } DBL_DBL CalculateMassAcceptanceRate(const double p_DonorMassRate, const double p_AccretorMassRate, - const bool p_IsHeRich) { return CalculateMassAcceptanceRate(p_DonorMassRate, p_IsHeRich); } // Ignore the input accretion rate for WDs + const bool p_IsHeRich) { return CalculateMassAcceptanceRate(p_DonorMassRate, p_IsHeRich); } double CalculateXRitter() const { return (m_Metallicity > 0.01) ? 0.7 : 0.8 ; } // Assumed Hydrogen-mass fraction diff --git a/src/changelog.h b/src/changelog.h index 5e2d4490f..40cdd9cdc 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1395,6 +1395,10 @@ // - Fixed behavior for core spin to be retained after envelope loss // 03.08.04 IM - Nov 25, 2024 - Defect repair: // - Recalculate timescales when updating stellar age after mass loss (addresses issue #1231) -const std::string VERSION_STRING = "03.08.04"; +// 03.09.00 IM - Nov 25, 2024 - Defect repair, enhancement +// - The nuclear timescale mass transfer rate is now set by the requirement that the star ends the time step just filling its Roche lobe (addresses issue #1285) +// - Fix an issue with the root finder for fitting into the RL that led to artificial failures to find a root +// - Fix issue (likely introduced in 03.08.00) with the accretor not gaining mass appropriately +const std::string VERSION_STRING = "03.09.00"; # endif // __changelog_h__ From 506b6c72be560fd7bd12316ea8085063827b5a11 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Mon, 25 Nov 2024 20:59:16 +1100 Subject: [PATCH 2/2] Cleaning extra space --- src/BaseBinaryStar.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 81059e914..7c5135e8c 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -2451,7 +2451,7 @@ void BaseBinaryStar::ResolveMassChanges() { // no need to resolve mass changes if the mass has already been updated // calculate mass change due to winds and mass transfer if (utils::Compare(m_Star2->MassPrev(), m_Star2->Mass()) == 0) { // mass already updated? - // no - resolve mass changes + // no - resolve mass changes double massChange = m_Star2->MassLossDiff() + m_Star2->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