From 6557a809753d9ce049fc3c943c70a8427ebe9e59 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Tue, 2 Sep 2025 22:43:37 +1000 Subject: [PATCH 1/3] First implementation of HAMSTARS mass transfer efficiency --- .../preprocessing/compasConfigDefault.yaml | 4 ++-- .../Program options/program-options-list-defaults.rst | 2 +- online-docs/pages/whats-new.rst | 4 ++++ src/BaseBinaryStar.cpp | 2 +- src/BaseStar.cpp | 9 +++++++-- src/Options.cpp | 2 +- src/changelog.h | 4 +++- src/typedefs.h | 5 +++-- 8 files changed, 22 insertions(+), 10 deletions(-) diff --git a/compas_python_utils/preprocessing/compasConfigDefault.yaml b/compas_python_utils/preprocessing/compasConfigDefault.yaml index fdc1e7c68..a603d6f13 100644 --- a/compas_python_utils/preprocessing/compasConfigDefault.yaml +++ b/compas_python_utils/preprocessing/compasConfigDefault.yaml @@ -1,5 +1,5 @@ ##~!!~## COMPAS option values -##~!!~## File Created Tue Aug 19 15:39:23 2025 by COMPAS v03.25.00 +##~!!~## File Created Tue Sep 2 22:42:52 2025 by COMPAS v03.26.00 ##~!!~## ##~!!~## The default COMPAS YAML file (``compasConfigDefault.yaml``), as distributed, has ##~!!~## all COMPAS option entries commented so that the COMPAS default value for the @@ -264,7 +264,7 @@ stringChoices: # --critical-mass-ratio-prescription: 'NONE' # Default: 'NONE' # Options: ['HURLEY_HJELLMING_WEBBINK','GE','GE_IC','CLAEYS','NONE'] # --stellar-zeta-prescription: 'SOBERMAN' # Default: 'SOBERMAN' # Options: ['ARBITRARY','HURLEY','SOBERMAN'] # --mass-transfer-angular-momentum-loss-prescription: 'ISOTROPIC' # Default: 'ISOTROPIC' # Options: ['ARBITRARY','KLENCKI_LINEAR','MACLEOD_LINEAR','CIRCUMBINARY','ISOTROPIC','JEANS'] -# --mass-transfer-accretion-efficiency-prescription: 'THERMAL' # Default: 'THERMAL' # Options: ['FIXED','THERMAL'] +# --mass-transfer-accretion-efficiency-prescription: 'THERMAL' # Default: 'THERMAL' # Options: ['HAMSTARS','FIXED','THERMAL'] # --mass-transfer-rejuvenation-prescription: 'STARTRACK' # Default: 'STARTRACK' # Options: ['STARTRACK','HURLEY'] # --mass-transfer-thermal-limit-accretor-multiplier: 'CFACTOR' # Default: 'CFACTOR' # Options: ['ROCHELOBE','CFACTOR'] # --response-to-spin-up: 'TRANSFER_TO_ORBIT' # Default: 'TRANSFER_TO_ORBIT' # Options: ['NO_LIMIT','KEPLERIAN_LIMIT','TRANSFER_TO_ORBIT'] diff --git a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst index 6c5df47a5..6c1d05b00 100644 --- a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst +++ b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst @@ -887,7 +887,7 @@ DEPRECATION NOTICE: this option has been deprecated and will soon be removed. Pl **--mass-transfer-accretion-efficiency-prescription** |br| Mass transfer accretion efficiency prescription. |br| -Options: { THERMAL, FIXED } |br| +Options: { THERMAL, FIXED, HAMSTARS } |br| Default = THERMAL **--mass-transfer-angular-momentum-loss-prescription** |br| diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index e20ac701e..221a28393 100644 --- a/online-docs/pages/whats-new.rst +++ b/online-docs/pages/whats-new.rst @@ -3,6 +3,10 @@ What's new Following is a brief list of important updates to the COMPAS code. A complete record of changes can be found in the file ``changelog.h``. +**03.26.00 September 2, 2025** + +* Added HAMSTARS mass transfer efficiency prescription + **03.25.00 August 19, 2025** * Added KLENCKI_LINEAR AM loss, which is linear in the specific AM gamma instead of the orbital separation (as in MACLEOD_LINEAR). diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index c6046f4a9..51c13d18d 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -2092,7 +2092,7 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { // can the mass transfer happen on a nuclear timescale? if (m_Donor->IsOneOf(NON_COMPACT_OBJECTS)) { // if MT_ACCRETION_EFFICIENCY_PRESCRIPTION::FIXED_FRACTION, then CalculateMassAcceptanceRate() already computed the correct m_FractionAccreted and massDiffDonor (no difference between nuclear and thermal timescale MT) - if (OPTIONS->MassTransferAccretionEfficiencyPrescription() == MT_ACCRETION_EFFICIENCY_PRESCRIPTION::THERMALLY_LIMITED) { + if (OPTIONS->MassTransferAccretionEfficiencyPrescription() == MT_ACCRETION_EFFICIENCY_PRESCRIPTION::THERMALLY_LIMITED || OPTIONS->MassTransferAccretionEfficiencyPrescription() == MT_ACCRETION_EFFICIENCY_PRESCRIPTION::HAMSTARS) { // 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 massDiffDonor = MassLossToFitInsideRocheLobe(this, m_Donor, m_Accretor, -1.0, m_Dt); // use root solver to determine how much mass should be lost from the donor to allow it to fit within the Roche lobe, estimating accretion efficiency based on a mass donation rate of massDiffDonor/m_Dt for self-consistency diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index bda5dec0b..4fc20e175 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -2713,12 +2713,17 @@ DBL_DBL BaseStar::CalculateMassAcceptanceRate(const double p_DonorMassRate, cons switch (OPTIONS->MassTransferAccretionEfficiencyPrescription()) { - case MT_ACCRETION_EFFICIENCY_PRESCRIPTION::THERMALLY_LIMITED: // thermally limited mass transfer: - + case MT_ACCRETION_EFFICIENCY_PRESCRIPTION::THERMALLY_LIMITED: // thermally limited mass transfer acceptanceRate = min(OPTIONS->MassTransferCParameter() * p_AccretorMassRate, p_DonorMassRate); fractionAccreted = acceptanceRate / p_DonorMassRate; break; + case MT_ACCRETION_EFFICIENCY_PRESCRIPTION::HAMSTARS: // thermally limited mass transfer, following Lau+, 2024 + // the mass transfer C parameter is fit using the data from Figure (1) of Lau+, 2024 + acceptanceRate = min(PPOW(10.0, (4.0 / PPOW((Mass() + 0.2), 0.3) - 0.6)) * p_AccretorMassRate, p_DonorMassRate); + fractionAccreted = acceptanceRate / p_DonorMassRate; + break; + case MT_ACCRETION_EFFICIENCY_PRESCRIPTION::FIXED_FRACTION: // fixed fraction of mass accreted, as in StarTrack fractionAccreted = OPTIONS->MassTransferFractionAccreted(); acceptanceRate = min(p_DonorMassRate, fractionAccreted * p_DonorMassRate); diff --git a/src/Options.cpp b/src/Options.cpp index e00917889..bc2d3c2a1 100644 --- a/src/Options.cpp +++ b/src/Options.cpp @@ -2442,7 +2442,7 @@ std::string Options::OptionValues::CheckAndSetOptions() { if (m_UseMassTransfer && !DEFAULTED("mass-transfer-accretion-efficiency-prescription")) { // mass transfer accretion efficiency prescription std::tie(found, m_MassTransferAccretionEfficiencyPrescription.type) = utils::GetMapKey(m_MassTransferAccretionEfficiencyPrescription.typeString, MT_ACCRETION_EFFICIENCY_PRESCRIPTION_LABEL, m_MassTransferAccretionEfficiencyPrescription.type); - COMPLAIN_IF(!found, "Unknown Mass Transfer Angular Momentum Loss Prescription"); + COMPLAIN_IF(!found, "Unknown Mass Transfer Efficiency Prescription"); } if (m_UseMassTransfer && !DEFAULTED("mass-transfer-angular-momentum-loss-prescription")) { // mass transfer angular momentum loss prescription diff --git a/src/changelog.h b/src/changelog.h index 026aa60d2..b75a33529 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1664,6 +1664,8 @@ // - Corrected calculation for Hurley Gamma constant C (C_GAMMA - see Hurley et al. 2000, just after eq 23, should use a(75) <= 1.0, not a(75) == 1.0 - confirmed in BSE Fortran source) // - Added abs() to gamma calculation in Mainsequence.cpp::CalculateGamma() (per BSE Fortran source) // - Clamped gamma to [0.0, gamma] in Mainsequence.cpp::CalculateGamma() (per discussion just after eq 23 - confirmed in BSE Fortran source) +// 03.26.00 IM - September 2, 2025 - Enhancement: +// - First (simplified) implementation of the Lau+ (2024) Hamstars thermally limited accretion prescription // // Version string format is MM.mm.rr, where // @@ -1674,7 +1676,7 @@ // if MM is incremented, set mm and rr to 00, even if defect repairs and minor enhancements were also made // if mm is incremented, set rr to 00, even if defect repairs were also made -const std::string VERSION_STRING = "03.25.02"; +const std::string VERSION_STRING = "03.26.00"; # endif // __changelog_h__ diff --git a/src/typedefs.h b/src/typedefs.h index cb2ae6db9..eb3e82907 100755 --- a/src/typedefs.h +++ b/src/typedefs.h @@ -643,10 +643,11 @@ const COMPASUnorderedMap METALLICITY_DIST }; // mass transfer accretion efficiency prescriptions -enum class MT_ACCRETION_EFFICIENCY_PRESCRIPTION: int { THERMALLY_LIMITED, FIXED_FRACTION }; +enum class MT_ACCRETION_EFFICIENCY_PRESCRIPTION: int { THERMALLY_LIMITED, FIXED_FRACTION, HAMSTARS }; const COMPASUnorderedMap MT_ACCRETION_EFFICIENCY_PRESCRIPTION_LABEL = { { MT_ACCRETION_EFFICIENCY_PRESCRIPTION::THERMALLY_LIMITED, "THERMAL" }, - { MT_ACCRETION_EFFICIENCY_PRESCRIPTION::FIXED_FRACTION, "FIXED" } + { MT_ACCRETION_EFFICIENCY_PRESCRIPTION::FIXED_FRACTION, "FIXED" }, + { MT_ACCRETION_EFFICIENCY_PRESCRIPTION::HAMSTARS, "HAMSTARS"} }; // mass transfer angular momentum loss prescriptions From d3790b5244d8ee70113b2510ad98b740c817b0e2 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Tue, 9 Sep 2025 21:22:37 +1000 Subject: [PATCH 2/3] Corrected errors in combining OB and WR winds in CH::CalculateMassLossRateBelczynski2010() and CH::CalculateMassLossRateFractionOB() [previously CalculateMassLossRateWeightOB()] --- src/CH.cpp | 57 +++++++++++++++++++++++-------------------------- src/CH.h | 2 +- src/changelog.h | 3 ++- 3 files changed, 30 insertions(+), 32 deletions(-) diff --git a/src/CH.cpp b/src/CH.cpp index fa01bf707..f75163634 100755 --- a/src/CH.cpp +++ b/src/CH.cpp @@ -286,36 +286,33 @@ void CH::UpdateAgeAfterMassLoss() { /* - * Calculate the weight for the OB and WR mass loss prescriptions + * CalculateMassLossFractionOB * - * According to the prescription described in Yoon et al. 2006 (and Szecsi et al. 2015) - * Use OB mass loss rates until Y1 = 0.55, WR mass loss rates above Y2 = 0.7, and - * linearly interpolate for Y1 < Ys < Y2 where Ys is the surface helium fraction. - * - * Since we don't have Ys by default in our models, and detailed models show Ys - * rises ~linearly from 0.2 to 1.0 across the main-sequence, here we simply - * use tau = t / tMS as a proxy for Ys. + * @brief + * Calculate the fraction of mass loss attributable to OB mass loss, per Yoon et al. 2006 + * + * The model described in Yoon et al. 2006 (also Szecsi et al. 2015) uses OB mass loss while the + * He surface abundance is below 0.55, WR mass loss when the surface He abundance is above 0.7, + * and linearly interpolate when the He surface abundance is between those limits. + * + * This function calculates the fraction of mass loss attributable to OB mass loss, based on + * the He surface abundance and the abundance limits described in Yoon et al. 2006. The value + * returned will be 1.0 if 100% of the mass loss is attributable to OB mass lass, 0.0 if 100% of + * the mass loss is attributable to WR mass loss, and in the range (0.0, 1.0) if the mass loss is + * a mix of OB and WR. * - * CalculateMassLossRateWeightOB(const double p_HeliumAbundanceSurface) * - * @param [IN] p_HeliumAbundanceSurface Surface helium abundance Ys + * double CalculateMassLossFractionOB(const double p_HeAbundanceSurface) const * - * @return weight for OB mass loss rate + * @param p_HeAbundanceSurface Helium abundance at the surface of the star + * @return Fraction of mass loss attributable to OB mass loss */ -double CH::CalculateMassLossRateWeightOB(const double p_HeliumAbundanceSurface) { +double CH::CalculateMassLossFractionOB(const double p_HeAbundanceSurface) const { - const double y1 = 0.55; - const double y2 = 0.70; - double Ys = p_HeliumAbundanceSurface; // Get surface helium abundance - double weight = 0.0; + constexpr double limOB = 0.55; // per Yoon et al. 2006 + constexpr double limWR = 0.70; // per Yoon et al. 2006 - - // Calculate the weight as a function of Ys according to the prescription from Yoon et al. 2006 - if (Ys < y1) weight = 1.0; - else if (Ys > y2) weight = 0.0; - else weight = (y2 - Ys) / y2 - y1; - - return weight; + return std::min(1.0, std::max (0.0, (limWR - p_HeAbundanceSurface) / (limWR - limOB))); } @@ -337,7 +334,7 @@ double CH::CalculateMassLossRateBelczynski2010() { double Mdot = 0.0; double MdotOB = 0.0; double MdotWR = 0.0; - double weight = 1.0; // Initialised to 1.0 to allow us to use the OB mass loss rate by default + double fractionOB = 1.0; // Initialised to 1.0 to allow us to use the OB mass loss rate by default // Calculate OB mass loss rate according to the Vink et al. formalism MdotOB = BaseStar::CalculateMassLossRateBelczynski2010(); @@ -348,16 +345,16 @@ double CH::CalculateMassLossRateBelczynski2010() { // Calculate WR mass loss rate MdotWR = BaseStar::CalculateMassLossRateWolfRayetZDependent(0.0); - // Calculate weight for combining these into total mass-loss rate - weight = CalculateMassLossRateWeightOB(m_HeliumAbundanceSurface); + // Calculate fraction for combining these into total mass-loss rate + fractionOB = CalculateMassLossRateFractionOB(m_HeliumAbundanceSurface); } + // Finally, combine each of these prescriptions according to the OB wind fraction + Mdot = (fractionOB * MdotOB) + ((1.0 - fractionOB) * MdotWR); + // Set dominant mass loss rate - m_DominantMassLossRate = weight == 0 ? MASS_LOSS_TYPE::WR : MASS_LOSS_TYPE::OB; - - // Finally, combine each of these prescriptions according to the weight - Mdot = (weight * MdotOB) + ((1.0 - weight) * MdotWR); + m_DominantMassLossRate = (fractionOB * MdotOB) > ((1.0 - fractionOB) * MdotWR) ? MASS_LOSS_TYPE::OB : MASS_LOSS_TYPE::WR; // Enhance mass loss rate due to rotation Mdot *= CalculateMassLossRateEnhancementRotation(); diff --git a/src/CH.h b/src/CH.h index dd6ad1cb5..26cc4d3e0 100755 --- a/src/CH.h +++ b/src/CH.h @@ -74,7 +74,7 @@ class CH: virtual public BaseStar, public MS_gt_07 { // Mass loss rate double CalculateMassLossRateBelczynski2010(); double CalculateMassLossRateMerritt2025(); - double CalculateMassLossRateWeightOB(const double p_HeliumAbundanceSurface); + double CalculateMassLossRateFractionOB(const double p_HeliumAbundanceSurface); // Radius double CalculateRadiusOnPhase() const { return m_RZAMS; } // Constant from birth diff --git a/src/changelog.h b/src/changelog.h index b75a33529..ecc053979 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1664,8 +1664,9 @@ // - Corrected calculation for Hurley Gamma constant C (C_GAMMA - see Hurley et al. 2000, just after eq 23, should use a(75) <= 1.0, not a(75) == 1.0 - confirmed in BSE Fortran source) // - Added abs() to gamma calculation in Mainsequence.cpp::CalculateGamma() (per BSE Fortran source) // - Clamped gamma to [0.0, gamma] in Mainsequence.cpp::CalculateGamma() (per discussion just after eq 23 - confirmed in BSE Fortran source) -// 03.26.00 IM - September 2, 2025 - Enhancement: +// 03.26.00 IM - September 2, 2025 - Enhancement, defect repairs: // - First (simplified) implementation of the Lau+ (2024) Hamstars thermally limited accretion prescription +// - Corrected errors in combining OB and WR winds in CH::CalculateMassLossRateBelczynski2010() and CH::CalculateMassLossRateFractionOB() [previously CalculateMassLossRateWeightOB()] // // Version string format is MM.mm.rr, where // From 326930cebe6f7f16167bb0f81ec01cb5a22fca47 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Tue, 9 Sep 2025 21:28:24 +1000 Subject: [PATCH 3/3] Typo fixes --- src/CH.cpp | 14 +++++++------- src/CH.h | 2 +- src/changelog.h | 2 +- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/CH.cpp b/src/CH.cpp index f75163634..ab2b59b6c 100755 --- a/src/CH.cpp +++ b/src/CH.cpp @@ -346,7 +346,7 @@ double CH::CalculateMassLossRateBelczynski2010() { MdotWR = BaseStar::CalculateMassLossRateWolfRayetZDependent(0.0); // Calculate fraction for combining these into total mass-loss rate - fractionOB = CalculateMassLossRateFractionOB(m_HeliumAbundanceSurface); + fractionOB = CalculateMassLossFractionOB(m_HeliumAbundanceSurface); } @@ -381,7 +381,7 @@ double CH::CalculateMassLossRateMerritt2025() { double Mdot = 0.0; double MdotOB = 0.0; double MdotWR = 0.0; - double weight = 1.0; // Initialised to 1.0 to allow us to use the OB mass loss rate by default + double fractionOB = 1.0; // Initialised to 1.0 to allow us to use the OB mass loss rate by default // Calculate OB mass loss rate according to the chosen prescription MdotOB = BaseStar::CalculateMassLossRateOB(OPTIONS->OBMassLossPrescription()); @@ -397,14 +397,14 @@ double CH::CalculateMassLossRateMerritt2025() { delete clone; clone = nullptr; // return the memory allocated for the clone // Calculate weight for combining these into total mass-loss rate - weight = CalculateMassLossRateWeightOB(m_HeliumAbundanceSurface); + fractionOB = CalculateMassLossFractionOB(m_HeliumAbundanceSurface); } + // Finally, combine each of these prescriptions according to the OB wind fraction + Mdot = (fractionOB * MdotOB) + ((1.0 - fractionOB) * MdotWR); + // Set dominant mass loss rate - m_DominantMassLossRate = weight == 0 ? MASS_LOSS_TYPE::WR : MASS_LOSS_TYPE::OB; - - // Finally, combine each of these prescriptions according to the weight - Mdot = (weight * MdotOB) + ((1.0 - weight) * MdotWR); + m_DominantMassLossRate = (fractionOB * MdotOB) > ((1.0 - fractionOB) * MdotWR) ? MASS_LOSS_TYPE::OB : MASS_LOSS_TYPE::WR; // Enhance mass loss rate due to rotation Mdot *= CalculateMassLossRateEnhancementRotation(); diff --git a/src/CH.h b/src/CH.h index 26cc4d3e0..9a8ff7300 100755 --- a/src/CH.h +++ b/src/CH.h @@ -74,7 +74,7 @@ class CH: virtual public BaseStar, public MS_gt_07 { // Mass loss rate double CalculateMassLossRateBelczynski2010(); double CalculateMassLossRateMerritt2025(); - double CalculateMassLossRateFractionOB(const double p_HeliumAbundanceSurface); + double CalculateMassLossFractionOB(const double p_HeAbundanceSurface) const; // Radius double CalculateRadiusOnPhase() const { return m_RZAMS; } // Constant from birth diff --git a/src/changelog.h b/src/changelog.h index ecc053979..787d40af5 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1666,7 +1666,7 @@ // - Clamped gamma to [0.0, gamma] in Mainsequence.cpp::CalculateGamma() (per discussion just after eq 23 - confirmed in BSE Fortran source) // 03.26.00 IM - September 2, 2025 - Enhancement, defect repairs: // - First (simplified) implementation of the Lau+ (2024) Hamstars thermally limited accretion prescription -// - Corrected errors in combining OB and WR winds in CH::CalculateMassLossRateBelczynski2010() and CH::CalculateMassLossRateFractionOB() [previously CalculateMassLossRateWeightOB()] +// - Corrected errors in combining OB and WR winds in CH::CalculateMassLossRateBelczynski2010(), CalculateMassLossRateMerritt2025() and CH::CalculateMassLossFractionOB() [previously CalculateMassLossRateWeightOB()] // // Version string format is MM.mm.rr, where //