From 91164505222023c45c78bda75b9021957a06b201 Mon Sep 17 00:00:00 2001 From: Jeff Riley Date: Fri, 5 Jan 2024 14:35:04 +1100 Subject: [PATCH 1/4] BSE vs SSE investigation and changes --- src/BaseBinaryStar.cpp | 84 +++++++++++++++++++++++++++--- src/BaseBinaryStar.h | 3 +- src/BaseStar.cpp | 86 ++++++++++++++++--------------- src/BaseStar.h | 7 ++- src/CHeB.cpp | 18 +++---- src/EAGB.cpp | 22 ++++---- src/GiantBranch.cpp | 11 ++-- src/HG.cpp | 2 +- src/HG.h | 10 ++-- src/HeHG.cpp | 14 ++--- src/HeMS.cpp | 4 +- src/Log.h | 6 +-- src/MainSequence.cpp | 13 ++--- src/Options.cpp | 18 +++++-- src/Options.h | 11 +++- src/Remnants.h | 2 +- src/Star.cpp | 94 +++++++++++++++++++++++++--------- src/Star.h | 4 +- src/TPAGB.cpp | 16 +++--- src/changelog.h | 3 +- src/constants.h | 113 +++++++++++++++++++++++++---------------- src/utils.cpp | 105 ++++++++++++++++++++++++++++++++++++++ src/utils.h | 2 + src/vector3d.h | 1 - src/yaml.h | 1 + 25 files changed, 457 insertions(+), 193 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 4713e990a..4ca3f701a 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -2407,10 +2407,37 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() { SAY("RandomSeed " << m_RandomSeed); } + bool usingProvidedTimesteps = false; // using user-provided timesteps? + DBL_VECTOR timesteps; + if (!OPTIONS->TimestepsFileName().empty()) { // have timesteps filename? + // yes + ERROR error; + std::tie(error, timesteps) = utils::ReadTimesteps(OPTIONS->TimestepsFileName()); // read timesteps from file + if (error != ERROR::NONE) { // ok? + evolutionStatus = EVOLUTION_STATUS::NO_TIMESTEPS; // no - set status + SHOW_WARN(error); // show warning + } + else usingProvidedTimesteps = true; // have user-provided timesteps + } + if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // continue evolution // evolve the current binary up to the maximum evolution time (and number of steps) - double dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep()) / 1000.0; // initialise the timestep - int stepNum = 1; // initialise step number + + double dt; // timestep + if (usingProvidedTimesteps) { // user-provided timesteps? + // get new timestep + // - don't quantise + // - don't apply timestep multiplier + // (we assume user wants the timesteps in the file) + dt = timesteps[0]; + } + else { // no - not using user-provided timesteps + dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep()) * OPTIONS->TimestepMultiplier() / 1000.0; // calculate timestep - make first step small + dt = std::max(std::round(dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, NUCLEAR_MINIMUM_TIMESTEP); // quantised + } + + long unsigned int stepNum = 1; // initialise step number + while (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // perform binary evolution - iterate over timesteps until told to stop EvolveOneTimestep(dt); // evolve the binary system one timestep @@ -2419,6 +2446,11 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() { if (m_Error != ERROR::NONE) { // SSE error for either constituent star? evolutionStatus = EVOLUTION_STATUS::SSE_ERROR; // yes - stop evolution } + + // check for reasons to not continue evolution + else if (!OPTIONS->EvolveDoubleWhiteDwarfs() && IsWDandWD()) { // double WD and their evolution is not enabled? + evolutionStatus = EVOLUTION_STATUS::WD_WD; // yes - do not evolve double WD systems + } else if (HasOneOf({ STELLAR_TYPE::MASSLESS_REMNANT })) { // at least one massless remnant? evolutionStatus = EVOLUTION_STATUS::MASSLESS_REMNANT; // yes - stop evolution } @@ -2444,7 +2476,7 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() { (void)PrintRLOFParameters(); // print (log) RLOF parameters - // check for problems + // check for problems/reasons to not continue evolution if (StellarMerger()) { // have stars merged? evolutionStatus = EVOLUTION_STATUS::STELLAR_MERGER; // for now, stop evolution } @@ -2486,20 +2518,37 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() { } } - // check for problems + // check whether to continue evolution if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // continue evolution? + + // check for problems if (m_Error != ERROR::NONE) { // error in binary evolution? evolutionStatus = EVOLUTION_STATUS::BINARY_ERROR; // yes - stop evolution } - else if (!OPTIONS->EvolveDoubleWhiteDwarfs() && IsWDandWD()) { // double WD and their evolution is not enabled? - evolutionStatus = EVOLUTION_STATUS::WD_WD; // yes - do not evolve double WD systems - } + + // check for reasons to stop evolution else if (IsDCO() && m_Time > (m_DCOFormationTime + m_TimeToCoalescence) && !IsUnbound()) { // evolution time exceeds DCO merger time? evolutionStatus = EVOLUTION_STATUS::DCO_MERGER_TIME; // yes - stop evolution } else if (m_Time > OPTIONS->MaxEvolutionTime()) { // evolution time exceeds maximum? evolutionStatus = EVOLUTION_STATUS::TIMES_UP; // yes - stop evolution } + else if (!OPTIONS->EvolveDoubleWhiteDwarfs() && IsWDandWD()) { // double WD and their evolution is not enabled? + evolutionStatus = EVOLUTION_STATUS::WD_WD; // yes - do not evolve double WD systems + } + else if (HasOneOf({ STELLAR_TYPE::MASSLESS_REMNANT })) { // at least one massless remnant? + evolutionStatus = EVOLUTION_STATUS::MASSLESS_REMNANT; // yes - stop evolution + } + else if (StellarMerger()) { // have stars merged? + evolutionStatus = EVOLUTION_STATUS::STELLAR_MERGER; // 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 + } + else if (IsUnbound() && !OPTIONS->EvolveUnboundSystems()) { // binary is unbound and we don't want unbound systems? + m_Unbound = true; // yes - set the unbound flag (should already be set) + evolutionStatus = EVOLUTION_STATUS::UNBOUND; // stop evolution + } } } } @@ -2507,16 +2556,35 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() { (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 (usingProvidedTimesteps && stepNum >= timesteps.size()) { + evolutionStatus = EVOLUTION_STATUS::TIMESTEPS_EXHAUSTED; // using user-provided timesteps and all consumed + SHOW_WARN(ERROR::TIMESTEPS_EXHAUSTED); // show warning + } if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // continue evolution? + if (usingProvidedTimesteps) { // user-provided timesteps + // get new timestep + // - don't quantise + // - don't apply timestep multiplier + // (we assume user wants the timesteps in the file) + dt = timesteps[stepNum]; + } + else { // not using user-provided timesteps + dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep()) * OPTIONS->TimestepMultiplier(); // calculate new timestep + dt = std::max(std::round(dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, NUCLEAR_MINIMUM_TIMESTEP); // quantised + } - dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep()) * OPTIONS->TimestepMultiplier(); // new timestep if ((m_Star1->IsOneOf({ STELLAR_TYPE::MASSLESS_REMNANT }) || m_Star2->IsOneOf({ STELLAR_TYPE::MASSLESS_REMNANT })) || dt < NUCLEAR_MINIMUM_TIMESTEP) { dt = NUCLEAR_MINIMUM_TIMESTEP; // but not less than minimum } stepNum++; // increment stepNum } } + + if (usingProvidedTimesteps && timesteps.size() > stepNum) { // all user-defined timesteps consumed? + evolutionStatus = EVOLUTION_STATUS::TIMESTEPS_NOT_CONSUMED; // no - set status + SHOW_WARN(ERROR::TIMESTEPS_NOT_CONSUMED); // show warning + } } m_EvolutionStatus = evolutionStatus; diff --git a/src/BaseBinaryStar.h b/src/BaseBinaryStar.h index 63cd97bc8..868f4d8b7 100644 --- a/src/BaseBinaryStar.h +++ b/src/BaseBinaryStar.h @@ -495,8 +495,7 @@ class BaseBinaryStar { return LOGGING->LogBSESystemParameters(this, p_RecordType); } - bool PrintDetailedOutput(const long int p_Id, - const BSE_DETAILED_RECORD_TYPE p_RecordType) const { + bool PrintDetailedOutput(const long int p_Id, const BSE_DETAILED_RECORD_TYPE p_RecordType) const { return OPTIONS->DetailedOutput() ? LOGGING->LogBSEDetailedOutput(this, p_Id, p_RecordType) : true; } diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index 21308dc3c..195028f5b 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -2587,7 +2587,8 @@ double BaseStar::CalculateMassLossRateBelczynski2010() { /* - * Calculate the mass loss rate according to the updated prescription. The structure is similar to the Vink wrapper (previous default), which should be called Belczynski. + * Calculate the mass loss rate according to the updated prescription. + * The structure is similar to the Vink wrapper (previous default), which should be called Belczynski. * Mass loss rates for hot, massive OB stars are given by Bjorklund et al. 2022 * Mass loss rates for helium rich Wolf--Rayet stars are given by Sander (not yet implemented) * Mass loss rates for red supergiants are given by Beasor and Davies (not yet implemented) @@ -2604,15 +2605,14 @@ double BaseStar::CalculateMassLossRateFlexible2023() { double LBVRate = CalculateMassLossRateLBV(OPTIONS->LuminousBlueVariablePrescription()); // start with LBV winds (can be, and is often, 0.0) double otherWindsRate = 0.0; - double teff = TSOL * m_Temperature; if (m_DominantMassLossRate != MASS_LOSS_TYPE::LBV || OPTIONS->LuminousBlueVariablePrescription() == LBV_PRESCRIPTION::HURLEY_ADD ) { // check whether we should add other winds to the LBV winds (always for HURLEY_ADD prescription, only if not in LBV regime for others) - - if ((utils::Compare(teff, RSG_MAXIMUM_TEMP) < 0) && (utils::Compare(m_MZAMS, 8.0) >= 0) && - IsOneOf(GIANTS)) { // RSG criteria, below 8kK, above 8Msol, and core helium burning giant(CHeB, FGB, EAGB, TPAGB) + if ((utils::Compare(teff, RSG_MAXIMUM_TEMP) < 0) && + (utils::Compare(m_MZAMS, 8.0) >= 0) && // JR: where does this 8.0 come from? We shouldn't have undocumented magic numbers in the code... If it's an accepted threshold, let's document it and give it a name. + IsOneOf(GIANTS)) { // RSG criteria, below 8kK, above 8Msol, and core helium burning giant(CHeB, FGB, EAGB, TPAGB) otherWindsRate = CalculateMassLossRateRSG(OPTIONS->RSGMassLoss()); m_DominantMassLossRate = MASS_LOSS_TYPE::RSG; } @@ -2621,18 +2621,16 @@ double BaseStar::CalculateMassLossRateFlexible2023() { } // change to Kelvin so it can be compared with values as stated in Vink prescription else if (utils::Compare(m_MZAMS, VERY_MASSIVE_MINIMUM_MASS) >= 0) { otherWindsRate = CalculateMassLossRateVMS(OPTIONS->VMSMassLoss()); - m_DominantMassLossRate = MASS_LOSS_TYPE::VMS; // massive MS, >100 Msol. Alternately could use Luminosity or Gamma and Mass threshold + m_DominantMassLossRate = MASS_LOSS_TYPE::VMS; // massive MS, >100 Msol. Alternatively could use Luminosity or Gamma and Mass threshold } - else { otherWindsRate = CalculateMassLossRateOB(OPTIONS->OBMassLoss()); m_DominantMassLossRate = MASS_LOSS_TYPE::OB; } if (utils::Compare(LBVRate, otherWindsRate) > 0) { - m_DominantMassLossRate = MASS_LOSS_TYPE::LBV; // set LBV dominant again in case Hurley or OB overwrote it + m_DominantMassLossRate = MASS_LOSS_TYPE::LBV; // set LBV dominant again in case Hurley or OB overwrote it } - } return LBVRate + otherWindsRate; @@ -2682,13 +2680,13 @@ double BaseStar::CalculateMassLossRate() { break; default: // unknown mass-loss prescription - SHOW_WARN(ERROR::UNKNOWN_MASS_LOSS_PRESCRIPTION, "Using HURLEY"); // show warning + SHOW_WARN(ERROR::UNKNOWN_MASS_LOSS_PRESCRIPTION, "DEFAULT!"); // show warning LBVRate = CalculateMassLossRateLBV(LBV_PRESCRIPTION::HURLEY_ADD); - otherWindsRate = CalculateMassLossRateHurley(); + otherWindsRate = CalculateMassLossRateHurley(); // use HURLEY if (utils::Compare(LBVRate, otherWindsRate) > 0) { m_DominantMassLossRate = MASS_LOSS_TYPE::LBV; } - mDot = LBVRate + otherWindsRate; // use HURLEY + mDot = LBVRate + otherWindsRate; } mDot = mDot * OPTIONS->OverallWindMassLossMultiplier(); // apply overall wind mass loss multiplier } @@ -2698,7 +2696,7 @@ double BaseStar::CalculateMassLossRate() { /* - * Calculate mass loss given mass loss rate - uses current timestep (m_Dt) + * Calculate mass loss given mass loss rate - uses timestep passed (p_Dt) * Returned mass loss is limited to MAXIMUM_MASS_LOSS_FRACTION (e.g., 1%) of current mass * * @@ -2727,7 +2725,10 @@ double BaseStar::CalculateMassLoss_Static(const double p_Mass, const double p_Md * Returns existing value for mass if mass loss not being used (program option) * * - * double CalculateMassLossValues() // JR: todo: pick a better name for this... + * double CalculateMassLossValues(const bool p_UpdateMDot, const bool p_UpdateMDt) // JR: todo: pick a better name for this... + * + * @param [IN] p_UpdateMDot flag to indicate whether the class member variable m_Mdot should be updated (default is false) + * @param [IN] p_UpdateMDt flag to indicate whether the class member variable m_Dt should be updated (default is false) * @return calculated mass (mSol) */ double BaseStar::CalculateMassLossValues(const bool p_UpdateMDot, const bool p_UpdateMDt) { @@ -2738,7 +2739,7 @@ double BaseStar::CalculateMassLossValues(const bool p_UpdateMDot, const bool p_U if (OPTIONS->UseMassLoss()) { // only if using mass loss (program option) - mDot = CalculateMassLossRate(); // calculate mass loss rate + mDot = CalculateMassLossRate(); // calculate mass loss rate double massLoss = CalculateMassLoss_Static(mass, mDot, dt); // calculate mass loss - limited to (mass * MAXIMUM_MASS_LOSS_FRACTION) if (OPTIONS->CheckPhotonTiringLimit()) { @@ -2753,6 +2754,7 @@ double BaseStar::CalculateMassLossValues(const bool p_UpdateMDot, const bool p_U } else { dt = massLoss / (mDot * 1.0E6); // new timestep to match limited mass loss + dt = std::round(dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM; mass -= massLoss; // new mass based on limited mass loss if (p_UpdateMDt) m_Dt = dt; // update class member variable if necessary @@ -2775,18 +2777,21 @@ double BaseStar::CalculateMassLossValues(const bool p_UpdateMDot, const bool p_U * - applies mass rejuvenation factor and calculates new age * * - * double ResolveMassLoss() + * double ResolveMassLoss(const bool p_UpdateMDt) + * + * @param [IN] p_UpdateMDt flag to indicate whether the class member variable m_Dt should be updated in BaseStar::CalculateMassLossValues() + * (default is true) */ -void BaseStar::ResolveMassLoss() { +void BaseStar::ResolveMassLoss(const bool p_UpdateMDt) { if (OPTIONS->UseMassLoss()) { - m_Mass = CalculateMassLossValues(true, true); // calculate new values assuming mass loss applied - - m_HeCoreMass = std::min(m_HeCoreMass, m_Mass); // update He mass if mass loss is happening from He stars - - m_COCoreMass = std::min(m_COCoreMass, m_Mass); // not expected, only a precaution to avoid inconsistencies - m_CoreMass = std::min(m_CoreMass, m_Mass); - + double mass = CalculateMassLossValues(true, p_UpdateMDt); // calculate new values assuming mass loss applied + + // JR: this is here (and the change to the line above) to keep attributes in sync BSE vs SSE + // Supernovae are caught in UpdateAttributesAndAgeOneTimestep() (hence the need to move the + // call to PrintStashedSupernovaDetails() in Star:EvolveOneTimestep()) + UpdateAttributesAndAgeOneTimestep(mass - m_Mass, 0.0, 0.0, false); + UpdateInitialMass(); // update effective initial mass (MS, HG & HeMS) UpdateAgeAfterMassLoss(); // update age (MS, HG & HeMS) ApplyMassTransferRejuvenationFactor(); // apply age rejuvenation factor @@ -3843,6 +3848,8 @@ double BaseStar::LimitTimestep(const double p_Dt) { if (utils::Compare(massLoss, 0.0) > 0) { // No change if no mass loss double dtWind = massLoss / (mDot * 1.0E6); // Calculate timestep to match (possibly limited) mass loss + dtWind = std::max(std::round(dtWind / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, NUCLEAR_MINIMUM_TIMESTEP); + dt = std::max(std::round(dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, NUCLEAR_MINIMUM_TIMESTEP); dt = min(dt, dtWind); // choose dt } } @@ -3868,12 +3875,12 @@ double BaseStar::CalculateTimestep() { // calls this functiom, the GBParams and Timescale functions // are called here - CalculateGBParams(); // calculate giant branch parameters - CalculateTimescales(); // calculate timescales + CalculateGBParams(); // calculate giant branch parameters + CalculateTimescales(); // calculate timescales - double dt = ChooseTimestep(m_Age); + double dt = std::max(std::round(ChooseTimestep(m_Age) / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, NUCLEAR_MINIMUM_TIMESTEP); - return max(LimitTimestep(dt), NUCLEAR_MINIMUM_TIMESTEP); // clamp timestep to minimum NUCLEAR_MINIMUM_TIMESTEP + return max(LimitTimestep(dt), NUCLEAR_MINIMUM_TIMESTEP); // clamp timestep to minimum NUCLEAR_MINIMUM_TIMESTEP } @@ -3944,7 +3951,7 @@ void BaseStar::UpdateAttributesAndAgeOneTimestepPreamble(const double p_DeltaMas * in which case no change is made to the star's age or the physical time attribute. * * - * Checks whether the star: + * Checks whether the star: (JR: the first two comments below are not true - they probably were once, but not now... check why, and fix the code or fix the comments) * - is a massless remnant (checked after applying p_DeltaMass and p_DeltaMass0, but before applying p_DeltaTime) * - has become a supernova (checked after applying p_DeltaMass and p_DeltaMass0, but before applying p_DeltaTime) * - should skip this phase for this timestep (checked after applying p_DeltaMass, p_DeltaMass0 and p_DeltaTime) @@ -3982,27 +3989,23 @@ STELLAR_TYPE BaseStar::UpdateAttributesAndAgeOneTimestep(const double p_DeltaMas stellarType = STELLAR_TYPE::MASSLESS_REMNANT; } else { - stellarType = ResolveSupernova(); // handle supernova - + stellarType = ResolveSupernova(); // handle supernova if (stellarType == m_StellarType) { // still on phase? - UpdateAttributesAndAgeOneTimestepPreamble(p_DeltaMass, p_DeltaMass0, p_DeltaTime); // apply mass changes and save current values if required + UpdateAttributesAndAgeOneTimestepPreamble(p_DeltaMass, p_DeltaMass0, p_DeltaTime); // apply mass changes and save current values if required if (p_ForceRecalculate || // force recalculate? utils::Compare(p_DeltaMass, 0.0) != 0 || // mass change? or... utils::Compare(p_DeltaMass0, 0.0) != 0 || // mass0 change? or... - utils::Compare(p_DeltaTime, 0.0) > 0) { // age/time advance? + p_DeltaTime > 0) { // age/time advance? (don't use utils::Compare() here) // yes - update attributes - AgeOneTimestepPreamble(p_DeltaTime); // advance dt, age, simulation time + if (p_DeltaTime > 0.0) AgeOneTimestepPreamble(p_DeltaTime); // advance dt, age, simulation time if necessary (don't use utils::Compare() here) - if (ShouldSkipPhase()) { // skip phase? - stellarType = ResolveSkippedPhase(); // phase skipped - do what's required - } + if (ShouldSkipPhase()) stellarType = ResolveSkippedPhase(); // skip phase if required else { // not skipped - execute phase - stellarType = EvolveOnPhase(); // evolve on phase - - if (stellarType == m_StellarType) { // still on phase? - stellarType = ResolveEndOfPhase(); // check for need to move off phase + if (p_DeltaTime > 0.0) {; // evolve on phase if necessary (don't use utils::Compare() here) + stellarType = EvolveOnPhase(); // evolve on phase + if (stellarType == m_StellarType) stellarType = ResolveEndOfPhase();// check for need to move off phase } } } @@ -4028,7 +4031,6 @@ void BaseStar::AgeOneTimestepPreamble(const double p_DeltaTime) { m_Age += p_DeltaTime; // advance age of star m_Dt = p_DeltaTime; // set timestep } - EvolveOneTimestepPreamble(); } diff --git a/src/BaseStar.h b/src/BaseStar.h index 665462be8..3b1483e39 100644 --- a/src/BaseStar.h +++ b/src/BaseStar.h @@ -236,7 +236,7 @@ class BaseStar { virtual STELLAR_TYPE ResolveEnvelopeLoss(bool p_NoCheck = false) { return m_StellarType; } - virtual void ResolveMassLoss(); + virtual void ResolveMassLoss(const bool p_UpdateMDt = true); void SetStellarTypePrev(const STELLAR_TYPE p_StellarTypePrev) { m_StellarTypePrev = p_StellarTypePrev; } @@ -264,8 +264,7 @@ class BaseStar { // printing functions - bool PrintDetailedOutput(const int p_Id, - const SSE_DETAILED_RECORD_TYPE p_RecordType = SSE_DETAILED_RECORD_TYPE::DEFAULT) const { + bool PrintDetailedOutput(const int p_Id, const SSE_DETAILED_RECORD_TYPE p_RecordType) const { return OPTIONS->DetailedOutput() ? LOGGING->LogSSEDetailedOutput(this, p_Id, p_RecordType) : true; // Write record to SSE Detailed Output log file } @@ -622,7 +621,7 @@ class BaseStar { virtual void SetSNHydrogenContent() { m_SupernovaDetails.isHydrogenPoor = false; } // Default is false - bool ShouldBeMasslessRemnant() const { return (m_Mass <= 0.0 || m_StellarType==STELLAR_TYPE::MASSLESS_REMNANT); } + bool ShouldBeMasslessRemnant() const { return (m_Mass <= 0.0 || m_StellarType == STELLAR_TYPE::MASSLESS_REMNANT); } virtual bool ShouldEvolveOnPhase() const { return true; } virtual bool ShouldSkipPhase() const { return false; } // Default is false diff --git a/src/CHeB.cpp b/src/CHeB.cpp index 94346478c..025f51b05 100755 --- a/src/CHeB.cpp +++ b/src/CHeB.cpp @@ -1432,11 +1432,6 @@ double CHeB::ChooseTimestep(const double p_Time) const { double dtk = 2.0E-3 * timescales(tHe); double dte = timescales(tHeI) + timescales(tHe) - p_Time; -// JR: todo: I took this next IF out - was it left in by mistake after some debugging? -// if (dtk < dte) { -// dtk = 2E-3 * m_Timescales[TS::tHe]; -// } - double timestep = std::max(std::min(dtk, dte), NUCLEAR_MINIMUM_TIMESTEP); return timestep; @@ -1471,12 +1466,13 @@ double CHeB::ChooseTimestep(const double p_Time) const { */ STELLAR_TYPE CHeB::ResolveEnvelopeLoss(bool p_NoCheck) { #define timescales(x) m_Timescales[static_cast(TIMESCALE::x)] // for convenience and readability - undefined at end of function - STELLAR_TYPE stellarType = m_StellarType; - if(ShouldEnvelopeBeExpelledByPulsations()) { m_EnvelopeJustExpelledByPulsations = true; } + if(ShouldEnvelopeBeExpelledByPulsations()) { + m_EnvelopeJustExpelledByPulsations = true; + } - if (p_NoCheck || utils::Compare(m_CoreMass, m_Mass) >= 0 || m_EnvelopeJustExpelledByPulsations ) { // Envelope loss + if (p_NoCheck || utils::Compare(m_CoreMass, m_Mass) >= 0 || m_EnvelopeJustExpelledByPulsations ) { // Envelope loss m_Mass = std::min(m_CoreMass, m_Mass); m_CoreMass = m_Mass; @@ -1488,15 +1484,15 @@ STELLAR_TYPE CHeB::ResolveEnvelopeLoss(bool p_NoCheck) { double tHeIPrime = timescales(tHeI); double tHePrime = timescales(tHe); - m_Tau = ((m_Age - tHeIPrime) / tHePrime); // Hurley et al. 2000, just after eq 81 - m_Age = m_Tau * HeMS::CalculateLifetimeOnPhase_Static(m_Mass0); // Hurley et al. 2000, eq 76 and following discussion + m_Tau = (m_Age - tHeIPrime) / tHePrime; // Hurley et al. 2000, just after eq 81 + m_Age = m_Tau * HeMS::CalculateLifetimeOnPhase_Static(m_Mass0); // Hurley et al. 2000, eq 76 and following discussion CalculateTimescales(m_Mass0, m_Timescales); CalculateGBParams(m_Mass0, m_GBParams); m_Luminosity = HeMS::CalculateLuminosityOnPhase_Static(m_Mass, m_Tau); m_Radius = HeMS::CalculateRadiusOnPhase_Static(m_Mass, m_Tau); - stellarType = STELLAR_TYPE::NAKED_HELIUM_STAR_MS; // will evolve to an evolved helium star + stellarType = STELLAR_TYPE::NAKED_HELIUM_STAR_MS; // will evolve to an evolved helium star } return stellarType; diff --git a/src/EAGB.cpp b/src/EAGB.cpp index 548717b47..44d9958d8 100755 --- a/src/EAGB.cpp +++ b/src/EAGB.cpp @@ -1000,17 +1000,17 @@ STELLAR_TYPE EAGB::ResolveEnvelopeLoss(bool p_NoCheck) { if(ShouldEnvelopeBeExpelledByPulsations()) { m_EnvelopeJustExpelledByPulsations = true; } - if (p_NoCheck || utils::Compare(m_HeCoreMass, m_Mass) >= 0 || m_EnvelopeJustExpelledByPulsations) { // Envelope lost, form an evolved naked helium giant + if (p_NoCheck || utils::Compare(m_HeCoreMass, m_Mass) >= 0 || m_EnvelopeJustExpelledByPulsations) { // Envelope lost, form an evolved naked helium giant - m_Mass = std::min(m_HeCoreMass, m_Mass); - m_HeCoreMass = m_Mass; - m_Mass0 = m_Mass; - m_CoreMass = m_COCoreMass; + m_Mass = std::min(m_HeCoreMass, m_Mass); + m_HeCoreMass = m_Mass; + m_Mass0 = m_Mass; + m_CoreMass = m_COCoreMass; - double p1 = gbParams(p) - 1.0; - double q1 = gbParams(q) - 1.0; - double p1_p = p1 / gbParams(p); - double q1_q = q1 / gbParams(q); + double p1 = gbParams(p) - 1.0; + double q1 = gbParams(q) - 1.0; + double p1_p = p1 / gbParams(p); + double q1_q = q1 / gbParams(q); timescales(tHeMS) = HeMS::CalculateLifetimeOnPhase_Static(m_Mass); // calculate common values @@ -1020,7 +1020,7 @@ STELLAR_TYPE EAGB::ResolveEnvelopeLoss(bool p_NoCheck) { timescales(tx_HeGB) = timescales(tinf1_HeGB) - (timescales(tinf1_HeGB) - timescales(tHeMS)) * PPOW((LTHe / gbParams(Lx)), p1_p); timescales(tinf2_HeGB) = timescales(tx_HeGB) + ((1.0 / (q1 * gbParams(AHe) * gbParams(B))) * PPOW((gbParams(B) / gbParams(Lx)), q1_q)); - m_Age = HeGB::CalculateAgeOnPhase_Static(m_Mass, m_COCoreMass, timescales(tHeMS), m_GBParams); + m_Age = HeGB::CalculateAgeOnPhase_Static(m_Mass, m_COCoreMass, timescales(tHeMS), m_GBParams); HeHG::CalculateGBParams_Static(m_Mass0, m_Mass, LogMetallicityXi(), m_MassCutoffs, m_AnCoefficients, m_BnCoefficients, m_GBParams); // IM: order of type change and parameter updates to be revisited (e.g., why not just CalculateGBParams(m_Mass0, m_GBParams)?) JR: static function has no access to class variables m_Luminosity = HeGB::CalculateLuminosityOnPhase_Static(m_COCoreMass, gbParams(B), gbParams(D)); @@ -1032,7 +1032,7 @@ STELLAR_TYPE EAGB::ResolveEnvelopeLoss(bool p_NoCheck) { } else { m_Radius = R2; - stellarType = STELLAR_TYPE::NAKED_HELIUM_STAR_GIANT_BRANCH; // Has a deep convective envelope + stellarType = STELLAR_TYPE::NAKED_HELIUM_STAR_GIANT_BRANCH; // Has a deep convective envelope } } diff --git a/src/GiantBranch.cpp b/src/GiantBranch.cpp index 05878c085..a3b2565b2 100644 --- a/src/GiantBranch.cpp +++ b/src/GiantBranch.cpp @@ -376,7 +376,7 @@ double GiantBranch::CalculatePerturbationMu() const { * function is called (and does nothing). (So far FGB is the only class that * defines this function where it actually does anything) * - * If DEBUG_PERTURB is defined then perturbation is not disabled while debbuging. + * If DEBUG_PERTURB is defined then perturbation is not disabled while debugging. * To enable perturbation while DEBUG is enabled, define DEBUG_PERTURB. * * @@ -1874,10 +1874,10 @@ STELLAR_TYPE GiantBranch::ResolveElectronCaptureSN() { */ STELLAR_TYPE GiantBranch::ResolveTypeIIaSN() { - m_Mass = 0.0; - m_Radius = 0.0; - m_Luminosity = 0.0; - m_Temperature = 0.0; + m_Mass = 0.0; + m_Radius = 0.0; + m_Luminosity = 0.0; + m_Temperature = 0.0; m_SupernovaDetails.drawnKickMagnitude = 0.0; m_SupernovaDetails.kickMagnitude = 0.0; @@ -1907,6 +1907,7 @@ STELLAR_TYPE GiantBranch::ResolveTypeIIaSN() { */ STELLAR_TYPE GiantBranch::ResolvePairInstabilitySN() { + m_Mass = 0.0; m_Luminosity = 0.0; m_Radius = 0.0; m_Temperature = 0.0; diff --git a/src/HG.cpp b/src/HG.cpp index d1dc9fb70..8240b6b9f 100755 --- a/src/HG.cpp +++ b/src/HG.cpp @@ -1108,7 +1108,7 @@ STELLAR_TYPE HG::ResolveEnvelopeLoss(bool p_NoCheck) { if (p_NoCheck || utils::Compare(m_CoreMass, m_Mass) >= 0) { // envelope loss - m_Mass = std::min(m_CoreMass, m_Mass); + m_Mass = std::min(m_CoreMass, m_Mass); if (utils::Compare(m_Mass0, massCutoffs(MHeF)) < 0) { // star evolves to Helium White Dwarf diff --git a/src/HG.h b/src/HG.h index 676125d63..d17b464dc 100755 --- a/src/HG.h +++ b/src/HG.h @@ -49,11 +49,11 @@ class HG: virtual public BaseStar, public GiantBranch { m_Age = m_Timescales[static_cast(TIMESCALE::tMS)]; CalculateGBParams(); } - m_CoreMass = CalculateCoreMassOnPhase(); - m_COCoreMass = CalculateCOCoreMassOnPhase(); - m_HeCoreMass = CalculateHeCoreMassOnPhase(); - m_Luminosity = CalculateLuminosityOnPhase(); - std::tie(m_Radius, std::ignore) = CalculateRadiusAndStellarTypeOnPhase(); // Update radius + m_CoreMass = CalculateCoreMassOnPhase(); + m_COCoreMass = CalculateCOCoreMassOnPhase(); + m_HeCoreMass = CalculateHeCoreMassOnPhase(); + m_Luminosity = CalculateLuminosityOnPhase(); + std::tie(m_Radius, std::ignore) = CalculateRadiusAndStellarTypeOnPhase(); // Update radius } diff --git a/src/HeHG.cpp b/src/HeHG.cpp index 2ea079985..5c8f60863 100755 --- a/src/HeHG.cpp +++ b/src/HeHG.cpp @@ -448,13 +448,13 @@ STELLAR_TYPE HeHG::ResolveEnvelopeLoss(bool p_NoCheck) { if (p_NoCheck || utils::Compare(m_COCoreMass, m_Mass) >= 0 || m_EnvelopeJustExpelledByPulsations) { // Envelope lost - determine what type of star to form - m_Mass = std::min(m_COCoreMass, m_Mass); - m_CoreMass = m_Mass; - m_HeCoreMass= m_Mass; - m_COCoreMass= m_Mass; - m_Mass0 = m_Mass; - m_Radius = COWD::CalculateRadiusOnPhase_Static(m_Mass); - m_Age = 0.0; + m_Mass = std::min(m_COCoreMass, m_Mass); + m_CoreMass = m_Mass; + m_HeCoreMass = m_Mass; + m_COCoreMass = m_Mass; + m_Mass0 = m_Mass; + m_Radius = COWD::CalculateRadiusOnPhase_Static(m_Mass); + m_Age = 0.0; if (!IsSupernova()) { stellarType = (utils::Compare(m_COCoreMass, OPTIONS->MCBUR1() ) < 0) ? STELLAR_TYPE::CARBON_OXYGEN_WHITE_DWARF : STELLAR_TYPE::OXYGEN_NEON_WHITE_DWARF; //Note that this uses the CO core mass, rather than the core mass at base of AGB or He mass at He star birth suggested by Hurley+, 2000 } diff --git a/src/HeMS.cpp b/src/HeMS.cpp index 8febed40a..14535b337 100755 --- a/src/HeMS.cpp +++ b/src/HeMS.cpp @@ -464,8 +464,8 @@ STELLAR_TYPE HeMS::ResolveEnvelopeLoss(bool p_NoCheck) { if (p_NoCheck || utils::Compare(m_Mass, 0.0) <= 0) { stellarType = STELLAR_TYPE::MASSLESS_REMNANT; - m_Radius = 0.0; // massless remnant - m_Mass = 0.0; + m_Radius = 0.0; // massless remnant + m_Mass = 0.0; } return stellarType; diff --git a/src/Log.h b/src/Log.h index 01e2d75cf..2610dec71 100755 --- a/src/Log.h +++ b/src/Log.h @@ -378,9 +378,9 @@ class FormatVariantValue: public boost::static_visitor { string operator()(const unsigned short int v, const string fmtStr) const { string fmt = fmtStr; fmt = "%" + fmt + "u"; return utils::vFormat(fmt.c_str(), v); } string operator()(const unsigned long int v, const string fmtStr) const { string fmt = fmtStr; fmt = "%" + fmt + "u"; return utils::vFormat(fmt.c_str(), v); } // also handles OBJECT_ID (typedef) string operator()(const unsigned long long int v, const string fmtStr) const { string fmt = fmtStr; fmt = "%" + fmt + "u"; return utils::vFormat(fmt.c_str(), v); } - string operator()(const float v, const string fmtStr) const { string fmt = fmtStr; fmt = "%" + fmt + "e"; return utils::vFormat(fmt.c_str(), v); } - string operator()(const double v, const string fmtStr) const { string fmt = fmtStr; fmt = "%" + fmt + "e"; return utils::vFormat(fmt.c_str(), v); } - string operator()(const long double v, const string fmtStr) const { string fmt = fmtStr; fmt = "%" + fmt + "e"; return utils::vFormat(fmt.c_str(), v); } + string operator()(const float v, const string fmtStr) const { string fmt = fmtStr; fmt = "%" + fmt + "f"; return utils::vFormat(fmt.c_str(), v); } + string operator()(const double v, const string fmtStr) const { string fmt = fmtStr; fmt = "%" + fmt + "f"; return utils::vFormat(fmt.c_str(), v); } + string operator()(const long double v, const string fmtStr) const { string fmt = fmtStr; fmt = "%" + fmt + "f"; return utils::vFormat(fmt.c_str(), v); } string operator()(const string v, const string fmtStr) const { string fmt = fmtStr; fmt = "%-" + fmt + "s"; return utils::vFormat(fmt.c_str(), v.c_str()); } string operator()(const ERROR v, const string fmtStr) const { string fmt = fmtStr; fmt = "%" + fmt + "d"; return utils::vFormat(fmt.c_str(), static_cast(v)); } string operator()(const STELLAR_TYPE v, const string fmtStr) const { string fmt = fmtStr; fmt = "%" + fmt + "d"; return utils::vFormat(fmt.c_str(), static_cast(v)); } diff --git a/src/MainSequence.cpp b/src/MainSequence.cpp index d5b9f2f22..58b4b3659 100644 --- a/src/MainSequence.cpp +++ b/src/MainSequence.cpp @@ -637,13 +637,14 @@ STELLAR_TYPE MainSequence::ResolveEnvelopeLoss(bool p_NoCheck) { if (p_NoCheck || utils::Compare(m_Mass, 0.0) <= 0) { stellarType = STELLAR_TYPE::MASSLESS_REMNANT; - m_Radius = 0.0; // massless remnant - m_Mass = 0.0; + m_Radius = 0.0; // massless remnant + m_Mass = 0.0; } return stellarType; } + /* * Update the minimum core mass of a main sequence star that loses mass through Case A mass transfer by * setting it equal to the core mass of a TAMS star, scaled by the fractional age. @@ -658,9 +659,9 @@ STELLAR_TYPE MainSequence::ResolveEnvelopeLoss(bool p_NoCheck) { void MainSequence::UpdateMinimumCoreMass() { if (OPTIONS->RetainCoreMassDuringCaseAMassTransfer()) { - double fractionalAge=CalculateTauOnPhase(); - HG clone = *this; //create an HG star clone to query its core mass just after TAMS - double TAMSCoreMass = clone.CoreMass(); - m_MinimumCoreMass = std::max(m_MinimumCoreMass, fractionalAge * TAMSCoreMass); + double fractionalAge =CalculateTauOnPhase(); + HG clone = *this; //create an HG star clone to query its core mass just after TAMS + double TAMSCoreMass = clone.CoreMass(); + m_MinimumCoreMass = std::max(m_MinimumCoreMass, fractionalAge * TAMSCoreMass); } } diff --git a/src/Options.cpp b/src/Options.cpp index 2b3267487..c299fd8d3 100644 --- a/src/Options.cpp +++ b/src/Options.cpp @@ -194,6 +194,7 @@ void Options::OptionValues::Initialise() { m_MaxEvolutionTime = 13700.0; m_MaxNumberOfTimestepIterations = 99999; m_TimestepMultiplier = 1.0; + m_TimestepsFileName = ""; // Initial mass options m_InitialMass = 5.0; @@ -365,8 +366,8 @@ void Options::OptionValues::Initialise() { m_LuminousBlueVariablePrescription.type = LBV_PRESCRIPTION::HURLEY_ADD; m_LuminousBlueVariablePrescription.typeString = LBV_PRESCRIPTION_LABEL.at(m_LuminousBlueVariablePrescription.type); - m_OBMassLoss.type = OB_MASS_LOSS::VINK2021; - m_OBMassLoss.typeString = OB_MASS_LOSS_LABEL.at(m_OBMassLoss.type); + m_OBMassLoss.type = OB_MASS_LOSS::VINK2021; + m_OBMassLoss.typeString = OB_MASS_LOSS_LABEL.at(m_OBMassLoss.type); m_VMSMassLoss.type = VMS_MASS_LOSS::SABHAHIT2023; m_VMSMassLoss.typeString = VMS_MASS_LOSS_LABEL.at(m_VMSMassLoss.type); @@ -374,8 +375,8 @@ void Options::OptionValues::Initialise() { m_RSGMassLoss.type = RSG_MASS_LOSS::DECIN2023; m_RSGMassLoss.typeString = RSG_MASS_LOSS_LABEL.at(m_RSGMassLoss.type); - m_WRMassLoss.type = WR_MASS_LOSS::SANDERVINK2023; - m_WRMassLoss.typeString = WR_MASS_LOSS_LABEL.at(m_WRMassLoss.type); + m_WRMassLoss.type = WR_MASS_LOSS::SANDERVINK2023; + m_WRMassLoss.typeString = WR_MASS_LOSS_LABEL.at(m_WRMassLoss.type); // Wind mass loss multiplicitive constants m_CoolWindMassLossMultiplier = 1.0; @@ -446,6 +447,7 @@ void Options::OptionValues::Initialise() { m_MassTransferCriticalMassRatioWhiteDwarfNonDegenerateAccretor = 0.0; // Claeys+ 2014 = unspecified m_MassTransferCriticalMassRatioWhiteDwarfDegenerateAccretor = 1.6; // Claeys+ 2014 = 1.6 + // Common Envelope options m_CommonEnvelopeAlpha = 1.0; m_CommonEnvelopeLambda = 0.1; @@ -969,7 +971,7 @@ bool Options::AddOptions(OptionValues *p_Options, po::options_description *p_Opt ( "maximum-number-timestep-iterations", - po::value(&p_Options->m_MaxNumberOfTimestepIterations)->default_value(p_Options->m_MaxNumberOfTimestepIterations), + po::value(&p_Options->m_MaxNumberOfTimestepIterations)->default_value(p_Options->m_MaxNumberOfTimestepIterations), ("Maximum number of timesteps to evolve binary before giving up (default = " + std::to_string(p_Options->m_MaxNumberOfTimestepIterations) + ")").c_str() ) @@ -1804,6 +1806,12 @@ bool Options::AddOptions(OptionValues *p_Options, po::options_description *p_Opt ("Prescription for stellar zeta (" + AllowedOptionValuesFormatted("stellar-zeta-prescription") + ", default = '" + p_Options->m_StellarZetaPrescription.typeString + "')").c_str() ) + ( + "timesteps-filename", + po::value(&p_Options->m_TimestepsFileName)->default_value(p_Options->m_TimestepsFileName), + ("Filename for file to provide timesteps to be used for evolution (SSE and BSE) (default = " + p_Options->m_TimestepsFileName + ")").c_str() + ) + ( "YAML-template", po::value(&p_Options->m_YAMLtemplate)->default_value(p_Options->m_YAMLtemplate), diff --git a/src/Options.h b/src/Options.h index e17b1f034..ec4c7c3e1 100755 --- a/src/Options.h +++ b/src/Options.h @@ -529,6 +529,8 @@ class Options { "store-input-files", "switch-log", + "timesteps-filename", + "use-mass-loss", "VMS-mass-loss", @@ -608,6 +610,8 @@ class Options { "store-input-files", "switch-log", + "timesteps-filename", + "version", "v", "yaml-template" @@ -679,12 +683,14 @@ class Options { unsigned long int m_RandomSeed; // Random seed to use double m_MaxEvolutionTime; // Maximum time to evolve a binary by - int m_MaxNumberOfTimestepIterations; // Maximum number of timesteps to evolve binary for before giving up + long unsigned int m_MaxNumberOfTimestepIterations; // Maximum number of timesteps to evolve binary for before giving up double m_TimestepMultiplier; // Multiplier for time step size (<1 -- shorter timesteps, >1 -- longer timesteps) std::streamsize m_GridStartLine; // The grid file line to start processing (0-based) std::streamsize m_GridLinesToProcess; // The number of grid file lines to process (starting at m_GridStartLine) + std::string m_TimestepsFileName; // The name of the timesteps file + // Initial distribution variables double m_InitialMass; // Initial mass of single star (SSE) @@ -1369,7 +1375,7 @@ class Options { MT_THERMALLY_LIMITED_VARIATION MassTransferThermallyLimitedVariation() const { return OPT_VALUE("mass-transfer-thermal-limit-accretor", m_MassTransferThermallyLimitedVariation.type, true); } double MaxEvolutionTime() const { return OPT_VALUE("maximum-evolution-time", m_MaxEvolutionTime, true); } double MaximumNeutronStarMass() const { return OPT_VALUE("maximum-neutron-star-mass", m_MaximumNeutronStarMass, true); } - int MaxNumberOfTimestepIterations() const { return OPT_VALUE("maximum-number-timestep-iterations", m_MaxNumberOfTimestepIterations, true); } + long unsigned int MaxNumberOfTimestepIterations() const { return OPT_VALUE("maximum-number-timestep-iterations", m_MaxNumberOfTimestepIterations, true); } double MaximumDonorMass() const { return OPT_VALUE("maximum-mass-donor-nandez-ivanova", m_MaximumMassDonorNandezIvanova, true); } double MCBUR1() const { return OPT_VALUE("mcbur1", m_mCBUR1, true); } @@ -1474,6 +1480,7 @@ class Options { ZETA_PRESCRIPTION StellarZetaPrescription() const { return OPT_VALUE("stellar-zeta-prescription", m_StellarZetaPrescription.type, true); } + std::string TimestepsFileName() const { return OPT_VALUE("timesteps-filename", m_TimestepsFileName, true); } double TimestepMultiplier() const { return m_CmdLine.optionValues.m_TimestepMultiplier; } bool UseFixedUK() const { return (m_GridLine.optionValues.m_UseFixedUK || m_CmdLine.optionValues.m_UseFixedUK); } diff --git a/src/Remnants.h b/src/Remnants.h index 24cdb94aa..284315988 100644 --- a/src/Remnants.h +++ b/src/Remnants.h @@ -95,7 +95,7 @@ class Remnants: virtual public BaseStar, public HeGB { void ResolveEnvelopeMassAtPhaseEnd(const double p_Tau) const { ResolveEnvelopeMassOnPhase(p_Tau); } // Same as on phase void ResolveEnvelopeMassOnPhase(const double p_Tau) const { } // NO-OP - void ResolveMassLoss() { } // NO-OP + void ResolveMassLoss(const bool p_UpdateMDt = true) { } // NO-OP STELLAR_TYPE ResolveSkippedPhase() { return BaseStar::ResolveSkippedPhase(); } // Default to BaseStar // diff --git a/src/Star.cpp b/src/Star.cpp index c63bc9cb1..2801db625 100644 --- a/src/Star.cpp +++ b/src/Star.cpp @@ -330,8 +330,8 @@ STELLAR_TYPE Star::UpdateAttributesAndAgeOneTimestep(const double p_DeltaMass, * * * Checks whether the star: - * - is a massless remnant (checked after applying p_DeltaMass and p_DeltaMass0) - * - has become a supernova (checked after applying p_DeltaMass and p_DeltaMass0) + * - is a massless remnant (checked after applying p_DeltaMass and p_DeltaMass0) <--- this is no longer ture (it may have been once) + * - has become a supernova (checked after applying p_DeltaMass and p_DeltaMass0) <--- this is no longer ture (it may have been once) JR: this causes different behaviour SSE vs BSE * - should skip this phase for this timestep (checked after applying p_DeltaMass and p_DeltaMass0) * * If none of the above are true the star's attributes are updated based on the mass changes required and the star's @@ -397,12 +397,19 @@ STELLAR_TYPE Star::AgeOneTimestep(const double p_DeltaTime, bool p_Switch) { * * The functional return is the timestep actually taken (in Myr) * - * double EvolveOneTimestep(const double p_Dt) + * double EvolveOneTimestep(const double p_Dt, const bool p_Force) * * @param [IN] p_Dt The suggested timestep to evolve + * @param [IN] p_Force Flag to force using the timestep passed (default = FALSE) + * If p_Force is TRUE + * - radial change is not checked and the timestep will not be shortened + * - the class member m_Dt will not be updated in BaseStar::CalculateMassLossValues() + * If p_Force is FALSE + * - radial change is checked and the timestep may be shortened + * - the class member m_Dt may be updated in BaseStar::CalculateMassLossValues() * @return The timestep actually taken */ -double Star::EvolveOneTimestep(const double p_Dt) { +double Star::EvolveOneTimestep(const double p_Dt, const bool p_Force) { double dt = p_Dt; @@ -415,7 +422,7 @@ double Star::EvolveOneTimestep(const double p_Dt) { SaveState(); // save the state of the star - in case we want to revert - double minTimestep = std::max(m_Star->CalculateDynamicalTimescale(), ABSOLUTE_MINIMUM_TIMESTEP); // calculate the minimum timestep - maximum of dynamical timescale for this star and the absolute minimum timestep + double minTimestep = std::max(m_Star->CalculateDynamicalTimescale(), NUCLEAR_MINIMUM_TIMESTEP); // calculate the minimum timestep - maximum of dynamical timescale for this star and the absolute minimum timestep // evolve the star a single timestep @@ -425,7 +432,7 @@ double Star::EvolveOneTimestep(const double p_Dt) { // don't take the timestep if we stepped too far takeTimestep = true; // flag to determine if the timestep should be taken - if (utils::Compare(m_Star->CalculateRadialChange(), MAXIMUM_RADIAL_CHANGE) >= 0) { // too much change? + if (!p_Force && utils::Compare(m_Star->CalculateRadialChange(), MAXIMUM_RADIAL_CHANGE) >= 0) { // too much change? if (utils::Compare(dt, minTimestep) <= 0) { // yes - already at or below minimum timestep? takeTimestep = true; // yes - just take the last timestep SHOW_WARN(ERROR::TIMESTEP_BELOW_MINIMUM); // announce the problem if required and plough on regardless... @@ -452,11 +459,15 @@ double Star::EvolveOneTimestep(const double p_Dt) { // take the timestep - (void)m_Star->PrintStashedSupernovaDetails(); // print stashed SSE Supernova log record if necessary - (void)SwitchTo(stellarType); // switch phase if required JR: whether this goes before or after the log record is a little problematic, but in the end probably doesn't matter too much - (void)m_Star->ResolveMassLoss(); // apply wind mass loss if required JR: should this really be before the call to SwitchTo()? It isn't in the original code + (void)m_Star->PrintDetailedOutput(m_Id, SSE_DETAILED_RECORD_TYPE::PRE_MASS_LOSS); // log record - pre mass loss + + (void)m_Star->ResolveMassLoss(!p_Force); // apply wind mass loss if required - after the switch is consistent with the BSE code + + (void)m_Star->PrintStashedSupernovaDetails(); // print stashed SSE Supernova log record if necessary + + (void)m_Star->PrintDetailedOutput(m_Id, SSE_DETAILED_RECORD_TYPE::POST_MASS_LOSS); // log record - post mass loss return dt; // return the timestep actually taken } @@ -475,44 +486,81 @@ EVOLUTION_STATUS Star::Evolve(const long int p_Id) { EVOLUTION_STATUS evolutionStatus = EVOLUTION_STATUS::CONTINUE; - m_Id = p_Id; // store the id + m_Id = p_Id; // store the id // evolve the star - m_Star->CalculateGBParams(); // calculate giant branch parameters - in case for some reason star is initially not MS + m_Star->CalculateGBParams(); // calculate giant branch parameters - in case for some reason star is initially not MS double dt = 0.0; + (void)m_Star->PrintDetailedOutput(m_Id, SSE_DETAILED_RECORD_TYPE::INITIAL_STATE); // log record + + bool usingProvidedTimesteps = false; // using user-provided timesteps? + DBL_VECTOR timesteps; + if (!OPTIONS->TimestepsFileName().empty()) { // have timesteps filename? + // yes + ERROR error; + std::tie(error, timesteps) = utils::ReadTimesteps(OPTIONS->TimestepsFileName()); // read timesteps from file + if (error != ERROR::NONE) { // ok? + evolutionStatus = EVOLUTION_STATUS::NO_TIMESTEPS; // no - set status + SHOW_WARN(error); // show warning + } + else usingProvidedTimesteps = true; // have user-provided timesteps + } + // JR: todo: // m_Error seems to be set ad hoc for SSE, and doesn't actually stop the evolution // we should be more rigorous in checking/setting error conditions, and stop the evolution for catastrophic errors - int stepNum = 0; + long unsigned int stepNum = 0; while (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { - if (m_Star->Time() > OPTIONS->MaxEvolutionTime()) { - evolutionStatus = EVOLUTION_STATUS::TIMES_UP; // out of time... + evolutionStatus = EVOLUTION_STATUS::TIMES_UP; // out of time... } else if (stepNum >= OPTIONS->MaxNumberOfTimestepIterations()) { - evolutionStatus = EVOLUTION_STATUS::STEPS_UP; // out of steps... + evolutionStatus = EVOLUTION_STATUS::STEPS_UP; // out of steps... } else if (!m_Star->IsOneOf({ STELLAR_TYPE::MS_LTE_07, STELLAR_TYPE::MS_GT_07, STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS, STELLAR_TYPE::HERTZSPRUNG_GAP, STELLAR_TYPE::FIRST_GIANT_BRANCH, STELLAR_TYPE::CORE_HELIUM_BURNING, STELLAR_TYPE::EARLY_ASYMPTOTIC_GIANT_BRANCH, STELLAR_TYPE::THERMALLY_PULSING_ASYMPTOTIC_GIANT_BRANCH, STELLAR_TYPE::NAKED_HELIUM_STAR_MS, STELLAR_TYPE::NAKED_HELIUM_STAR_HERTZSPRUNG_GAP, STELLAR_TYPE::NAKED_HELIUM_STAR_GIANT_BRANCH })) { - evolutionStatus = EVOLUTION_STATUS::DONE; // we're done + evolutionStatus = EVOLUTION_STATUS::DONE; // we're done + } + else if (usingProvidedTimesteps && stepNum >= timesteps.size()) { + evolutionStatus = EVOLUTION_STATUS::TIMESTEPS_EXHAUSTED; // using user-provided timesteps and all consumed + SHOW_WARN(ERROR::TIMESTEPS_EXHAUSTED); // show warning } - else { - stepNum++; // increment step number - dt = m_Star->CalculateTimestep() * OPTIONS->TimestepMultiplier(); // calculate new timestep - EvolveOneTimestep(dt); // evolve for timestep - (void)m_Star->PrintDetailedOutput(m_Id); // log record + else { // evolve one timestep + if (usingProvidedTimesteps) { // user-provided timesteps + // get new timestep + // - don't quantise + // - don't apply timestep multiplier + // (we assume user wants the timesteps in the file) + dt = timesteps[stepNum]; + } + else { // not using user-provided timesteps + dt = m_Star->CalculateTimestep() * OPTIONS->TimestepMultiplier(); // calculate new timestep + dt = std::round(dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM; // quantised + } + stepNum++; // increment step number + + EvolveOneTimestep(dt, true); // evolve for timestep + UpdateAttributes(0.0, 0.0, true); // JR: if this is not included, BSE and SSE are out of sync by 1 timestep. If we remove this, we have to change BS accordingly. Not sure which one is right yet... (or if that actually matters) + (void)m_Star->PrintDetailedOutput(m_Id, SSE_DETAILED_RECORD_TYPE::TIMESTEP_COMPLETED); // log record } } - m_Star->SetEvolutionStatus(evolutionStatus); // set evolution final outcome for star + if (usingProvidedTimesteps && timesteps.size() > stepNum) { // all user-defined timesteps consumed? + evolutionStatus = EVOLUTION_STATUS::TIMESTEPS_NOT_CONSUMED; // no - set status + SHOW_WARN(ERROR::TIMESTEPS_NOT_CONSUMED); // show warning + } + + (void)m_Star->PrintDetailedOutput(m_Id, SSE_DETAILED_RECORD_TYPE::FINAL_STATE); // log record + + m_Star->SetEvolutionStatus(evolutionStatus); // set evolution final outcome for star - (void)m_Star->PrintSystemParameters(); // log system parameters + (void)m_Star->PrintSystemParameters(); // log system parameters return evolutionStatus; } diff --git a/src/Star.h b/src/Star.h index 674aa2c91..fc1c79416 100755 --- a/src/Star.h +++ b/src/Star.h @@ -204,7 +204,7 @@ class Star { EVOLUTION_STATUS Evolve(const long int p_Id); - double EvolveOneTimestep(const double p_Dt); + double EvolveOneTimestep(const double p_Dt, const bool p_Force = false); double InterpolateGe20QCrit(const QCRIT_PRESCRIPTION p_qCritPrescription) { return m_Star->InterpolateGe20QCrit(p_qCritPrescription); } void HaltWinds() { m_Star->HaltWinds(); } @@ -268,6 +268,8 @@ class Star { BaseStar *m_Star; // pointer to current star BaseStar *m_SaveStar; // pointer to saved star + std::vector m_Timesteps; // timesteps vector - for debugging/testing + }; #endif // __Star_h__ diff --git a/src/TPAGB.cpp b/src/TPAGB.cpp index 6e2f435a5..60da3aae0 100755 --- a/src/TPAGB.cpp +++ b/src/TPAGB.cpp @@ -898,14 +898,14 @@ STELLAR_TYPE TPAGB::ResolveEnvelopeLoss(bool p_NoCheck) { if (p_NoCheck || (utils::Compare(m_CoreMass, m_Mass)) >= 0 || m_EnvelopeJustExpelledByPulsations) { - m_Mass = std::min(m_CoreMass, m_Mass); - m_CoreMass = m_Mass; - m_HeCoreMass= m_Mass; - m_COCoreMass= m_Mass; - m_Mass0 = m_Mass; - m_Radius = COWD::CalculateRadiusOnPhase_Static(m_Mass); - m_Age = 0.0; - stellarType = (utils::Compare(gbParams(McBAGB), OPTIONS->MCBUR1() ) < 0) ? STELLAR_TYPE::CARBON_OXYGEN_WHITE_DWARF : STELLAR_TYPE::OXYGEN_NEON_WHITE_DWARF; + m_Mass = std::min(m_CoreMass, m_Mass); + m_CoreMass = m_Mass; + m_HeCoreMass = m_Mass; + m_COCoreMass = m_Mass; + m_Mass0 = m_Mass; + m_Radius = COWD::CalculateRadiusOnPhase_Static(m_Mass); + m_Age = 0.0; + stellarType = (utils::Compare(gbParams(McBAGB), OPTIONS->MCBUR1() ) < 0) ? STELLAR_TYPE::CARBON_OXYGEN_WHITE_DWARF : STELLAR_TYPE::OXYGEN_NEON_WHITE_DWARF; } return stellarType; diff --git a/src/changelog.h b/src/changelog.h index 04296b406..1a93b6ea2 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1088,8 +1088,9 @@ // - Change to functionality (noted above) noted in 'What's New' online documentation page // 02.41.04 JR - Dec 30, 2023 - Defect repair: // - Fix for issue #1048 +// 02.42.99 JR - Jan 04. 2024 - Enhancements etc. TBC -const std::string VERSION_STRING = "02.41.04"; +const std::string VERSION_STRING = "02.42.00"; # endif // __changelog_h__ diff --git a/src/constants.h b/src/constants.h index 064776f03..723959e6c 100755 --- a/src/constants.h +++ b/src/constants.h @@ -215,12 +215,14 @@ constexpr double KM_TO_RSOL = 1.0 / RSOL_TO_KM; constexpr double KM_TO_AU = 1.0 / AU_TO_KM; // convert km to Astronomical Units AU // time -constexpr double SECONDS_IN_YEAR = 31556926.0; // number of second in 1 year constexpr double DAYS_IN_QUAD = 1461.0; // number of days in any 4-year period -constexpr double DAYS_IN_YEAR = DAYS_IN_QUAD / 4.0; // mean days per year (given DAYS_IN_QUAD) -constexpr double SECONDS_IN_DAY = SECONDS_IN_YEAR / DAYS_IN_YEAR; // number of second in 1 day -constexpr double SECONDS_IN_MS = 1.0E-3; // number of second in 1 millisecond -constexpr double SECONDS_IN_MYR = 31556926.0 * 1.0E6; // number of second in 1 Myr +constexpr double DAYS_IN_YEAR = DAYS_IN_QUAD / 4.0; // mean days per year, given DAYS_IN_QUAD +constexpr double SECONDS_IN_DAY = 24.0 * 60.0 * 60.0; // number of seconds in 1 day +constexpr double SECONDS_IN_QUAD = DAYS_IN_QUAD * SECONDS_IN_DAY; // number of seconds in 1 quad +constexpr double SECONDS_IN_YEAR = SECONDS_IN_QUAD / 4.0; // number of seconds in 1 year +constexpr double SECONDS_IN_MYR = SECONDS_IN_YEAR * 1.0E6; // number of seconds in 1 Myr +constexpr double SECONDS_IN_MS = 1.0E-3; // number of seconds in 1 millisecond + constexpr double MYR_TO_YEAR = 1.0E6; // convert Myr to year constexpr double YEAR_TO_MYR = 1.0E-6; // convert year to Myr @@ -259,7 +261,7 @@ constexpr double G_SOLAR_YEAR = 3.14E7; constexpr double RSOL = 6.957E8; // Solar Radius (in m) constexpr double ZSOL = 0.02; // Solar Metallicity used in scalings constexpr double LOG10_ZSOL = -1.69897; // log10(ZSOL) - for performance -constexpr double ZSOL_ASPLUND = 0.0142; // Solar Metallicity (Asplund+ 2010) used in initial condition +constexpr double ZSOL_ASPLUND = 0.0142; // Solar Metallicity (Asplund+ 2010) used in initial condition constexpr double TSOL = 5778.0; // Solar Temperature in kelvin constexpr double LSOL = 3.844E33; // Solar Luminosity in erg/s constexpr double LSOLW = 3.844E26; // Solar luminosity (in W) @@ -301,8 +303,11 @@ constexpr double LBV_LUMINOSITY_LIMIT_VANBEVEREN = 3.0E5; constexpr double CONVECTIVE_BOUNDARY_TEMPERATURE_BELCZYNSKI = 5.37E3; // Threshold temperature for the star to develop a convective envelope, in Kelvin (10^3.73 K, from Belczynski+, 2008) -constexpr double ABSOLUTE_MINIMUM_TIMESTEP = 100.0 / SECONDS_IN_MYR; // 100 seconds expressed in Myr (3.1688765E-12 Myr) -constexpr double NUCLEAR_MINIMUM_TIMESTEP = 1.0E-6; // Minimum time step for nuclear evolution = 1 year expressed in Myr +constexpr double TIMESTEP_QUANTUM = 1.0E-12; // Timestep quantum in Myr (=31.5576 seconds, given DAYS_IN_QUAD) +constexpr double ABSOLUTE_MINIMUM_TIMESTEP = 3.0 * TIMESTEP_QUANTUM; // In Myr (=94.6728 seconds, given TIMESTEP QUANTUM) +constexpr double NUCLEAR_MINIMUM_TIMESTEP = 1.0E6 * TIMESTEP_QUANTUM; // Minimum time step for nuclear evolution in My (= 1 year = 31577600 seconds, given TIMESTEP_QUANTUM) + +constexpr unsigned int ABSOLUTE_MAXIMUM_TIMESTEPS = 1000000; // Absolute maximum number of timesteps constexpr int MAX_BSE_INITIAL_CONDITIONS_ITERATIONS = 100; // Maximum loop iterations looking for initial conditions for binary systems constexpr int MAX_TIMESTEP_RETRIES = 30; // Maximum retries to find a good timestep for stellar evolution @@ -395,11 +400,15 @@ enum class BSE_DETAILED_RECORD_TYPE: unsigned int { POST_CEE, // 12 - record was logged immediately following a common envelope event POST_SN, // 13 - record was logged immediately following a supernova event POST_MASS_RESOLUTION, // 14 - record was logged immediately following mass resolution (i.e. after winds mass loss & mass transfer complete) - POST_MASS_RESOLUTION_MERGER // 15 - record was logged immediately following a merger after mass resolution + POST_MASS_RESOLUTION_MERGER }; enum class SSE_DETAILED_RECORD_TYPE: unsigned int { // SSE_DETAILED_OUTPUT file record type - DEFAULT = 1 // 1 - default SSE_DETAILED_OUTPUT record type + INITIAL_STATE = 1, // 1 - record describes the initial state of the star + PRE_MASS_LOSS, // 2 - record was logged after timestep taken, but before mass loss resolution + POST_MASS_LOSS, // 3 - record was logged after after mass loss resolution + TIMESTEP_COMPLETED, // 4 - record was logged immediately following the completion of the timestep (after all changes to the star) + FINAL_STATE // 5 - record describes the final state of the star }; enum class BSE_SN_RECORD_TYPE: unsigned int { // BSE_SUPERNOVAE file record type @@ -556,6 +565,8 @@ enum class ERROR: int { BOOST_OPTION_CMDLINE, // failed to initialise Boost options descriptions for commandline options BOOST_OPTION_GRIDLINE, // failed to initialise Boost options descriptions for grid line options BOOST_OPTION_INTERNAL_ERROR, // Boost option internal error + EMPTY_FILE, // filename is empty (contains no content) + EMPTY_FILENAME, // filename is an empty string ERROR, // unspecified error ERROR_PROCESSING_CMDLINE_OPTIONS, // an error occurred while processing commandline options ERROR_PROCESSING_GRIDLINE_OPTIONS, // an error occurred while processing grid file options @@ -570,7 +581,6 @@ enum class ERROR: int { EXPECTED_PROPERTY_SPECIFIER, // expected a valid property specifier EXPECTED_SN_EVENT, // expected a supernova event EXPECTED_STELLAR_PROPERTY, // expected a stellar property (STAR_PROPERTY) - EMPTY_FILENAME, // filename is an empty string FILE_DOES_NOT_EXIST, // file does not exist FILE_NOT_CLOSED, // error closing file - file not closed FILE_OPEN_ERROR, // error opening file @@ -587,7 +597,8 @@ enum class ERROR: int { INVALID_TYPE_MT_MASS_RATIO, // invalid stellar type for mass ratio calculation INVALID_TYPE_MT_THERMAL_TIMESCALE, // invalid stellar type for thermal timescale calculation INVALID_TYPE_ZETA_CALCULATION, // invalid stellar type for Zeta calculation - INVALID_VALUE_FOR_BOOLEAN_OPTION, // invalid values specified for boolean option + INVALID_VALUE_FOR_BOOLEAN_OPTION, // invalid value specified for boolean option + INVALID_VALUE_IN_TIMESTEPS_FILE, // invalid value in timesteps file LAMBDA_NOT_POSITIVE, // lambda is <= 0.0 - invalid LOW_GAMMA, // very massive mass-loss prescription being extrapolated to low gamma (<0.5) LOW_TEFF_WINDS, // winds being used at low temperature @@ -611,9 +622,12 @@ enum class ERROR: int { STELLAR_SIMULATION_STOPPED, // stellar simulation stopped SUGGEST_HELP, // suggest using --help TIMESTEP_BELOW_MINIMUM, // timestep too small - below minimum + TIMESTEPS_EXHAUSTED, // timesteps provided exhausted, but evolution not complete + TIMESTEPS_NOT_CONSUMED, // evolution complete, but provided timesteps not consumed TOO_MANY_MASS0_ITERATIONS, // too many iterations in MASS0 root finder TOO_MANY_OMEGA_ITERATIONS, // too many iterations in OMEGA root finder TOO_MANY_RLOF_ITERATIONS, // too many iterations in RLOF root finder + TOO_MANY_TIMESTEPS_IN_TIMESTEPS_FILE, // too many timesteps in timesteps file (exceeds maximum) UNEXPECTED_END_OF_FILE, // unexpected end of file UNEXPECTED_LOG_FILE_TYPE, // unexpected log file type UNEXPECTED_SN_EVENT, // unexpected supernova event in this context @@ -698,6 +712,8 @@ const COMPASUnorderedMap> ERROR_CATA { ERROR::BOOST_OPTION_CMDLINE, { ERROR_SCOPE::ALWAYS, "Failed to initialise Boost options descriptions for commandline options" }}, { ERROR::BOOST_OPTION_GRIDLINE, { ERROR_SCOPE::ALWAYS, "Failed to initialise Boost options descriptions for grid line options" }}, { ERROR::BOOST_OPTION_INTERNAL_ERROR, { ERROR_SCOPE::ALWAYS, "Internal error: Boost vm, option" }}, + { ERROR::EMPTY_FILE, { ERROR_SCOPE::ALWAYS, "File is empty" }}, + { ERROR::EMPTY_FILENAME, { ERROR_SCOPE::ALWAYS, "Filename is an empty string" }}, { ERROR::ERROR, { ERROR_SCOPE::ALWAYS, "Error!" }}, { ERROR::ERROR_PROCESSING_CMDLINE_OPTIONS, { ERROR_SCOPE::ALWAYS, "An error occurred while processing commandline options" }}, { ERROR::ERROR_PROCESSING_GRIDLINE_OPTIONS, { ERROR_SCOPE::ALWAYS, "An error occurred while processing grid file options" }}, @@ -712,7 +728,6 @@ const COMPASUnorderedMap> ERROR_CATA { ERROR::EXPECTED_PROPERTY_SPECIFIER, { ERROR_SCOPE::ALWAYS, "Expected a property specifier or close brace '}'" }}, { ERROR::EXPECTED_SN_EVENT, { ERROR_SCOPE::ALWAYS, "Expected a supernova event" }}, { ERROR::EXPECTED_STELLAR_PROPERTY, { ERROR_SCOPE::ALWAYS, "Expected stellar logfile property: one of { STAR_PROPERTY, PROGRAM_OPTION }" }}, - { ERROR::EMPTY_FILENAME, { ERROR_SCOPE::ALWAYS, "Filename is an empty string" }}, { ERROR::FILE_DOES_NOT_EXIST, { ERROR_SCOPE::ALWAYS, "File does not exist" }}, { ERROR::FILE_NOT_CLOSED, { ERROR_SCOPE::ALWAYS, "Error closing file - file not closed" }}, { ERROR::FILE_OPEN_ERROR, { ERROR_SCOPE::ALWAYS, "Error opening file" }}, @@ -730,6 +745,7 @@ const COMPASUnorderedMap> ERROR_CATA { ERROR::INVALID_TYPE_MT_THERMAL_TIMESCALE, { ERROR_SCOPE::ALWAYS, "Invalid stellar type for thermal timescale calculation" }}, { ERROR::INVALID_TYPE_ZETA_CALCULATION, { ERROR_SCOPE::ALWAYS, "Invalid stellar type for Zeta calculation" }}, { ERROR::INVALID_VALUE_FOR_BOOLEAN_OPTION, { ERROR_SCOPE::ALWAYS, "Invalid value specified for BOOLEAN option" }}, + { ERROR::INVALID_VALUE_IN_TIMESTEPS_FILE, { ERROR_SCOPE::ALWAYS, "Invalid value in timesteps file" }}, { ERROR::LAMBDA_NOT_POSITIVE, { ERROR_SCOPE::ALWAYS, "Lambda <= 0.0" }}, { ERROR::LOW_GAMMA, { ERROR_SCOPE::ALWAYS, "Very massive prescription being extrapolated to low gamma (<0.5)" }}, { ERROR::LOW_TEFF_WINDS, { ERROR_SCOPE::ALWAYS, "Winds being used at low temperature" }}, @@ -754,9 +770,12 @@ const COMPASUnorderedMap> ERROR_CATA { ERROR::STELLAR_SIMULATION_STOPPED, { ERROR_SCOPE::ALWAYS, "Stellar simulation stopped" }}, { ERROR::SUGGEST_HELP, { ERROR_SCOPE::ALWAYS, "Use option '-h' (or '--help') to see (descriptions of) available options" }}, { ERROR::TIMESTEP_BELOW_MINIMUM, { ERROR_SCOPE::ALWAYS, "Timestep below minimum - timestep taken" }}, + { ERROR::TIMESTEPS_EXHAUSTED, { ERROR_SCOPE::ALWAYS, "Provided timesteps exhausted, but evolution not complete" }}, + { ERROR::TIMESTEPS_NOT_CONSUMED, { ERROR_SCOPE::ALWAYS, "Evolution complete, but provided timesteps not consumed" }}, { ERROR::TOO_MANY_MASS0_ITERATIONS, { ERROR_SCOPE::ALWAYS, "Reached maximum number of iterations when looking for effective initial mass Mass_0 to match desired stellar core of HG star following case A mass transfer" }}, { ERROR::TOO_MANY_OMEGA_ITERATIONS, { ERROR_SCOPE::ALWAYS, "Reached maximum number of iterations when looking for omega when circularising and synchronising for tides" }}, { ERROR::TOO_MANY_RLOF_ITERATIONS, { ERROR_SCOPE::ALWAYS, "Reached maximum number of iterations when fitting star inside Roche Lobe in RLOF" }}, + { ERROR::TOO_MANY_TIMESTEPS_IN_TIMESTEPS_FILE, { ERROR_SCOPE::ALWAYS, "Number of timesteps in timestpes file exceeds maximum timesteps" }}, { ERROR::UNEXPECTED_END_OF_FILE, { ERROR_SCOPE::ALWAYS, "Unexpected end of file" }}, { ERROR::UNEXPECTED_LOG_FILE_TYPE, { ERROR_SCOPE::ALWAYS, "Unexpected log file type" }}, { ERROR::UNEXPECTED_SN_EVENT, { ERROR_SCOPE::ALWAYS, "Unexpected supernova event in this context" }}, @@ -813,6 +832,9 @@ enum class EVOLUTION_STATUS: int { ERROR, TIMES_UP, STEPS_UP, + NO_TIMESTEPS, + TIMESTEPS_EXHAUSTED, + TIMESTEPS_NOT_CONSUMED, SSE_ERROR, BINARY_ERROR, DCO_MERGER_TIME, @@ -832,6 +854,9 @@ const COMPASUnorderedMap EVOLUTION_STATUS_LABEL = { EVOLUTION_STATUS::ERROR, "Evolution stopped because an error occurred" }, { EVOLUTION_STATUS::TIMES_UP, "Allowed time exceeded" }, { EVOLUTION_STATUS::STEPS_UP, "Allowed timesteps exceeded" }, + { EVOLUTION_STATUS::NO_TIMESTEPS, "No user-provided timesteps read" }, + { EVOLUTION_STATUS::TIMESTEPS_EXHAUSTED, "User-provided timesteps exhausted" }, + { EVOLUTION_STATUS::TIMESTEPS_NOT_CONSUMED, "User-provided timesteps not consumed" }, { EVOLUTION_STATUS::SSE_ERROR, "SSE error for one of the constituent stars" }, { EVOLUTION_STATUS::BINARY_ERROR, "Error evolving binary" }, { EVOLUTION_STATUS::DCO_MERGER_TIME, "Time exceeded DCO merger (formation + coalescence) time" }, @@ -2808,8 +2833,8 @@ class VariantPropertyType: public boost::static_visitor { // data type is, as the name suggest, the data type of the property, from the TYPENAME enum class (above) // header string is the string to be printed as the column header for the property // units string is the string to be printed as the units header for the property -// field width is the printf() field width of the property (meaning varies per the data type, refer to printf() documentation) -// field precision is the printf() field precision of the property (meaning varies per the data type, refer to printf() documentation) +// field width is the printf() field width of the property (does not apply to HDF5 files; meaning varies per the data type, refer to printf() documentation) +// field precision is the printf() field precision of the property (does not apply to HDF5 files; meaning varies per the data type, refer to printf() documentation) typedef std::tuple PROPERTY_DETAILS; @@ -2822,7 +2847,7 @@ typedef std::tuple PROPERTY_DETAIL // the logfiles - all keys present here should also be in the STAR_PROPERTIES #define and // STAR_PROPERTIES_LABEL const std::map ANY_STAR_PROPERTY_DETAIL = { - { ANY_STAR_PROPERTY::AGE, { TYPENAME::DOUBLE, "Age", "Myr", 16, 8 }}, + { ANY_STAR_PROPERTY::AGE, { TYPENAME::DOUBLE, "Age", "Myr", 24, 15}}, { ANY_STAR_PROPERTY::ANGULAR_MOMENTUM, { TYPENAME::DOUBLE, "Ang_Momentum", "Msol AU^2 yr^-1", 14, 6 }}, { ANY_STAR_PROPERTY::BINDING_ENERGY_AT_COMMON_ENVELOPE, { TYPENAME::DOUBLE, "Binding_Energy@CE", "erg", 14, 6 }}, { ANY_STAR_PROPERTY::BINDING_ENERGY_FIXED, { TYPENAME::DOUBLE, "BE_Fixed", "erg", 14, 6 }}, @@ -2840,10 +2865,10 @@ const std::map ANY_STAR_PROPERTY_DETAIL = { { ANY_STAR_PROPERTY::CORE_MASS_AT_COMPACT_OBJECT_FORMATION, { TYPENAME::DOUBLE, "Mass_Core@CO", "Msol", 14, 6 }}, { ANY_STAR_PROPERTY::DRAWN_KICK_MAGNITUDE, { TYPENAME::DOUBLE, "Drawn_Kick_Magnitude", "kms^-1", 14, 6 }}, { ANY_STAR_PROPERTY::DOMINANT_MASS_LOSS_RATE, { TYPENAME::INT, "Dominant_Mass_Loss_Rate", "-", 4, 1 }}, - { ANY_STAR_PROPERTY::DT, { TYPENAME::DOUBLE, "dT", "Myr", 16, 8 }}, - { ANY_STAR_PROPERTY::DYNAMICAL_TIMESCALE, { TYPENAME::DOUBLE, "Tau_Dynamical", "Myr", 16, 8 }}, - { ANY_STAR_PROPERTY::DYNAMICAL_TIMESCALE_POST_COMMON_ENVELOPE, { TYPENAME::DOUBLE, "Tau_Dynamical>CE", "Myr", 16, 8 }}, - { ANY_STAR_PROPERTY::DYNAMICAL_TIMESCALE_PRE_COMMON_ENVELOPE, { TYPENAME::DOUBLE, "Tau_DynamicalCE", "Myr", 24, 15}}, + { ANY_STAR_PROPERTY::DYNAMICAL_TIMESCALE_PRE_COMMON_ENVELOPE, { TYPENAME::DOUBLE, "Tau_Dynamical ANY_STAR_PROPERTY_DETAIL = { { ANY_STAR_PROPERTY::LAMBDA_LOVERIDGE_WINDS, { TYPENAME::DOUBLE, "Lambda_Loveridge_Winds", "-", 14, 6 }}, { ANY_STAR_PROPERTY::LAMBDA_NANJING, { TYPENAME::DOUBLE, "Lambda_Nanjing", "-", 14, 6 }}, { ANY_STAR_PROPERTY::LBV_PHASE_FLAG, { TYPENAME::BOOL, "LBV_Phase_Flag", "Event", 0, 0 }}, - { ANY_STAR_PROPERTY::LUMINOSITY, { TYPENAME::DOUBLE, "Luminosity", "Lsol", 14, 6 }}, + { ANY_STAR_PROPERTY::LUMINOSITY, { TYPENAME::DOUBLE, "Luminosity", "Lsol", 24, 15}}, { ANY_STAR_PROPERTY::LUMINOSITY_POST_COMMON_ENVELOPE, { TYPENAME::DOUBLE, "Luminosity>CE", "Lsol", 14, 6 }}, { ANY_STAR_PROPERTY::LUMINOSITY_PRE_COMMON_ENVELOPE, { TYPENAME::DOUBLE, "Luminosity ANY_STAR_PROPERTY_DETAIL = { { ANY_STAR_PROPERTY::METALLICITY, { TYPENAME::DOUBLE, "Metallicity@ZAMS", "-", 14, 6 }}, { ANY_STAR_PROPERTY::MOMENT_OF_INERTIA, { TYPENAME::DOUBLE, "Moment_Of_Inertia", "Msol Rsol^2", 14, 6 }}, { ANY_STAR_PROPERTY::MZAMS, { TYPENAME::DOUBLE, "Mass@ZAMS", "Msol", 14, 6 }}, - { ANY_STAR_PROPERTY::NUCLEAR_TIMESCALE, { TYPENAME::DOUBLE, "Tau_Nuclear", "Myr", 16, 8 }}, - { ANY_STAR_PROPERTY::NUCLEAR_TIMESCALE_POST_COMMON_ENVELOPE, { TYPENAME::DOUBLE, "Tau_Nuclear>CE", "Myr", 16, 8 }}, - { ANY_STAR_PROPERTY::NUCLEAR_TIMESCALE_PRE_COMMON_ENVELOPE, { TYPENAME::DOUBLE, "Tau_NuclearCE", "Myr", 24, 15}}, + { ANY_STAR_PROPERTY::NUCLEAR_TIMESCALE_PRE_COMMON_ENVELOPE, { TYPENAME::DOUBLE, "Tau_Nuclear ANY_STAR_PROPERTY_DETAIL = { { ANY_STAR_PROPERTY::PULSAR_BIRTH_SPIN_DOWN_RATE, { TYPENAME::DOUBLE, "Pulsar_Birth_Spin_Down", "s/s", 14, 6 }}, { ANY_STAR_PROPERTY::PULSAR_SPIN_FREQUENCY, { TYPENAME::DOUBLE, "Pulsar_Spin_Freq", "rad/s", 14, 6 }}, { ANY_STAR_PROPERTY::PULSAR_SPIN_PERIOD, { TYPENAME::DOUBLE, "Pulsar_Spin_Period", "ms", 14, 6 }}, - { ANY_STAR_PROPERTY::RADIAL_EXPANSION_TIMESCALE, { TYPENAME::DOUBLE, "Tau_Radial", "Myr", 16, 8 }}, - { ANY_STAR_PROPERTY::RADIAL_EXPANSION_TIMESCALE_POST_COMMON_ENVELOPE, { TYPENAME::DOUBLE, "Tau_Radial>CE", "Myr", 16, 8 }}, - { ANY_STAR_PROPERTY::RADIAL_EXPANSION_TIMESCALE_PRE_COMMON_ENVELOPE, { TYPENAME::DOUBLE, "Tau_RadialCE", "Myr", 24, 15}}, + { ANY_STAR_PROPERTY::RADIAL_EXPANSION_TIMESCALE_PRE_COMMON_ENVELOPE, { TYPENAME::DOUBLE, "Tau_RadialNS", "Event", 0, 0 }}, @@ -2934,11 +2959,11 @@ const std::map ANY_STAR_PROPERTY_DETAIL = { { ANY_STAR_PROPERTY::TEMPERATURE, { TYPENAME::DOUBLE, "Teff", "K", 14, 6 }}, { ANY_STAR_PROPERTY::TEMPERATURE_POST_COMMON_ENVELOPE, { TYPENAME::DOUBLE, "Teff>CE", "K", 14, 6 }}, { ANY_STAR_PROPERTY::TEMPERATURE_PRE_COMMON_ENVELOPE, { TYPENAME::DOUBLE, "TeffCE", "Myr", 16, 8 }}, - { ANY_STAR_PROPERTY::THERMAL_TIMESCALE_PRE_COMMON_ENVELOPE, { TYPENAME::DOUBLE, "Tau_ThermalCE", "Myr", 24, 15}}, + { ANY_STAR_PROPERTY::THERMAL_TIMESCALE_PRE_COMMON_ENVELOPE, { TYPENAME::DOUBLE, "Tau_Thermal BINARY_PROPERTY_DETAIL = { { BINARY_PROPERTY::BE_BINARY_CURRENT_COMPANION_MASS, { TYPENAME::DOUBLE, "Companion_Mass", "Msol", 14, 6 }}, { BINARY_PROPERTY::BE_BINARY_CURRENT_COMPANION_RADIUS, { TYPENAME::DOUBLE, "Companion_Radius", "Rsol", 14, 6 }}, { BINARY_PROPERTY::BE_BINARY_CURRENT_COMPANION_TEFF, { TYPENAME::DOUBLE, "Companion_Teff", "K", 14, 6 }}, - { BINARY_PROPERTY::BE_BINARY_CURRENT_DT, { TYPENAME::DOUBLE, "dT", "Myr", 16, 8 }}, + { BINARY_PROPERTY::BE_BINARY_CURRENT_DT, { TYPENAME::DOUBLE, "dT", "Myr", 24, 15}}, { BINARY_PROPERTY::BE_BINARY_CURRENT_ECCENTRICITY, { TYPENAME::DOUBLE, "Eccentricity", "-", 14, 6 }}, { BINARY_PROPERTY::BE_BINARY_CURRENT_ID, { TYPENAME::OBJECT_ID, "ID", "-", 12, 1 }}, { BINARY_PROPERTY::BE_BINARY_CURRENT_NS_MASS, { TYPENAME::DOUBLE, "NS_Mass", "Msol", 14, 6 }}, { BINARY_PROPERTY::BE_BINARY_CURRENT_SEMI_MAJOR_AXIS, { TYPENAME::DOUBLE, "SemiMajorAxis", "Rsol", 14, 6 }}, - { BINARY_PROPERTY::BE_BINARY_CURRENT_TOTAL_TIME, { TYPENAME::DOUBLE, "Total_Time", "Myr", 16, 8 }}, - { BINARY_PROPERTY::CIRCULARIZATION_TIMESCALE, { TYPENAME::DOUBLE, "Tau_Circ", "Myr", 16, 8 }}, + { BINARY_PROPERTY::BE_BINARY_CURRENT_TOTAL_TIME, { TYPENAME::DOUBLE, "Total_Time", "Myr", 24, 15}}, + { BINARY_PROPERTY::CIRCULARIZATION_TIMESCALE, { TYPENAME::DOUBLE, "Tau_Circ", "Myr", 24, 15}}, { BINARY_PROPERTY::COMMON_ENVELOPE_AT_LEAST_ONCE, { TYPENAME::BOOL, "CEE", "Event", 0, 0 }}, { BINARY_PROPERTY::COMMON_ENVELOPE_EVENT_COUNT, { TYPENAME::UINT, "CE_Event_Counter", "Count", 6, 1 }}, { BINARY_PROPERTY::DIMENSIONLESS_KICK_MAGNITUDE, { TYPENAME::DOUBLE, "Kick_Magnitude(uK)", "-", 14, 6 }}, { BINARY_PROPERTY::DOUBLE_CORE_COMMON_ENVELOPE, { TYPENAME::BOOL, "Double_Core_CE", "Event", 0, 0 }}, - { BINARY_PROPERTY::DT, { TYPENAME::DOUBLE, "dT", "Myr", 16, 8 }}, + { BINARY_PROPERTY::DT, { TYPENAME::DOUBLE, "dT", "Myr", 24, 15}}, { BINARY_PROPERTY::ECCENTRICITY, { TYPENAME::DOUBLE, "Eccentricity", "-", 14, 6 }}, { BINARY_PROPERTY::ECCENTRICITY_AT_DCO_FORMATION, { TYPENAME::DOUBLE, "Eccentricity@DCO", "-", 14, 6 }}, { BINARY_PROPERTY::ECCENTRICITY_INITIAL, { TYPENAME::DOUBLE, "Eccentricity@ZAMS", "-", 14, 6 }}, @@ -3030,8 +3055,8 @@ const std::map BINARY_PROPERTY_DETAIL = { { BINARY_PROPERTY::RLOF_PRE_STEP_STAR_TO_ROCHE_LOBE_RADIUS_RATIO_1, { TYPENAME::DOUBLE, "Radius(1)|RLCE", "Event", 0, 0 }}, - { BINARY_PROPERTY::RLOF_TIME_POST_MT, { TYPENAME::DOUBLE, "Time>MT", "Myr", 16, 8 }}, - { BINARY_PROPERTY::RLOF_TIME_PRE_MT, { TYPENAME::DOUBLE, "TimeMT", "Myr", 24, 15}}, + { BINARY_PROPERTY::RLOF_TIME_PRE_MT, { TYPENAME::DOUBLE, "TimeCE", "Rsol", 14, 6 }}, { BINARY_PROPERTY::ROCHE_LOBE_RADIUS_1_PRE_COMMON_ENVELOPE, { TYPENAME::DOUBLE, "RocheLobe(1) BINARY_PROPERTY_DETAIL = { { BINARY_PROPERTY::SUPERNOVA_ORBIT_INCLINATION_VECTOR_Y, { TYPENAME::DOUBLE, "Orbital_AM_Vector>SN_Y", "-", 14, 6 }}, { BINARY_PROPERTY::SUPERNOVA_ORBIT_INCLINATION_VECTOR_Z, { TYPENAME::DOUBLE, "Orbital_AM_Vector>SN_Z", "-", 14, 6 }}, { BINARY_PROPERTY::SUPERNOVA_STATE, { TYPENAME::SN_STATE, "Supernova_State", "State", 4, 1 }}, - { BINARY_PROPERTY::SYNCHRONIZATION_TIMESCALE, { TYPENAME::DOUBLE, "Tau_Sync", "Myr", 16, 8 }}, + { BINARY_PROPERTY::SYNCHRONIZATION_TIMESCALE, { TYPENAME::DOUBLE, "Tau_Sync", "Myr", 24, 15}}, { BINARY_PROPERTY::SYSTEMIC_SPEED, { TYPENAME::DOUBLE, "SystemicSpeed", "kms^-1", 14, 6 }}, - { BINARY_PROPERTY::TIME, { TYPENAME::DOUBLE, "Time", "Myr", 16, 8 }}, - { BINARY_PROPERTY::TIME_TO_COALESCENCE, { TYPENAME::DOUBLE, "Coalescence_Time", "Myr", 16, 8 }}, + { BINARY_PROPERTY::TIME, { TYPENAME::DOUBLE, "Time", "Myr", 24, 15}}, + { BINARY_PROPERTY::TIME_TO_COALESCENCE, { TYPENAME::DOUBLE, "Coalescence_Time", "Myr", 24, 15}}, { BINARY_PROPERTY::TOTAL_ANGULAR_MOMENTUM, { TYPENAME::DOUBLE, "Ang_Momentum_Total", "Msol AU^2 yr^-1", 14, 6 }}, { BINARY_PROPERTY::TOTAL_ENERGY, { TYPENAME::DOUBLE, "Energy_Total", "Msol AU^2 yr^-2", 14, 6 }}, { BINARY_PROPERTY::UNBOUND, { TYPENAME::BOOL, "Unbound", "State", 0, 0 }}, @@ -3176,7 +3201,7 @@ const std::map PROGRAM_OPTION_DETAIL = { { PROGRAM_OPTION::MASS_RATIO_DISTRIBUTION_MAX, { TYPENAME::DOUBLE, "Mass_Ratio_Dstrbtn_Max", "-", 14, 6 }}, { PROGRAM_OPTION::MASS_RATIO_DISTRIBUTION_MIN, { TYPENAME::DOUBLE, "Mass_Ratio_Dstrbtn_Min", "-", 14, 6 }}, - { PROGRAM_OPTION::MAXIMUM_EVOLUTION_TIME, { TYPENAME::DOUBLE, "Max_Evolution_Time", "Myr", 14, 6 }}, + { PROGRAM_OPTION::MAXIMUM_EVOLUTION_TIME, { TYPENAME::DOUBLE, "Max_Evolution_Time", "Myr", 24, 15}}, { PROGRAM_OPTION::MAXIMUM_DONOR_MASS, { TYPENAME::DOUBLE, "Max_Donor_Mass", "Msol", 14, 6 }}, { PROGRAM_OPTION::MAXIMUM_NEUTRON_STAR_MASS, { TYPENAME::DOUBLE, "Max_NS_Mass", "Msol", 14, 6 }}, { PROGRAM_OPTION::MAXIMUM_TIMESTEPS, { TYPENAME::INT, "Max_Timesteps", "Count", 10, 1 }}, diff --git a/src/utils.cpp b/src/utils.cpp index 9c6507d04..a12dda643 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -1438,4 +1438,109 @@ namespace utils { return splashString; // return the splash string } + + /* + * Read timesteps from timesteps file + * + * Timesteps file is expected to be an ascii file with one timestep per record. + * Timesteps must be > 0.0 + * + * + * std::tuple ReadTimesteps(const std::string p_TimestepsFileName) + * + * @param [IN] p_TimestepsFileName Filename to be read - should be fully qualified + * @return Tuple containing error value and timesteps vector + * The error value returned will be: + * ERROR::NONE if no error occurred + * ERROR::EMPTY_FILENAME if the filename provided was an empty string + * ERROR::FILE_DOES_NOT_EXIST if the timesteps file does not exist + * ERROR::FILE_OPEN_ERROR if the timesteps file exists but could not be opened + * ERROR::FILE_READ_ERROR if the timesteps file could not be read + * ERROR::EMPTY_FILE if the timesteps file contains no content + * ERROR::INVALID_VALUE_IN_TIMESTEPS_FILE if the file contains an invalid value for timestep + * ERROR::TOO_MANY_TIMESTEPS_IN_TIMESTEPS_FILE if the file contains too many timesteps (> maximum per OPTIONS) + * If the error returned is not ERROR:NONE, the content of the timesteps vector returned is not defined + */ + std::tuple ReadTimesteps(const std::string p_TimestepsFileName) { + + ERROR error = ERROR::NONE; // error - initially NONE + + DBL_VECTOR timesteps; // timesteps vector + + if (p_TimestepsFileName.empty()) { // timesteps filename empty? + error = ERROR::EMPTY_FILENAME; // yes - fail + } + else { + + if (!FileExists(p_TimestepsFileName)) { // timesteps file exists? + error = ERROR::FILE_DOES_NOT_EXIST; // no - fail + } + else { // yes + std::ifstream timestepsFile; + timestepsFile.exceptions(std::ifstream::failbit | std::ifstream::badbit); + try { + timestepsFile.open(p_TimestepsFileName); // open the timesteps file + if (!timestepsFile.is_open()) { // open ok? + error = ERROR::FILE_OPEN_ERROR; // no - fail + } + else { // yes - file open + std::string rec; // record read from file + unsigned int numTimesteps = 0; + while (std::getline(timestepsFile, rec)) { // get next record from timesteps file + + if (rec.size() > 0 && (rec[rec.size() - 1] == '\n' || rec[rec.size() - 1] == '\r')) { // last character `\n` or `\r`? + rec.erase(rec.size() - 1); // yes - strip it + } + + if (!rec.empty()) { // blank record? + try { // no - process it + size_t lastChar; + long double v = std::stold(rec, &lastChar); // try conversion + if (lastChar != (rec.size())) { // conversion valid only if rec completely consumed + error = ERROR::INVALID_VALUE_IN_TIMESTEPS_FILE; // not a valid DOUBLE + break; // stop processing + } + + timesteps.push_back(v); // add timestep to timesteps vector + + numTimesteps++; // increment number of timesteps read + if (numTimesteps >= ABSOLUTE_MAXIMUM_TIMESTEPS) { // number of timesteps exceeds maximum? + error = ERROR::TOO_MANY_TIMESTEPS_IN_TIMESTEPS_FILE; // yes - fail + break; // stop processing + } + } + catch (const std::out_of_range& e) { // conversion failed + error = ERROR::INVALID_VALUE_IN_TIMESTEPS_FILE; // not a valid DOUBLE + break; // stop processing + } + catch (const std::invalid_argument& e) { // conversion failed + error = ERROR::INVALID_VALUE_IN_TIMESTEPS_FILE; // not a valid DOUBLE + break; // stop processing + } + } + } + try { + timestepsFile.close(); // close the timesteps file + } + catch (std::ifstream::failure& e) { // close failed + error = ERROR::FILE_NOT_CLOSED; // fail + } + } + } + catch (std::ifstream::failure& e) { // something was flagged... + if (timestepsFile.eof()) { // end-of-file? + if (timesteps.size() < 1) { // yes - at least one timestep read? + error = ERROR::EMPTY_FILE; // no - fail + } + } + else { // not end-of-file - error + error = ERROR::FILE_READ_ERROR; // fail + } + } + + } + } + return std::make_tuple(error, timesteps); + } + } diff --git a/src/utils.h b/src/utils.h index 459e78e26..1ec7bf997 100755 --- a/src/utils.h +++ b/src/utils.h @@ -138,6 +138,8 @@ namespace utils { std::string SplashScreen(const bool p_Print = true); + std::tuple ReadTimesteps(const std::string p_TimestepsFileName); + } #endif // __utils_h__ diff --git a/src/vector3d.h b/src/vector3d.h index c0f6594b9..78bb291c3 100644 --- a/src/vector3d.h +++ b/src/vector3d.h @@ -44,7 +44,6 @@ class Vector3d { const double p_PsiE); Vector3d UnitVector(); - /////////////////////////////// // Overload operators // Overload Indexing operator diff --git a/src/yaml.h b/src/yaml.h index 83a7def68..69fd6ea76 100644 --- a/src/yaml.h +++ b/src/yaml.h @@ -276,6 +276,7 @@ namespace yaml { " --notes", " --notes-hdrs", " --output-container", + " --timesteps-filename", "", " ### STELLAR PROPERTIES", " --chemically-homogeneous-evolution # chemically homogeneous evolution", From f2914dfc97616baae7a008b738e4cdfd73af6dfe Mon Sep 17 00:00:00 2001 From: Jeff Riley Date: Fri, 5 Jan 2024 14:48:52 +1100 Subject: [PATCH 2/4] Updates --- src/BaseStar.cpp | 2 +- src/GiantBranch.cpp | 1 - src/Star.cpp | 6 +++--- 3 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index 195028f5b..c440196be 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -2680,7 +2680,7 @@ double BaseStar::CalculateMassLossRate() { break; default: // unknown mass-loss prescription - SHOW_WARN(ERROR::UNKNOWN_MASS_LOSS_PRESCRIPTION, "DEFAULT!"); // show warning + SHOW_WARN(ERROR::UNKNOWN_MASS_LOSS_PRESCRIPTION, "Using HURLEY"); // show warning LBVRate = CalculateMassLossRateLBV(LBV_PRESCRIPTION::HURLEY_ADD); otherWindsRate = CalculateMassLossRateHurley(); // use HURLEY if (utils::Compare(LBVRate, otherWindsRate) > 0) { diff --git a/src/GiantBranch.cpp b/src/GiantBranch.cpp index a3b2565b2..936d76f4d 100644 --- a/src/GiantBranch.cpp +++ b/src/GiantBranch.cpp @@ -1907,7 +1907,6 @@ STELLAR_TYPE GiantBranch::ResolveTypeIIaSN() { */ STELLAR_TYPE GiantBranch::ResolvePairInstabilitySN() { - m_Mass = 0.0; m_Luminosity = 0.0; m_Radius = 0.0; m_Temperature = 0.0; diff --git a/src/Star.cpp b/src/Star.cpp index 2801db625..12a41a5d3 100644 --- a/src/Star.cpp +++ b/src/Star.cpp @@ -330,8 +330,8 @@ STELLAR_TYPE Star::UpdateAttributesAndAgeOneTimestep(const double p_DeltaMass, * * * Checks whether the star: - * - is a massless remnant (checked after applying p_DeltaMass and p_DeltaMass0) <--- this is no longer ture (it may have been once) - * - has become a supernova (checked after applying p_DeltaMass and p_DeltaMass0) <--- this is no longer ture (it may have been once) JR: this causes different behaviour SSE vs BSE + * - is a massless remnant (checked after applying p_DeltaMass and p_DeltaMass0) <--- JR: this is no longer true (it may have been once) + * - has become a supernova (checked after applying p_DeltaMass and p_DeltaMass0) <--- JR this is no longer true (it may have been once) This causes different behaviour SSE vs BSE * - should skip this phase for this timestep (checked after applying p_DeltaMass and p_DeltaMass0) * * If none of the above are true the star's attributes are updated based on the mass changes required and the star's @@ -422,7 +422,7 @@ double Star::EvolveOneTimestep(const double p_Dt, const bool p_Force) { SaveState(); // save the state of the star - in case we want to revert - double minTimestep = std::max(m_Star->CalculateDynamicalTimescale(), NUCLEAR_MINIMUM_TIMESTEP); // calculate the minimum timestep - maximum of dynamical timescale for this star and the absolute minimum timestep + double minTimestep = std::max(m_Star->CalculateDynamicalTimescale(), ABSOLUTE_MINIMUM_TIMESTEP); // calculate the minimum timestep - maximum of dynamical timescale for this star and the absolute minimum timestep // evolve the star a single timestep From 1fb0796c0e2853ed43a5ea5e8acdd59af1e427fd Mon Sep 17 00:00:00 2001 From: Jeff Riley Date: Sat, 6 Jan 2024 17:06:51 +1100 Subject: [PATCH 3/4] Updates --- src/BaseBinaryStar.cpp | 2 +- src/BaseStar.cpp | 17 ++++++++--------- src/BaseStar.h | 2 +- src/CHeB.cpp | 4 ++-- 4 files changed, 12 insertions(+), 13 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 4ca3f701a..74d6b6480 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1453,7 +1453,7 @@ void BaseBinaryStar::EvaluateSupernovae() { * void ResolveCommonEnvelopeEvent() */ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { - + double alphaCE = OPTIONS->CommonEnvelopeAlpha(); // CE efficiency parameter double eccentricity = Eccentricity(); // current eccentricity (before CEE) diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index c440196be..43d60a937 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -301,7 +301,7 @@ COMPAS_VARIABLE BaseStar::StellarPropertyValue(const T_ANY_PROPERTY p_Property) case ANY_STAR_PROPERTY::DT: value = Dt(); break; case ANY_STAR_PROPERTY::DYNAMICAL_TIMESCALE: value = CalculateDynamicalTimescale(); break; case ANY_STAR_PROPERTY::ECCENTRIC_ANOMALY: value = SN_EccentricAnomaly(); break; - case ANY_STAR_PROPERTY::ENV_MASS: value = Mass()-CoreMass(); break; + case ANY_STAR_PROPERTY::ENV_MASS: value = Mass() - CoreMass(); break; case ANY_STAR_PROPERTY::ERROR: value = Error(); break; case ANY_STAR_PROPERTY::EVOL_STATUS: value = EvolutionStatus(); break; case ANY_STAR_PROPERTY::EXPERIENCED_AIC: value = ExperiencedAIC(); break; @@ -2611,7 +2611,7 @@ double BaseStar::CalculateMassLossRateFlexible2023() { OPTIONS->LuminousBlueVariablePrescription() == LBV_PRESCRIPTION::HURLEY_ADD ) { // check whether we should add other winds to the LBV winds (always for HURLEY_ADD prescription, only if not in LBV regime for others) if ((utils::Compare(teff, RSG_MAXIMUM_TEMP) < 0) && - (utils::Compare(m_MZAMS, 8.0) >= 0) && // JR: where does this 8.0 come from? We shouldn't have undocumented magic numbers in the code... If it's an accepted threshold, let's document it and give it a name. + (utils::Compare(m_MZAMS, 8.0) >= 0) && // JR: We shouldn't have undocumented magic numbers in the code... This is an accepted threshold, so let's document it and give it a name. IsOneOf(GIANTS)) { // RSG criteria, below 8kK, above 8Msol, and core helium burning giant(CHeB, FGB, EAGB, TPAGB) otherWindsRate = CalculateMassLossRateRSG(OPTIONS->RSGMassLoss()); m_DominantMassLossRate = MASS_LOSS_TYPE::RSG; @@ -3999,14 +3999,12 @@ STELLAR_TYPE BaseStar::UpdateAttributesAndAgeOneTimestep(const double p_DeltaMas utils::Compare(p_DeltaMass0, 0.0) != 0 || // mass0 change? or... p_DeltaTime > 0) { // age/time advance? (don't use utils::Compare() here) // yes - update attributes - if (p_DeltaTime > 0.0) AgeOneTimestepPreamble(p_DeltaTime); // advance dt, age, simulation time if necessary (don't use utils::Compare() here) + AgeOneTimestepPreamble(p_DeltaTime); // advance dt, age, simulation time if necessary (don't use utils::Compare() here) if (ShouldSkipPhase()) stellarType = ResolveSkippedPhase(); // skip phase if required else { // not skipped - execute phase - if (p_DeltaTime > 0.0) {; // evolve on phase if necessary (don't use utils::Compare() here) - stellarType = EvolveOnPhase(); // evolve on phase - if (stellarType == m_StellarType) stellarType = ResolveEndOfPhase();// check for need to move off phase - } + stellarType = EvolveOnPhase(p_DeltaTime); // evolve on phase + if (stellarType == m_StellarType) stellarType = ResolveEndOfPhase();// check for need to move off phase } } } @@ -4087,12 +4085,11 @@ void BaseStar::UpdateMassTransferDonorHistory() { * * @return Stellar Type to which star should evolve - unchanged if not moving off current phase */ -STELLAR_TYPE BaseStar::EvolveOnPhase() { +STELLAR_TYPE BaseStar::EvolveOnPhase(const double p_DeltaTime) { STELLAR_TYPE stellarType = m_StellarType; if (ShouldEvolveOnPhase()) { // evolve timestep on phase - m_Tau = CalculateTauOnPhase(); m_COCoreMass = CalculateCOCoreMassOnPhase(); @@ -4109,10 +4106,12 @@ STELLAR_TYPE BaseStar::EvolveOnPhase() { m_Temperature = CalculateTemperatureOnPhase(); + if (p_DeltaTime > 0.0) { STELLAR_TYPE thisStellarType = ResolveEnvelopeLoss(); // resolve envelope loss if it occurs - possibly new stellar type if (thisStellarType != m_StellarType) { // thisStellarType overrides stellarType (from CalculateRadiusAndStellarTypeOnPhase()) stellarType = thisStellarType; } + } } return stellarType; diff --git a/src/BaseStar.h b/src/BaseStar.h index 3b1483e39..5564c3b33 100644 --- a/src/BaseStar.h +++ b/src/BaseStar.h @@ -582,7 +582,7 @@ class BaseStar { virtual void EvolveOneTimestepPreamble() { }; // Default is NO-OP - STELLAR_TYPE EvolveOnPhase(); + STELLAR_TYPE EvolveOnPhase(const double p_DeltaTime); virtual STELLAR_TYPE EvolveToNextPhase() { return m_StellarType; } diff --git a/src/CHeB.cpp b/src/CHeB.cpp index 025f51b05..6d8662727 100755 --- a/src/CHeB.cpp +++ b/src/CHeB.cpp @@ -100,7 +100,7 @@ double CHeB::CalculateLambdaDewi() const { double lambdaCE; if (utils::Compare(envMass, 1.0) >= 0) lambdaCE = 2.0 * lambda1; // (A.1) Bottom, Claeys+2014 - else if (utils::Compare(envMass, 0.0) > 0) lambdaCE = 2.0 * (lambda2 + (std::sqrt(envMass) * (lambda1 - lambda2))); // (A.1) Mid, Claeys+2014 + else if (utils::Compare(envMass, 0.0) > 0) lambdaCE = 2.0 * (lambda2 + (std::sqrt(envMass) * (lambda1 - lambda2))); // (A.1) Mid, Claeys+2014 else lambdaCE = 2.0 * lambda2; // (A.1) Top, Claeys+2014 return lambdaCE; @@ -1214,7 +1214,7 @@ double CHeB::CalculateRemnantRadius() const { * @return Core mass on the First Giant Branch in Msol */ double CHeB::CalculateCoreMassOnPhase(const double p_Mass, const double p_Tau) const { - return std::min(((1.0 - p_Tau) * CalculateCoreMassAtHeIgnition(p_Mass)) + (p_Tau * CalculateCoreMassAtBAGB(p_Mass)), m_Mass); //He mass capped at total mass (should become HeMS star) + return std::min(((1.0 - p_Tau) * CalculateCoreMassAtHeIgnition(p_Mass)) + (p_Tau * CalculateCoreMassAtBAGB(p_Mass)), m_Mass); // We mass capped at total mass (should become HeMS star) } From 00ceff94735a4ae07e46039e42a09c3a4d2a7d49 Mon Sep 17 00:00:00 2001 From: Jeff Riley Date: Wed, 10 Jan 2024 10:11:40 +1100 Subject: [PATCH 4/4] Update documention; allow comment lines in timesteps files --- .../program-options-list-defaults.rst | 6 +++- .../pages/User guide/timestep-files.rst | 28 +++++++++++++++++++ online-docs/pages/User guide/user-guide.rst | 1 + online-docs/pages/whats-new.rst | 6 ++++ src/changelog.h | 7 ++++- src/utils.cpp | 12 ++++++-- 6 files changed, 56 insertions(+), 4 deletions(-) create mode 100644 online-docs/pages/User guide/timestep-files.rst 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 ed88ee079..7235db69e 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 @@ -413,7 +413,7 @@ Default = 5.75 :ref:`Back to Top ` **--grid** |br| -Grid filename. |br| +Grid filename. (See :doc:`Grid files <../grid-files>`) |br| Default = ’’ (None) **--grid-lines-to-process** |br| @@ -1161,6 +1161,10 @@ Default = FALSE :ref:`Back to Top ` +**--timestep-filename** |br| +User-defined timesteps filename. (See :doc:`Timestep files <../timestep-files>`) |br| +Default = ’’ (None) + **--timestep-multiplier** |br| Multiplicative factor for timestep duration. |br| Default = 1.0 diff --git a/online-docs/pages/User guide/timestep-files.rst b/online-docs/pages/User guide/timestep-files.rst new file mode 100644 index 000000000..a97af938a --- /dev/null +++ b/online-docs/pages/User guide/timestep-files.rst @@ -0,0 +1,28 @@ +Timestep files +============== + +A timestep file allows users to specify the timesteps to be taken during Single Star Evolution (SSE) or Binary Star Evolution (BSE). + +If a timestep file is provided (via the ``timesteps-filename`` program option), the COMPAS code uses the user-provided timeteps to +evolve stars/systems instead of calculating a new timestep at each iteration. The COMPAS code will read the first non-blank and +non-comment line of the timesteps file and set the first timestep to that value, then for each iteration during evolution, each +subsequent non-blank and non-comment line of the timesteps file is read and the timestep for that iteration set to the value read. + +Note that COMPAS will use the timestep exactly as read - the timesteps read from the timesteps file are not quantised by COMPAS, and +neither will they be truncated to any limit (e.g. ``ABSOLUTE_MINIMUM_TIMESTEP`` from ``constants.h``) or multiplied by any timestep +multiplier set by the ``timestep-multiplier`` option. Furthermore, no check will be made after the timestep is taken to limit the +radial expansion of the star being evolved. + +Timesteps files can have blank lines and comments. Comments begin with a hash/pound character ('#') - the hash/pound character and text +beyond it are ignored by COMPAS. + +The ``timesteps-filename`` option is valid both on the command line and in a grid file. If the ``timesteps-filename`` option is specified +in conjunction with the ``--number-of-systems`` (``-n``) option, or with option ranges and/or sets, the same timeteps file specified will +be used for each star or system being evolved. On the other hand, if the ``timesteps-filename`` option is specified in a grid file, the +timesteps file specified on a grid file line will be read and used for the star or system evolved by that grid file line. + +Each line of a timesteps file which is not a blank line or comment line must contain a single, non-negative, floating-point number (in +plain text - COMPAS will read the plain text and convert it to a floating-point number). + +COMPAS imposes a hard limit of ``1,000,000`` timesteps in a timesteps file. + diff --git a/online-docs/pages/User guide/user-guide.rst b/online-docs/pages/User guide/user-guide.rst index 5f40c53cf..1f0cbe64b 100644 --- a/online-docs/pages/User guide/user-guide.rst +++ b/online-docs/pages/User guide/user-guide.rst @@ -10,6 +10,7 @@ This section contains the basic user guide for COMPAS. ./configuration ./Program options/program-options ./grid-files + ./timestep-files ./random-seed ./Running COMPAS/running-compas ./COMPAS output/output diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index ee4360849..13e4d7952 100644 --- a/online-docs/pages/whats-new.rst +++ b/online-docs/pages/whats-new.rst @@ -6,6 +6,12 @@ Following is a brief list of important updates to the COMPAS code. A complete r **LATEST RELEASE** |br| +**02.42.00 Jan 04, 2023** + +* Timesteps are now quantised to an integral multiple of 1e-12Myr. +* New option provided to allow user-defined timesteps: ``--timesteps-filename`` (See :doc:`Timestep files <../"User guide"/timestep-files>`). +* Code changes to make SSE and BSE evolution more consistent (See `PR 1052 `_). + **02.41.03 Dec 28, 2023** * The functions ``BaseBinaryStar::CalculateAngularMomentum()``, ``BaseBinaryStar::CalculateTotalEnergy()``, and ``BaseStar::AngularMomentum()`` changed to use moment of inertia instead of gyration radius. diff --git a/src/changelog.h b/src/changelog.h index 1a93b6ea2..8a7e73446 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1088,7 +1088,12 @@ // - Change to functionality (noted above) noted in 'What's New' online documentation page // 02.41.04 JR - Dec 30, 2023 - Defect repair: // - Fix for issue #1048 -// 02.42.99 JR - Jan 04. 2024 - Enhancements etc. TBC +// 02.42.99 JR - Jan 04. 2024 - Enhancements, defect repair, a little cleanup +// - added `timesteps-filename` option to allow users to provide preset timesteps for both SSE and BSE +// - updated documentation for new option; updated `What's New` +// - SSE vs BSE consistency: modified SSE to evolve a single star exactly as the primary in a wide binary with small companion +// - quantised timesteps to an integral multiple of 1E-12Myr - new constant `TIMESTEP_QUANTUM` in constants.h +// - little bit of code cleanup const std::string VERSION_STRING = "02.42.00"; diff --git a/src/utils.cpp b/src/utils.cpp index a12dda643..c0eb2a63e 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -1492,7 +1492,9 @@ namespace utils { rec.erase(rec.size() - 1); // yes - strip it } - if (!rec.empty()) { // blank record? + rec = trim(rec); // remove leading and trailing blanks + + if (!rec.empty() or rec[0] == '#') { // blank record or comment? try { // no - process it size_t lastChar; long double v = std::stold(rec, &lastChar); // try conversion @@ -1501,7 +1503,13 @@ namespace utils { break; // stop processing } - timesteps.push_back(v); // add timestep to timesteps vector + if (v < 0.0) { // timestep must be >= 0.0 + error = ERROR::INVALID_VALUE_IN_TIMESTEPS_FILE; // not a valid timestep + break; // stop processing + } + else { // ok - timestep >= 0.0 + timesteps.push_back(v); // add timestep to timesteps vector + } numTimesteps++; // increment number of timesteps read if (numTimesteps >= ABSOLUTE_MAXIMUM_TIMESTEPS) { // number of timesteps exceeds maximum?