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 d362c6d05..5a380bd6d 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 @@ -847,7 +847,7 @@ Options: { JEANS, ISOTROPIC, CIRCUMBINARY, MACLEOD_LINEAR, ARBITRARY } |br| Default = ISOTROPIC **--mass-transfer-fa** |br| -Mass Transfer fraction accreted. |br| +Mass Transfer fraction accreted (beta). |br| Used when ``--mass-transfer-accretion-efficiency-prescription = FIXED_FRACTION``. |br| Default = 0.5 @@ -1122,10 +1122,11 @@ Enable mass loss due to pulsational-pair-instability (PPI). |br| Default = TRUE **--pulsational-pair-instability-prescription** |br| -Pulsational pair instability prescription. |br| +Pulsational pair instability prescription (only relevant when using ``--pulsational-pair-instability``). |br| Options: { HENDRIKS, COMPAS, STARTRACK, MARCHANT, FARMER } |br| ``HENDRIKS`` implements the prescription from Hendriks et al. 2023 |br| -``COMPAS``, ``STARTRACK`` and ``MARCHANT`` follow Woosley 2017, Belczynski et al. 2016, and Marchant et al. 2018, all as implemented in Stevenson et al. 2019. |br| +``COMPAS``, ``STARTRACK`` and ``MARCHANT`` follow Woosley 2017, Belczynski et al. 2016, and Marchant et al. 2018, +all as implemented in Stevenson et al. 2019. |br| ``FARMER`` follows Farmer et al. 2019 |br| Default = MARCHANT diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index dfaef4514..a5fd6420a 100644 --- a/online-docs/pages/whats-new.rst +++ b/online-docs/pages/whats-new.rst @@ -3,6 +3,13 @@ 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.09.00 Nov 25, 2024** + +Improved nuclear timescale mass transfer: the nuclear timescale mass transfer rate is now set by the requirement that the star +ends the time step just filling its Roche lobe. +Fixed several significant mass-transfer issues, such as accretors not gaining mass appropriately and failures +in the root solver for fitting the star into the Roche lobe that were leading to artificial common envelopes and mergers. + **03.08.02 Nov 18, 2024** Updated implementation of the mass transfer stability critical mass ratio tables from the team of Hongwei Ge. diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 0025303b4..951daddd0 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1644,8 +1644,13 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { } } - if (!m_Flags.stellarMerger) { - + if (utils::Compare(m_SemiMajorAxis, 0.0) <= 0 || utils::Compare(m_Star1->CalculateRemnantRadius() + m_Star2->CalculateRemnantRadius(), m_SemiMajorAxis * AU_TO_RSOL) > 0) { // catch merger in CE here, do not update stars + m_MassTransferTrackerHistory = MT_TRACKING::MERGER; + m_Flags.stellarMerger = true; + } + + if (!m_Flags.stellarMerger) { // stellar merger? + // no - continue evolution STELLAR_TYPE stellarType1 = m_Star1->StellarType(); // star 1 stellar type before resolving envelope loss STELLAR_TYPE stellarType2 = m_Star2->StellarType(); // star 2 stellar type before resolving envelope loss @@ -1670,30 +1675,26 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { if (m_Star1->StellarType() != stellarType1 || m_Star2->StellarType() != stellarType2) { // stellar type change? (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::STELLAR_TYPE_CHANGE_DURING_CEE); // yes - print (log) detailed output } - } - if (utils::Compare(m_SemiMajorAxis, 0.0) <= 0 || utils::Compare(m_Star1->Radius() + m_Star2->Radius(), m_SemiMajorAxis * AU_TO_RSOL) > 0) { - m_Flags.stellarMerger = true; - } - - // if stars are evolving as CHE stars, update their rotational frequency under the assumption of tidal locking if tides are not enabled - if (OPTIONS->TidesPrescription() == TIDES_PRESCRIPTION::NONE) { - double omega = OrbitalAngularVelocity(); // orbital angular velocity - if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star1->SetOmega(omega); - if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star2->SetOmega(omega); - } + // if stars are evolving as CHE stars, update their rotational frequency under the assumption of tidal locking if tides are not enabled + if (!m_Flags.stellarMerger && OPTIONS->TidesPrescription() == TIDES_PRESCRIPTION::NONE) { + double omega = OrbitalAngularVelocity(); // orbital angular velocity + if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star1->SetOmega(omega); + if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star2->SetOmega(omega); + } - - m_Star1->SetPostCEEValues(); // squirrel away post CEE stellar values for star 1 - m_Star2->SetPostCEEValues(); // squirrel away post CEE stellar values for star 2 - SetPostCEEValues(m_SemiMajorAxis * AU_TO_RSOL, m_Eccentricity, rRLdfin1Rsol, rRLdfin2Rsol); // squirrel away post CEE binary values (checks for post-CE RLOF, so should be done at end) + m_Star1->SetPostCEEValues(); // squirrel away post CEE stellar values for star 1 + m_Star2->SetPostCEEValues(); // squirrel away post CEE stellar values for star 2 + SetPostCEEValues(m_SemiMajorAxis * AU_TO_RSOL, m_Eccentricity, rRLdfin1Rsol, rRLdfin2Rsol); // squirrel away post CEE binary values (checks for post-CE RLOF, so should be done at end) - if (m_RLOFDetails.immediateRLOFPostCEE == true && !OPTIONS->AllowImmediateRLOFpostCEToSurviveCommonEnvelope()) { // is there immediate post-CE RLOF which is not allowed? + if (m_RLOFDetails.immediateRLOFPostCEE == true && !OPTIONS->AllowImmediateRLOFpostCEToSurviveCommonEnvelope()) {// is there immediate post-CE RLOF which is not allowed? m_MassTransferTrackerHistory = MT_TRACKING::MERGER; m_Flags.stellarMerger = true; + } } - + (void)PrintCommonEnvelope(); // print (log) common envelope details + } @@ -1717,9 +1718,8 @@ void BaseBinaryStar::ResolveMainSequenceMerger() { double tau1 = m_Star1->Tau(); double tau2 = m_Star2->Tau(); - // /*ILYA*/ temporary solution, should use TAMS core mass - double TAMSCoreMass1 = 0.3 * mass1; - double TAMSCoreMass2 = 0.3 * mass2; + double TAMSCoreMass1 = m_Star1->TAMSCoreMass(); + double TAMSCoreMass2 = m_Star2->TAMSCoreMass(); double q = std::min(mass1 / mass2, mass2 / mass1); double phi = 0.3 * q / (1.0 + q) / (1.0 + q); // fraction of mass lost in merger, Wang+ 2022, https://www.nature.com/articles/s41550-021-01597-5 @@ -2891,24 +2891,27 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { } } - if ((m_Star1->IsSNevent() || m_Star2->IsSNevent())) { - EvaluateSupernovae(); // evaluate supernovae (both stars) if mass changes are responsible for a supernova - (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_SN); // print (log) detailed output - if (HasOneOf({ STELLAR_TYPE::NEUTRON_STAR })) { - (void)PrintPulsarEvolutionParameters(PULSAR_RECORD_TYPE::POST_SN); // print (log) pulsar evolution parameters + if (!StellarMerger()) { // stellar merger? + // no - continue evolution + if ((m_Star1->IsSNevent() || m_Star2->IsSNevent())) { + EvaluateSupernovae(); // evaluate supernovae (both stars) if mass changes are responsible for a supernova + (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_SN); // print (log) detailed output + if (HasOneOf({ STELLAR_TYPE::NEUTRON_STAR })) { + (void)PrintPulsarEvolutionParameters(PULSAR_RECORD_TYPE::POST_SN); // print (log) pulsar evolution parameters + } } - } - CalculateEnergyAndAngularMomentum(); // perform energy and angular momentum calculations + CalculateEnergyAndAngularMomentum(); // perform energy and angular momentum calculations - ProcessTides(p_Dt); // process tides if required + ProcessTides(p_Dt); // process tides if required - // assign new values to "previous" values, for following timestep - m_EccentricityPrev = m_Eccentricity; - m_SemiMajorAxisPrev = m_SemiMajorAxis; + // assign new values to "previous" values, for following timestep + m_EccentricityPrev = m_Eccentricity; + m_SemiMajorAxisPrev = m_SemiMajorAxis; - m_Star1->UpdateMagneticFieldAndSpin(m_CEDetails.CEEnow, m_Dt * MYR_TO_YEAR * SECONDS_IN_YEAR, EPSILON_PULSAR); // update pulsar parameters for star1 - m_Star2->UpdateMagneticFieldAndSpin(m_CEDetails.CEEnow, m_Dt * MYR_TO_YEAR * SECONDS_IN_YEAR, EPSILON_PULSAR); // update pulsar parameters for star2 + m_Star1->UpdateMagneticFieldAndSpin(m_CEDetails.CEEnow, m_Dt * MYR_TO_YEAR * SECONDS_IN_YEAR, EPSILON_PULSAR); // update pulsar parameters for star1 + m_Star2->UpdateMagneticFieldAndSpin(m_CEDetails.CEEnow, m_Dt * MYR_TO_YEAR * SECONDS_IN_YEAR, EPSILON_PULSAR); // update pulsar parameters for star2 + } } @@ -3081,8 +3084,12 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() { if (StellarMerger() && !HasOneOf({ STELLAR_TYPE::MASSLESS_REMNANT })) { // have stars merged without merger already being resolved? if (m_Star1->IsOneOf(MAIN_SEQUENCE) && m_Star2->IsOneOf(MAIN_SEQUENCE) && OPTIONS->EvolveMainSequenceMergers()) // yes - both MS and evolving MS merger products? ResolveMainSequenceMerger(); // yes - handle main sequence mergers gracefully; no need to change evolution status - else + else { + // make both stars massless remnants if merging during CE, so this is recorded in the Switch log; eventually, will want to implement a more careful prescription for the merger product, perhaps allowing further evolution of the merger product + m_Star1->SwitchTo(STELLAR_TYPE::MASSLESS_REMNANT); + m_Star2->SwitchTo(STELLAR_TYPE::MASSLESS_REMNANT); evolutionStatus = EVOLUTION_STATUS::STELLAR_MERGER; // no - for now, stop evolution + } } else if (HasStarsTouching()) { // binary components touching? (should usually be avoided as MT or CE or merger should happen prior to this) evolutionStatus = EVOLUTION_STATUS::STARS_TOUCHING; // yes - stop evolution @@ -3140,45 +3147,48 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() { (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::PRE_STELLAR_TIMESTEP); // print (log) detailed output - error = EvolveOneTimestep(dt); // evolve the binary system one timestep - if (error != ERROR::NONE) { // SSE error for either constituent star? - evolutionStatus = EVOLUTION_STATUS::SSE_ERROR; // yes - stop evolution - } - else { // continue evolution + if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // continue evolution? + // yes + error = EvolveOneTimestep(dt); // evolve the binary system one timestep + if (error != ERROR::NONE) { // SSE error for either constituent star? + evolutionStatus = EVOLUTION_STATUS::SSE_ERROR; // yes - stop evolution + } + else { // continue evolution - (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::TIMESTEP_COMPLETED); // print (log) detailed output: this is after all changes made in the timestep + (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::TIMESTEP_COMPLETED); // print (log) detailed output: this is after all changes made in the timestep - if (stepNum >= OPTIONS->MaxNumberOfTimestepIterations()) evolutionStatus = EVOLUTION_STATUS::STEPS_UP; // number of timesteps for evolution exceeds maximum - else if (evolutionStatus == EVOLUTION_STATUS::CONTINUE && usingProvidedTimesteps && stepNum >= timesteps.size()) { // using user-provided timesteps and all consumed - evolutionStatus = EVOLUTION_STATUS::TIMESTEPS_EXHAUSTED; // yes - set status - SHOW_WARN(ERROR::TIMESTEPS_EXHAUSTED); // show warning + if (stepNum >= OPTIONS->MaxNumberOfTimestepIterations()) evolutionStatus = EVOLUTION_STATUS::STEPS_UP; // number of timesteps for evolution exceeds maximum + else if (evolutionStatus == EVOLUTION_STATUS::CONTINUE && usingProvidedTimesteps && stepNum >= timesteps.size()) { // using user-provided timesteps and all consumed + evolutionStatus = EVOLUTION_STATUS::TIMESTEPS_EXHAUSTED; // yes - set status + SHOW_WARN(ERROR::TIMESTEPS_EXHAUSTED); // show warning + } } - } - if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // continue evolution? + if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // continue evolution? // yes - // if user selects to emit GWs, calculate the effects of radiation - // - note that this is placed before the call to ChooseTimestep() because when - // emitting GWs the timestep is a function of gravitational radiation - if (OPTIONS->EmitGravitationalRadiation()) CalculateGravitationalRadiation(); + // if user selects to emit GWs, calculate the effects of radiation + // - note that this is placed before the call to ChooseTimestep() because when + // emitting GWs the timestep is a function of gravitational radiation + if (OPTIONS->EmitGravitationalRadiation()) CalculateGravitationalRadiation(); - m_Star2->UpdatePreviousTimestepDuration(); // update stellar property for star2 - m_Star1->UpdatePreviousTimestepDuration(); // update stellar property for star1 + m_Star2->UpdatePreviousTimestepDuration(); // update stellar property for star2 + m_Star1->UpdatePreviousTimestepDuration(); // update stellar property for star1 - if (usingProvidedTimesteps) { // user-provided timesteps? - // select a timestep - // - don't quantise - // - don't apply timestep multiplier - // (we assume user wants the timesteps in the file) - // - // Open question: should we clamp this to NUCLEAR_MINIMUM_TIMESTEP? - dt = timesteps[stepNum]; - } - else { // no - not using user-provided timesteps - dt = ChooseTimestep(OPTIONS->TimestepMultiplier()); - } + if (usingProvidedTimesteps) { // user-provided timesteps? + // select a timestep + // - don't quantise + // - don't apply timestep multiplier + // (we assume user wants the timesteps in the file) + // + // Open question: should we clamp this to NUCLEAR_MINIMUM_TIMESTEP? + dt = timesteps[stepNum]; + } + else { // no - not using user-provided timesteps + dt = ChooseTimestep(OPTIONS->TimestepMultiplier()); + } - stepNum++; // increment stepNum + stepNum++; // increment stepNum + } } } diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index a550f7750..500fb76d3 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -2478,6 +2478,8 @@ double BaseStar::CalculateMassLossRateWolfRayetSanderVink2020(const double p_Mu) */ double BaseStar::CalculateMassLossRateWolfRayetTemperatureCorrectionSander2023(const double p_Mdot) const { + if (p_Mdot <= 0.0) return 0.0; // nothing to adjust + const double teffRef = 141.0E3; // reference effective temperature in Kelvin const double teffMin = 100.0E3; // minimum effective temperature in Kelvin to apply correction @@ -2493,6 +2495,8 @@ double BaseStar::CalculateMassLossRateWolfRayetTemperatureCorrectionSander2023(c logMdotCorrected = logMdotUncorrected; } + std::cout<StashSSESupernovaDetails(this, p_StellarType, p_RecordType); } + virtual double TAMSCoreMass() const { return 0.0; } // except MS stars + virtual void UpdateAfterMerger(double p_Mass, double p_HydrogenMass) { } // Default is NO-OP virtual void UpdateAgeAfterMassLoss() { } // Default is NO-OP diff --git a/src/GiantBranch.h b/src/GiantBranch.h index 93c432029..d07ec4fc5 100755 --- a/src/GiantBranch.h +++ b/src/GiantBranch.h @@ -20,9 +20,8 @@ class GiantBranch: virtual public BaseStar, public MainSequence { GiantBranch(){}; GiantBranch(const BaseStar &p_BaseStar) : BaseStar(p_BaseStar), MainSequence(p_BaseStar) {} - - virtual double CalculateRemnantRadius() const; - + + double CalculateRemnantRadius() const; protected: diff --git a/src/HG.cpp b/src/HG.cpp index 861de1f7f..5d028a53e 100755 --- a/src/HG.cpp +++ b/src/HG.cpp @@ -2,6 +2,7 @@ #include "HeMS.h" #include "HeWD.h" + /////////////////////////////////////////////////////////////////////////////////////// // // // COEFFICIENT AND CONSTANT CALCULATIONS ETC. // @@ -861,6 +862,7 @@ double HG::CalculateRadiusOnPhase(const double p_Mass, const double p_Tau, const } + /////////////////////////////////////////////////////////////////////////////////////// // // // MASS CALCULATIONS // diff --git a/src/HG.h b/src/HG.h index 1b7ea3e90..cf23a5fc6 100755 --- a/src/HG.h +++ b/src/HG.h @@ -11,52 +11,49 @@ #include "GiantBranch.h" -// JR: todo: revisit this one day - sometimes HG works better as GiantBranch, sometimes not... -// Right now it is GiantBranch - figure out which is more appropriate - class BaseStar; class GiantBranch; class HG: virtual public BaseStar, public GiantBranch { - + public: - + HG() { m_StellarType = STELLAR_TYPE::HERTZSPRUNG_GAP; }; HG(const BaseStar &p_BaseStar, const bool p_Initialise = true) : BaseStar(p_BaseStar), GiantBranch(p_BaseStar) { - m_StellarType = STELLAR_TYPE::HERTZSPRUNG_GAP; // Set stellar type + m_StellarType = STELLAR_TYPE::HERTZSPRUNG_GAP; // Set stellar type if (p_Initialise) Initialise(); // Initialise if required } - + HG* Clone(const OBJECT_PERSISTENCE p_Persistence, const bool p_Initialise = true) { - HG* clone = new HG(*this, p_Initialise); - clone->SetPersistence(p_Persistence); - return clone; + HG* clone = new HG(*this, p_Initialise); + clone->SetPersistence(p_Persistence); + return clone; } - + static HG* Clone(HG& p_Star, const OBJECT_PERSISTENCE p_Persistence, const bool p_Initialise = true) { - HG* clone = new HG(p_Star, p_Initialise); - clone->SetPersistence(p_Persistence); - return clone; + HG* clone = new HG(p_Star, p_Initialise); + clone->SetPersistence(p_Persistence); + return clone; } - + MT_CASE DetermineMassTransferTypeAsDonor() const { return MT_CASE::B; } // Always case B - - + + protected: - + void Initialise() { - + m_Tau = 0.0; // Start of phase - + // update stellar properties at start of HG phase (since core definition changes) CalculateGBParams(); CalculateTimescales(); - // Initialise timescales + // Initialise timescales m_Age = m_Timescales[static_cast(TIMESCALE::tMS)]; // Set age appropriately - + // update effective "initial" mass (m_Mass0) so that the core mass is at least equal to the minimum core mass but no more than total mass - // (only relevant if RetainCoreMassDuringCaseAMassTransfer()) + // (only relevant if RetainCoreMassDuringCaseAMassTransfer()) if (utils::Compare(CalculateCoreMassOnPhase(m_Mass0, m_Age), std::min(m_Mass, MinimumCoreMass())) < 0) { double desiredCoreMass = std::min(m_Mass, MinimumCoreMass()); // desired core mass m_Mass0 = Mass0ToMatchDesiredCoreMass(this, desiredCoreMass); // use root finder to find new core mass estimate @@ -69,56 +66,56 @@ class HG: virtual public BaseStar, public GiantBranch { } EvolveOnPhase(0.0); } - - + + // member functions - alphabetically double CalculateCOCoreMassAtPhaseEnd() const { return 0.0; } // McCO(HG) = 0.0 double CalculateCOCoreMassOnPhase() const { return 0.0; } // McCO(HG) = 0.0 - + double CalculateCoreMassAt2ndDredgeUp(const DBL_VECTOR &p_GBParams) { return p_GBParams[static_cast(GBP::McDU)]; } // NO-OP double CalculateCoreMassAtPhaseEnd(const double p_Mass) const; double CalculateCoreMassAtPhaseEnd() const { return CalculateCoreMassAtPhaseEnd(m_Mass0); } // Use class member variables double CalculateCoreMassOnPhase(const double p_Mass, const double p_Time) const; double CalculateCoreMassOnPhase() const { return CalculateCoreMassOnPhase(m_Mass0, m_Age); } // Use class member variables double CalculateCoreMassOnPhaseIgnoringPreviousCoreMass(const double p_Mass, const double p_Time) const; // Ignore previous core mass constraint when computing expected core mass - + double CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const; double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return 0.25; } // As coded in BSE. Using the inverse owing to how qCrit is defined in COMPAS. See Hurley et al. 2002 sect. 2.6.1 for additional details. - + double CalculateHeCoreMassAtPhaseEnd() const { return m_CoreMass; } // McHe(HG) = Core Mass double CalculateHeCoreMassOnPhase() const { return m_CoreMass; } // McHe(HG) = Core Mass - double CalculateHeliumAbundanceCoreAtPhaseEnd() const { return 1.0 - m_Metallicity; } - double CalculateHeliumAbundanceCoreOnPhase() const { return 1.0 - m_Metallicity; } // Use class member variables + double CalculateHeliumAbundanceCoreAtPhaseEnd() const { return 1.0 - m_Metallicity; } + double CalculateHeliumAbundanceCoreOnPhase() const { return 1.0 - m_Metallicity; } // Use class member variables double CalculateHeliumAbundanceSurfaceAtPhaseEnd() const { return CalculateHeliumAbundanceSurfaceOnPhase(); } - double CalculateHeliumAbundanceSurfaceOnPhase() const { return m_InitialHeliumAbundance; } // Use class member variables + double CalculateHeliumAbundanceSurfaceOnPhase() const { return m_InitialHeliumAbundance; } // Use class member variables - double CalculateHydrogenAbundanceCoreAtPhaseEnd() const { return CalculateHydrogenAbundanceCoreOnPhase(); } - double CalculateHydrogenAbundanceCoreOnPhase(const double p_Tau) const; - double CalculateHydrogenAbundanceCoreOnPhase() const { return 0.0; } // Star has exhausted hydrogen in its core + double CalculateHydrogenAbundanceCoreAtPhaseEnd() const { return CalculateHydrogenAbundanceCoreOnPhase(); } + double CalculateHydrogenAbundanceCoreOnPhase(const double p_Tau) const; + double CalculateHydrogenAbundanceCoreOnPhase() const { return 0.0; } // Star has exhausted hydrogen in its core - double CalculateHydrogenAbundanceSurfaceAtPhaseEnd() const { return CalculateHydrogenAbundanceSurfaceOnPhase(); } + double CalculateHydrogenAbundanceSurfaceAtPhaseEnd() const { return CalculateHydrogenAbundanceSurfaceOnPhase(); } double CalculateHydrogenAbundanceSurfaceOnPhase() const { return m_InitialHydrogenAbundance; } // Use class member variables - - + + double CalculateLambdaDewi() const; double CalculateLambdaNanjingStarTrack(const double p_Mass, const double p_Metallicity) const; double CalculateLambdaNanjingEnhanced(const int p_MassIndex, const STELLAR_POPULATION p_StellarPop) const; - + double CalculateLuminosityAtPhaseEnd(const double p_Mass) const; double CalculateLuminosityAtPhaseEnd() const { return CalculateLuminosityAtPhaseEnd(m_Mass0);} // Use class member variables double CalculateLuminosityOnPhase(const double p_Age, const double p_Mass) const; double CalculateLuminosityOnPhase() const { return CalculateLuminosityOnPhase(m_Age, m_Mass0); } // Use class member variables - + double CalculateMassTransferRejuvenationFactor() { return 1.0; } - + double CalculateRadiusAtPhaseEnd(const double p_Mass) const; double CalculateRadiusAtPhaseEnd() const { return CalculateRadiusAtPhaseEnd(m_Mass); } // Use class member variables double CalculateRadiusOnPhase(const double p_Mass, const double p_Tau, const double p_RZAMS) const; double CalculateRadiusOnPhase() const { return CalculateRadiusOnPhase(m_Mass0, m_Tau, m_RZAMS0); } // Use class member variables - + double CalculateRho(const double p_Mass) const; double CalculateTauAtPhaseEnd() const { return 1.0; } // tau = 1.0 at end of HG @@ -141,6 +138,7 @@ class HG: virtual public BaseStar, public GiantBranch { bool ShouldEvolveOnPhase() const { return (utils::Compare(m_Age, m_Timescales[static_cast(TIMESCALE::tBGB)]) < 0); } // Evolve on HG phase if age < Base Giant Branch timescale bool ShouldSkipPhase() const { return false; } // Never skip HG phase + double TAMSCoreMass() const { return 0.0; } void UpdateAfterMerger(double p_Mass, double p_HydrogenMass) { } // Nothing to do for stars beyond the Main Sequence for now void UpdateAgeAfterMassLoss(); // Per Hurley et al. 2000, section 7.1 diff --git a/src/HeHG.cpp b/src/HeHG.cpp index 56d411af7..34bfbd485 100755 --- a/src/HeHG.cpp +++ b/src/HeHG.cpp @@ -1,6 +1,7 @@ #include "HeHG.h" #include "HeGB.h" #include "COWD.h" +#include "HeWD.h" /* @@ -234,6 +235,21 @@ double HeHG::CalculatePerturbationMu() const { return std::max(5.0 * ((McMax - m_CoreMass) / McMax), 0.0); //return non-negative value to avoid round-off issues } +/* + * Calculate radius of the remnant the star would become if it lost all of its + * envelope immediately (i.e. M = Mc, coreMass) + * + * Hurley et al. 2000, end of section 6 + * + * + * double CalculateRemnantRadius() + * + * @return Radius of remnant core in Rsol + */ +double HeHG::CalculateRemnantRadius() const { + return HeWD::CalculateRadiusOnPhase_Static(m_CoreMass); +} + /////////////////////////////////////////////////////////////////////////////////////// // // diff --git a/src/HeHG.h b/src/HeHG.h index 038dea565..7956ef76e 100755 --- a/src/HeHG.h +++ b/src/HeHG.h @@ -40,6 +40,7 @@ class HeHG: virtual public BaseStar, public HeMS { double CalculateConvectiveCoreRadius () const { return std::min(5.0 * CalculateRemnantRadius(), m_Radius); } // Last paragraph of section 6 of Hurley+ 2000 static void CalculateGBParams_Static(const double p_Mass0, const double p_Mass, const double p_LogMetallicityXi, const DBL_VECTOR &p_MassCutoffs, const DBL_VECTOR &p_AnCoefficients, const DBL_VECTOR &p_BnCoefficients, DBL_VECTOR &p_GBParams); + double CalculateRemnantRadius() const; protected: diff --git a/src/HeMS.h b/src/HeMS.h index 00770d6fa..5b38a32e1 100644 --- a/src/HeMS.h +++ b/src/HeMS.h @@ -52,6 +52,8 @@ class HeMS: virtual public BaseStar, public TPAGB { static DBL_DBL CalculateRadiusAtPhaseEnd_Static(const double p_Mass, const double p_Luminosity); static double CalculateRadiusAtZAMS_Static(const double p_Mass); static double CalculateRadiusOnPhase_Static(const double p_Mass, const double p_Tau); + + double CalculateRemnantRadius() const { return Radius(); } MT_CASE DetermineMassTransferTypeAsDonor() const { return MT_CASE::OTHER; } // Not A, B, C, or NONE diff --git a/src/MainSequence.cpp b/src/MainSequence.cpp index fe69f6f3f..04b7657fe 100644 --- a/src/MainSequence.cpp +++ b/src/MainSequence.cpp @@ -492,10 +492,11 @@ double MainSequence::CalculateRadiusOnPhase(const double p_Mass, const double p_ double tau1 = std::min(1.0, (p_Time / tHook)); // Hurley et al. 2000, eq 14 double tau2 = std::max(0.0, std::min(1.0, (p_Time - ((1.0 - epsilon) * tHook)) / (epsilon * tHook))); // Hurley et al. 2000, eq 15 - // pow() is slow - use multipliaction where it makes sense + // pow() is slow - use multiplication where it makes sense double tau_3 = tau * tau * tau; - double tau_10 = tau_3 * tau_3 * tau_3 * tau; - double tau_40 = tau_10 * tau_10 * tau_10 * tau_10; + double tau_10 = tau < FLOAT_TOLERANCE_ABSOLUTE ? 0.0: tau_3 * tau_3 * tau_3 * tau; // direct comparison, to avoid underflow + double tau_40 = tau_10 < FLOAT_TOLERANCE_ABSOLUTE ? 0.0: tau_10 * tau_10 * tau_10 * tau_10; // direct comparison, to avoid underflow + double tau1_3 = tau1 * tau1 * tau1; double tau2_3 = tau2 * tau2 * tau2; @@ -636,20 +637,7 @@ double MainSequence::CalculateConvectiveCoreRadius() const { * @return Mass of convective core in Msol */ double MainSequence::CalculateConvectiveCoreMass() const { - - // We need TAMSCoreMass, which is just the core mass at the start of the HG phase. - // Since we are on the main sequence here, we can clone this object as an HG object - // and, as long as it is initialised (to correctly set Tau to 0.0 on the HG phase), - // we can query the cloned object for its core mass. - // - // The clone should not evolve, and so should not log anything, but to be sure the - // clone does not participate in logging, we set its persistence to EPHEMERAL. - - HG *clone = HG::Clone(static_cast(const_cast(*this)), OBJECT_PERSISTENCE::EPHEMERAL); - double TAMSCoreMass = clone->CoreMass(); // get core mass from clone - delete clone; clone = nullptr; // return the memory allocated for the clone - - double finalConvectiveCoreMass = TAMSCoreMass; + double finalConvectiveCoreMass = TAMSCoreMass(); // core mass at TAMS double initialConvectiveCoreMass = finalConvectiveCoreMass / 0.6; return (initialConvectiveCoreMass - m_Tau * (initialConvectiveCoreMass - finalConvectiveCoreMass)); } @@ -849,23 +837,34 @@ STELLAR_TYPE MainSequence::ResolveEnvelopeLoss(bool p_Force) { */ void MainSequence::UpdateMinimumCoreMass() { if (OPTIONS->RetainCoreMassDuringCaseAMassTransfer()) { - - // We need TAMSCoreMass, which is just the core mass at the start of the HG phase. - // Since we are on the main sequence here, we can clone this object as an HG object - // and, as long as it is initialised (to correctly set Tau to 0.0 on the HG phase), - // we can query the cloned object for its core mass. - // - // The clone should not evolve, and so should not log anything, but to be sure the - // clone does not participate in logging, we set its persistence to EPHEMERAL. - - HG *clone = HG::Clone(static_cast(const_cast(*this)), OBJECT_PERSISTENCE::EPHEMERAL); - double TAMSCoreMass = clone->CoreMass(); // get core mass from clone - delete clone; clone = nullptr; // return the memory allocated for the clone - - m_MinimumCoreMass = std::max(m_MinimumCoreMass, CalculateTauOnPhase() * TAMSCoreMass); // update minimum core mass + m_MinimumCoreMass = std::max(m_MinimumCoreMass, CalculateTauOnPhase() * TAMSCoreMass()); // update minimum core mass } } +/* + * Return the expected core mass at terminal age main sequence, i.e., at the start of the HG phase + * + * double TAMSCoreMass() const + * + * + * @return TAMS core Mass (Msol) + * + */ +double MainSequence::TAMSCoreMass() const { + // Since we are on the main sequence here, we can clone this object as an HG object + // and, as long as it is initialised (to correctly set Tau to 0.0 on the HG phase), + // we can query the cloned object for its core mass. + // + // The clone should not evolve, and so should not log anything, but to be sure the + // clone does not participate in logging, we set its persistence to EPHEMERAL. + + HG *clone = HG::Clone(static_cast(const_cast(*this)), OBJECT_PERSISTENCE::EPHEMERAL); + double TAMSCoreMass = clone->CoreMass(); // get core mass from clone + delete clone; clone = nullptr; // return the memory allocated for the clone + + return TAMSCoreMass; +} + /* * Sets the mass and age of a merge product of two main sequence stars diff --git a/src/MainSequence.h b/src/MainSequence.h index a83b4cbe7..36aa11b2c 100644 --- a/src/MainSequence.h +++ b/src/MainSequence.h @@ -107,6 +107,8 @@ class MainSequence: virtual public BaseStar { STELLAR_TYPE ResolveEnvelopeLoss(bool p_Force = false); bool ShouldEvolveOnPhase() const { return (m_Age < m_Timescales[static_cast(TIMESCALE::tMS)]); } // Evolve on MS phase if age in MS timescale + + double TAMSCoreMass() const; void UpdateInitialMass() { m_Mass0 = m_Mass; } // Per Hurley et al. 2000, section 7.1 diff --git a/src/Remnants.h b/src/Remnants.h index 762c154bf..587616564 100644 --- a/src/Remnants.h +++ b/src/Remnants.h @@ -24,6 +24,9 @@ class Remnants: virtual public BaseStar, public HeGB { // member functions static double CalculateRemnantMass_Static(const double p_COCoreMass) { return utils::Compare(p_COCoreMass, 7.0)< 0 ? 1.17 + (0.09 * p_COCoreMass) : p_COCoreMass; } // Hurley et al., eq 92; Hurley says that a >7.0 CO core leads to a collapse into a BH, but is ambiguous about the BH mass in this case -- we will assume it's just the CO core + + double CalculateRemnantRadius() const { return Radius(); } + protected: diff --git a/src/Star.h b/src/Star.h index 2b34e090e..c486cc1a9 100755 --- a/src/Star.h +++ b/src/Star.h @@ -197,6 +197,8 @@ class Star { double CalculateRadiusOnPhaseTau(const double p_Mass, const double p_Tau) const { return m_Star->CalculateRadiusOnPhaseTau(p_Mass, p_Tau); } + double CalculateRemnantRadius() { return m_Star->CalculateRemnantRadius(); } + void CalculateSNAnomalies(const double p_Eccentricity) { m_Star->CalculateSNAnomalies(p_Eccentricity); } double CalculateSNKickMagnitude(const double p_RemnantMass, @@ -260,6 +262,8 @@ class Star { STELLAR_TYPE SwitchTo(const STELLAR_TYPE p_StellarType, bool p_SetInitialType = false); + + double TAMSCoreMass() const { return m_Star->TAMSCoreMass(); } void UpdateAfterMerger(double p_Mass, double p_HydrogenMass) { m_Star->UpdateAfterMerger(p_Mass, p_HydrogenMass); } void UpdateAgeAfterMassLoss() { m_Star->UpdateAgeAfterMassLoss(); } diff --git a/src/changelog.h b/src/changelog.h index 3743c765a..b85e946e3 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1344,7 +1344,7 @@ // - Changed the prescription for Tonset in the Picker+ models to take advantage of improved metallicity-dependent fits // 03.05.02 IM - Oct 10, 2024 - Enhancement, defect repair: // - Reverted IsCCSN() to include USSN following a change in 3.00.00 that inadvertently led to no binary orbit updates following USSNe -// - Include a call to EvolveBinary(0.0) on initialisation of evolved stellar types; this ensures that Switch logs include consistent data for the new stellar type +// - Include a call to EvolveOnPhase(0.0) on initialisation of evolved stellar types; this ensures that Switch logs include consistent data for the new stellar type // 03.06.00 IM - Oct 14, 2024 - Enhancement, code cleanup: // - Incorporating the Maltsev+ (2024) prescription for supernova remnant masses // - Minor fixes, including in fallback fraction for Schneider SN prescription, documentation @@ -1404,6 +1404,11 @@ // 03.09.02 RTW - Nov 27, 2024 - Defect repair, enhancement: // - Fixed bugs in vector3d related to indexing and rotation // - Added tweak for circular systems at first SN, to fix the x-axis along the separation vector -const std::string VERSION_STRING = "03.09.02"; +// 03.09.03 IM - Nov 28, 2024 - Enhancement, defect repair: +// - Delay changing stellar types until after checking for whether remnant cores would touch in a common enevelope, use core radii instead (partial fix to #1286) +// - Define a new function MainSequence::TAMSCoreMass(); use it for determining the amount of He in a star during MS mergers +// - Switch both stars to Massless remnants during a CE merger, resolve #1265 +// - Minor fixes, including to #1255, #1258 +const std::string VERSION_STRING = "03.09.03"; # endif // __changelog_h__