diff --git a/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-specification-stellar.rst b/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-specification-stellar.rst index 87a61857e..b106ab1a9 100644 --- a/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-specification-stellar.rst +++ b/online-docs/pages/User guide/COMPAS output/standard-logfiles-record-specification-stellar.rst @@ -2050,6 +2050,22 @@ but not both. If both are printed then the file will contain two columns with th * - Header Strings: - Teff` **--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 @@ -1324,7 +1328,7 @@ Go to :ref:`the top of this page ` for the full alphabetical **Administrative** --mode, --number-of-systems, --evolve-double-white-dwarfs, --evolve-pulsars, --evolve-unbound-systems, --maximum-evolution-time, --maximum-number-timestep-iterations, ---random-seed, --timestep-multiplier +--random-seed, --timestep-multiplier, --timestep-filename --grid, --grid-start-line, --grid-lines-to-process 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..d0e5e172b --- /dev/null +++ b/online-docs/pages/User guide/timestep-files.rst @@ -0,0 +1,63 @@ +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. + +If the timesteps provided in the timesteps file are completely consumed (i.e. each timestep in the timesteps file has been taken +during evolution), but evolution is not complete (i.e. no stopping condition for evolution has been met), evolution is stopped and +a warning issued, and the status of the star/system being evolved is set to ``"User-provided timesteps exhausted"``. + +If the timesteps provided in the timesteps file are not completely consumed when evolution is complete (i.e. a stopping condition for +evolution has been met), a warning issued, and the status of the star/system evolved is set to ``"User-provided timesteps not consumed"``. + +Timesteps files can have blank lines and/or comment lines. 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). An excerpt from a timesteps file is shown below:: + + # timesteps files can have blank lines and/or comment lines + # comment lines begin with the '#' character + + # timesteps file for my experiment... + # description of file here + + 0.201713464000000 + 0.201935435000000 + + 0.020216053100000 + 0.020218312700000 + 0.020220604600000 + 0.002430922190000 + 0.002430922190000 + + # The following timestep is the NUCLEAR_MINIMUM_TIMESTEP (10^6 * TIMESTEP_QUANTUM = 10^6 * 10^-12 = 10^-6 Myr) + 0.000001000000000 + + 0.000171981412000 + 0.000168541784000 + 0.000165170948000 + 0.000161867529000 + 0.000099674841500 + 0.000097681344700 + 0.000095727717800 + 0.000093813163400 + + +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/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 4713e990a..87abf60d9 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) @@ -1879,7 +1879,7 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { bool accretorIsWD = m_Accretor->IsOneOf(WHITE_DWARFS); // Determine stability - bool isUnstable; + bool isUnstable = false; if (donorIsHeHGorHeGB && (caseBBAlwaysStable || caseBBAlwaysUnstable || (caseBBAlwaysUnstableOntoNSBH && accretorIsNSorBH))) { // Determine stability based on case BB isUnstable = (caseBBAlwaysUnstable || (caseBBAlwaysUnstableOntoNSBH && accretorIsNSorBH)); // Already established that donor is HeHG or HeGB - need to check if new case BB prescriptions are added } @@ -1917,32 +1917,43 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { double envMassDonor = m_Donor->Mass() - m_Donor->CoreMass(); bool isEnvelopeRemoved = false; - if (utils::Compare(m_Donor->CoreMass(), 0) > 0 && utils::Compare(envMassDonor, 0) > 0) { // donor has a core and an envelope - massDiffDonor = -envMassDonor; // set donor mass loss to (negative of) the envelope mass + if (utils::Compare(m_Donor->CoreMass(), 0) > 0 && utils::Compare(envMassDonor, 0) > 0) { // donor has a core and an envelope? + massDiffDonor = -envMassDonor; // yes - set donor mass loss to (negative of) the envelope mass isEnvelopeRemoved = true; } - else{ // donor has no envelope - massDiffDonor = -MassLossToFitInsideRocheLobe(this, m_Donor, m_Accretor, m_FractionAccreted); // use root solver to determine how much mass should be lost from the donor to allow it to fit within the Roche lobe - m_Donor->UpdateMinimumCoreMass(); // reset the minimum core mass following case A - } - double massGainAccretor = -massDiffDonor * m_FractionAccreted; // set accretor mass gain to mass loss * conservativeness + else { // donor has no envelope + massDiffDonor = MassLossToFitInsideRocheLobe(this, m_Donor, m_Accretor, m_FractionAccreted); // use root solver to determine how much mass should be lost from the donor to allow it to fit within the Roche lobe + if (massDiffDonor <= 0.0) { // no root found + // if donor cannot lose mass to fit inside Roche lobe, the only viable action is to enter CE phase + m_CEDetails.CEEnow = true; // flag CE + } + else { // have required mass loss + massDiffDonor = -massDiffDonor; // set mass difference + m_Donor->UpdateMinimumCoreMass(); // reset the minimum core mass following case A + } + } - m_Donor->SetMassTransferDiffAndResolveWDShellChange(massDiffDonor); // set new mass of donor - m_Accretor->SetMassTransferDiffAndResolveWDShellChange(massGainAccretor); // set new mass of accretor + if (!m_CEDetails.CEEnow) { // CE flagged? + // no + double massGainAccretor = -massDiffDonor * m_FractionAccreted; // set accretor mass gain to mass loss * conservativeness - aFinal = CalculateMassTransferOrbit(m_Donor->Mass(), massDiffDonor, *m_Accretor, m_FractionAccreted); // calculate new orbit - m_aMassTransferDiff = aFinal - aInitial; // set change in orbit (semi-major axis) + m_Donor->SetMassTransferDiffAndResolveWDShellChange(massDiffDonor); // set new mass of donor + m_Accretor->SetMassTransferDiffAndResolveWDShellChange(massGainAccretor); // set new mass of accretor + + aFinal = CalculateMassTransferOrbit(m_Donor->Mass(), massDiffDonor, *m_Accretor, m_FractionAccreted); // calculate new orbit + m_aMassTransferDiff = aFinal - aInitial; // set change in orbit (semi-major axis) - STELLAR_TYPE stellarTypeDonor = m_Donor->StellarType(); // donor stellar type before resolving envelope loss - if (isEnvelopeRemoved) m_Donor->ResolveEnvelopeLossAndSwitch(); // if this was an envelope stripping episode, resolve envelope loss - if (m_Donor->StellarType() != stellarTypeDonor) { // stellar type change? - (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::STELLAR_TYPE_CHANGE_DURING_MT); // yes - print (log) detailed output - } + STELLAR_TYPE stellarTypeDonor = m_Donor->StellarType(); // donor stellar type before resolving envelope loss + if (isEnvelopeRemoved) m_Donor->ResolveEnvelopeLossAndSwitch(); // if this was an envelope stripping episode, resolve envelope loss + if (m_Donor->StellarType() != stellarTypeDonor) { // stellar type change? + (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::STELLAR_TYPE_CHANGE_DURING_MT); // yes - print (log) detailed output + } - // Check if this was stable mass transfer after a CEE - if (m_CEDetails.CEEcount > 0 && !m_RLOFDetails.stableRLOFPostCEE) { - m_RLOFDetails.stableRLOFPostCEE = m_MassTransferTrackerHistory == MT_TRACKING::STABLE_2_TO_1_SURV || - m_MassTransferTrackerHistory == MT_TRACKING::STABLE_1_TO_2_SURV; + // Check if this was stable mass transfer after a CEE + if (m_CEDetails.CEEcount > 0 && !m_RLOFDetails.stableRLOFPostCEE) { + m_RLOFDetails.stableRLOFPostCEE = m_MassTransferTrackerHistory == MT_TRACKING::STABLE_2_TO_1_SURV || + m_MassTransferTrackerHistory == MT_TRACKING::STABLE_1_TO_2_SURV; + } } } @@ -2407,10 +2418,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 + } + + unsigned long 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 +2457,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 +2487,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 +2529,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 +2567,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 (evolutionStatus == EVOLUTION_STATUS::CONTINUE && 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..4f90e3f27 100644 --- a/src/BaseBinaryStar.h +++ b/src/BaseBinaryStar.h @@ -13,8 +13,8 @@ #include -#include // For float_distance. -#include // For boost::math::cbrt. +//#include // For float_distance. +//#include // For boost::math::cbrt. #include // for std::tuple and std::make_tuple. @@ -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; } @@ -519,37 +518,55 @@ class BaseBinaryStar { } - //Functor for the boost root finder to determine how much mass needs to be lost from a donor without an envelope in order to fit inside the Roche lobe + /* + * Functor for MassLossToFitInsideRocheLobe() + * + * + * Constructor: initialise the class + * template RadiusEqualsRocheLobeFunctor(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted) + * + * @param [IN] p_Binary (Pointer to) The binary star under examination + * @param [IN] p_Donor (Pointer to) The star donating mass + * @param [IN] p_Accretor (Pointer to) The star accreting mass + * @param [IN] p_FractionAccreted The fraction of the donated mass accreted by the accretor + * @param [IN] p_Error (Address of variable to record) Error encountered in functor + * + * Function: calculate radius difference after mass loss + * T RadiusEqualsRocheLobeFunctor(double const& p_dM) + * + * @param [IN] p_dM Mass to be donated + * @return Difference between star's Roche Lobe radius and radius after mass loss + */ template struct RadiusEqualsRocheLobeFunctor { - RadiusEqualsRocheLobeFunctor(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, ERROR *p_Error, double p_FractionAccreted) { + RadiusEqualsRocheLobeFunctor(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted, ERROR *p_Error) { m_Binary = p_Binary; m_Donor = p_Donor; m_Accretor = p_Accretor; m_Error = p_Error; m_FractionAccreted = p_FractionAccreted; } - T operator()(double const& dM) { + T operator()(double const& p_dM) { - if (dM >= m_Donor->Mass()) { // Can't remove more than the donor's mass - *m_Error = ERROR::TOO_MANY_RLOF_ITERATIONS; - return m_Donor->Radius(); + if (p_dM >= m_Donor->Mass()) { // Can't remove more than the donor's mass + *m_Error = ERROR::TOO_MANY_RLOF_ITERATIONS; // set error + return 1000.0 * ROOT_ABS_TOLERANCE; // arbitrary value to indicate no (sensible) solution found } double donorMass = m_Donor->Mass(); double accretorMass = m_Accretor->Mass(); BinaryConstituentStar* donorCopy = new BinaryConstituentStar(*m_Donor); - double semiMajorAxis = m_Binary->CalculateMassTransferOrbit(donorCopy->Mass(), -dM , *m_Accretor, m_FractionAccreted); - double RLRadius = semiMajorAxis * (1.0 - m_Binary->Eccentricity()) * CalculateRocheLobeRadius_Static(donorMass - dM, accretorMass + (m_Binary->FractionAccreted() * dM)) * AU_TO_RSOL; + double semiMajorAxis = m_Binary->CalculateMassTransferOrbit(donorCopy->Mass(), -p_dM , *m_Accretor, m_FractionAccreted); + double RLRadius = semiMajorAxis * (1.0 - m_Binary->Eccentricity()) * CalculateRocheLobeRadius_Static(donorMass - p_dM, accretorMass + (m_Binary->FractionAccreted() * p_dM)) * AU_TO_RSOL; - (void)donorCopy->UpdateAttributes(-dM, -dM * donorCopy->Mass0() / donorCopy->Mass()); + (void)donorCopy->UpdateAttributes(-p_dM, -p_dM * donorCopy->Mass0() / donorCopy->Mass()); // Modify donor Mass0 and Age for MS (including HeMS) and HG stars - donorCopy->UpdateInitialMass(); // update initial mass (MS, HG & HeMS) - donorCopy->UpdateAgeAfterMassLoss(); // update age (MS, HG & HeMS) + donorCopy->UpdateInitialMass(); // update initial mass (MS, HG & HeMS) + donorCopy->UpdateAgeAfterMassLoss(); // update age (MS, HG & HeMS) - (void)donorCopy->AgeOneTimestep(0.0); // recalculate radius of star - don't age - just update values + (void)donorCopy->AgeOneTimestep(0.0); // recalculate radius of star - don't age - just update values double thisRadiusAfterMassLoss = donorCopy->Radius(); @@ -564,45 +581,98 @@ class BaseBinaryStar { ERROR *m_Error; double m_FractionAccreted; }; - - - //Root solver to determine how much mass needs to be lost from a donor without an envelope in order to fit inside the Roche lobe - double MassLossToFitInsideRocheLobe(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted) { - using namespace std; // Help ADL of std functions. - using namespace boost::math::tools; // For bracket_and_solve_root. - - double guess = ADAPTIVE_RLOF_FRACTION_DONOR_GUESS * p_Donor->Mass(); // Rough guess at solution - double factor = ADAPTIVE_RLOF_SEARCH_FACTOR; // Size of search steps - - const boost::uintmax_t maxit = ADAPTIVE_RLOF_MAX_ITERATIONS; // Limit to maximum iterations. - boost::uintmax_t it = maxit; // Initially our chosen max iterations, but updated with actual. - bool is_rising = true; // So if result with guess is too low, then try increasing guess. - int digits = std::numeric_limits::digits; // Maximum possible binary digits accuracy for type T. - - // Some fraction of digits is used to control how accurate to try to make the result. - int get_digits = digits - 5; // We have to have a non-zero interval at each step, so - - // maximum accuracy is digits - 1. But we also have to - // allow for inaccuracy in f(x), otherwise the last few - // iterations just thrash around. - eps_tolerance tol(get_digits); // Set the tolerance. + + /* + * Root solver to determine how much mass needs to be lost from a donor without an envelope + * in order to fit inside the Roche lobe + * + * Uses boost::math::tools::bracket_and_solve_root() + * + * + * double MassLossToFitInsideRocheLobe(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted) + * + * @param [IN] p_Binary (Pointer to) The binary star under examination + * @param [IN] p_Donor (Pointer to) The star donating mass + * @param [IN] p_Accretor (Pointer to) The star accreting mass + * @param [IN] p_FractionAccreted The fraction of the donated mass accreted by the accretor + * @return Root found: will be -1.0 if no acceptable real root found + */ + double MassLossToFitInsideRocheLobe(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted) { - std::pair root(0.0, 0.0); - try { - ERROR error = ERROR::NONE; - root = bracket_and_solve_root(RadiusEqualsRocheLobeFunctor(p_Binary, p_Donor, p_Accretor, &error, p_FractionAccreted), guess, factor, is_rising, tol, it); - if (error != ERROR::NONE) SHOW_WARN(error); - } - catch(exception& e) { - SHOW_ERROR(ERROR::TOO_MANY_RLOF_ITERATIONS, e.what()); // Catch generic boost root finding error + const boost::uintmax_t maxit = ADAPTIVE_RLOF_MAX_ITERATIONS; // Limit to maximum iterations. + boost::uintmax_t it = maxit; // Initially our chosen max iterations, but updated with actual. + + // find root + // we use an iterative algorithm to find the root here: + // - if the root finder throws an exception, we stop and return a negative value for the root (indicating no root found) + // - if the root finder reaches the maximum number of (internal) iterations, we stop and return a negative value for the root (indicating no root found) + // - if the root finder returns a solution, we check that func(solution) = 0.0 +/ ROOT_ABS_TOLERANCE + // - if the solution is acceptable, we stop and return the solution + // - if the solution is not acceptable, we reduce the search step size and try again + // - if we reach the maximum number of search step reduction iterations, or the search step factor reduces to 1.0 (so search step size = 0.0), + // we stop and return a negative value for the root (indicating no root found) + + double guess = ADAPTIVE_RLOF_FRACTION_DONOR_GUESS * p_Donor->Mass(); // Rough guess at solution + + double factorFrac = ADAPTIVE_RLOF_SEARCH_FACTOR_FRAC; // search step size factor fractional part + double factor = 1.0 + factorFrac; // factor to determine search step size (size = guess * factor) + + std::pair root(-1.0, -1.0); // initialise root - default return + std::size_t tries = 0; // number of tries + bool done = false; // finished (found root or exceed maximum tries)? + ERROR error = ERROR::NONE; + RadiusEqualsRocheLobeFunctor func = RadiusEqualsRocheLobeFunctor(p_Binary, p_Donor, p_Accretor, p_FractionAccreted, &error); // no need to check error here + while (!done) { // while no error and acceptable root found + bool isRising = func((const double)guess) >= func((const double)guess * factor) ? false : true; // gradient direction from guess to upper search increment + + // run the root finder + // regardless of any exceptions or errors, display any problems as a warning, then + // check if the root returned is within tolerance - so even if the root finder + // bumped up against the maximum iterations, or couldn't bracket the root, use + // whatever value it ended with and check if it's good enough for us - not finding + // an acceptable root should be the exception rather than the rule, so this strategy + // shouldn't cause undue performance issues. + try { + error = ERROR::NONE; + root = boost::math::tools::bracket_and_solve_root(func, guess, factor, isRising, utils::BracketTolerance, it); // find root + // root finder returned without raising an exception + if (error != ERROR::NONE) { SHOW_WARN(error); } // root finder encountered an error + else if (it >= maxit) { SHOW_WARN(ERROR::TOO_MANY_RLOF_ITERATIONS); } // too many root finder iterations + } + catch(std::exception& e) { // catch generic boost root finding error + // root finder exception + // could be too many iterations, or unable to bracket root - it may not + // be a hard error - so no matter what the reason is that we are here, + // we'll just emit a warning and keep trying + if (it >= maxit) { SHOW_WARN(ERROR::TOO_MANY_RLOF_ITERATIONS); } // too many root finder iterations + else { SHOW_WARN(ERROR::ROOT_FINDER_FAILED, e.what()); } // some other problem - show it as a warning + } + + // we have a solution from the root finder - it may not be an acceptable solution + // so we check if it is within our preferred tolerance + if (fabs(func(root.first + (root.second - root.first) / 2.0)) <= ROOT_ABS_TOLERANCE) { // solution within tolerance? + done = true; // yes - we're done + } + else { // no - try again + // we don't have an acceptable solution - reduce search step size and try again + factorFrac /= 2.0; // reduce fractional part of factor + factor = 1.0 + factorFrac; // new search step size + tries++; // increment number of tries + if (tries > ADAPTIVE_RLOF_MAX_TRIES || fabs(factor - 1.0) <= ROOT_ABS_TOLERANCE) { // too many tries, or step size 0.0? + // we've tried as much as we can - fail here with -ve return value + root.first = -1.0; // yes - set error return + root.second = -1.0; + SHOW_WARN(ERROR::TOO_MANY_RLOF_TRIES); // show warning + done = true; // we're done + } + } } - SHOW_WARN_IF(it >= maxit, ERROR::TOO_MANY_RLOF_ITERATIONS); - return root.first + (root.second - root.first) / 2.0; // Midway between brackets is our result, if necessary we could return the result as an interval here. + return root.first + (root.second - root.first) / 2.0; // Midway between brackets is our result, if necessary we could return the result as an interval here. } - + /* * Root solver to determine rotational frequency after synchronisation for tides * @@ -617,19 +687,13 @@ class BaseBinaryStar { * @param [IN] p_I2 Moment of inertia of star 1 * @param [IN] p_Ltot Total angular momentum for binary * @param [IN] p_Guess Initial guess for value of root - * @return Root found: will be -1.0 if no real root found + * @return Root found: will be -1.0 if no acceptable real root found */ double OmegaAfterSynchronisation(const double p_M1, const double p_M2, const double p_I1, const double p_I2, const double p_Ltot, const double p_Guess) { - const boost::uintmax_t maxit = TIDES_OMEGA_MAX_ITERATIONS; // maximum iterations + const boost::uintmax_t maxit = TIDES_OMEGA_MAX_ITERATIONS; // maximum iterations for root finder boost::uintmax_t it = maxit; // initially max iterations, but updated with actual count - int digits = std::numeric_limits::digits; // maximum possible binary digits accuracy - int get_digits = digits - 5; // we have to have a non-zero interval at each step - - // maximum accuracy is digits - 1. But we also have to allow for inaccuracy - // in the functor, otherwise the last few iterations just thrash around. - boost::math::tools::eps_tolerance tol(get_digits); // tolerance - + // define functor // function: (I_1 + I_2) Omega + L(Omega) - p_Ltot = 0 // where L(Omega) = b*Omega(-1/3) @@ -640,25 +704,63 @@ class BaseBinaryStar { auto func = [a, b, c](double x) -> double { return (a * x) + (b / std::cbrt(x)) + c; }; // functor // find root - double factor = TIDES_OMEGA_SEARCH_FACTOR; // size of search steps - bool is_rising = func(p_Guess) > func(p_Guess * factor) ? false : true; // so bracket_and_solve_root() knows whether to increase or decrease guess per iteration - - std::pair root(-1.0, -1.0); // initialise root - try { - root = boost::math::tools::bracket_and_solve_root(func, p_Guess, factor, is_rising, tol, it); // iterate to find root - } - catch(std::exception& e) { // catch generic boost root finding error - root.first = -1.0; // set error return - root.second = -1.0; - if (it < maxit) { // not too many iterations? - SHOW_ERROR(ERROR::ROOT_FINDER_FAILED, e.what()); // no - some other error - show it + // we use an iterative algorithm to find the root here: + // - if the root finder throws an exception, we stop and return a negative value for the root (indicating no root found) + // - if the root finder reaches the maximum number of (internal) iterations, we stop and return a negative value for the root (indicating no root found) + // - if the root finder returns a solution, we check that func(solution) = 0.0 +/ ROOT_ABS_TOLERANCE + // - if the solution is acceptable, we stop and return the solution + // - if the solution is not acceptable, we reduce the search step size and try again + // - if we reach the maximum number of search step reduction iterations, or the search step factor reduces to 1.0 (so search step size = 0.0), + // we stop and return a negative value for the root (indicating no root found) + + double factorFrac = TIDES_OMEGA_SEARCH_FACTOR_FRAC; // search step size factor fractional part + double factor = 1.0 + factorFrac; // factor to determine search step size (size = guess * factor) + + std::pair root(-1.0, -1.0); // initialise root - default return + std::size_t tries = 0; // number of tries + bool done = false; // finished (found root or exceed maximum tries)? + while (!done) { // while no acceptable root found + bool isRising = func(p_Guess) >= func(p_Guess * factor) ? false : true; // gradient direction from guess to upper search increment + + // run the root finder + // regardless of any exceptions or errors, display any problems as a warning, then + // check if the root returned is within tolerance - so even if the root finder + // bumped up against the maximum iterations, or couldn't bracket the root, use + // whatever value it ended with and check if it's good enough for us - not finding + // an acceptable root should be the exception rather than the rule, so this strategy + // shouldn't cause undue performance issues. + try { + root = boost::math::tools::bracket_and_solve_root(func, p_Guess, factor, isRising, utils::BracketTolerance, it); // find root + // root finder returned without raising an exception + if (it >= maxit) { SHOW_WARN(ERROR::TOO_MANY_OMEGA_ITERATIONS); } // too many root finder iterations + } + catch(std::exception& e) { // catch generic boost root finding error + // root finder exception + // could be too many iterations, or unable to bracket root - it may not + // be a hard error - so no matter what the reason is that we are here, + // we'll just emit a warning and keep trying + if (it >= maxit) { SHOW_WARN(ERROR::TOO_MANY_OMEGA_ITERATIONS); } // too many root finder iterations + else { SHOW_WARN(ERROR::ROOT_FINDER_FAILED, e.what()); } // some other problem - show it as a warning } - } - if (it >= maxit) { // too many iterations? - root.first = -1.0; // yes - set error return - root.second = -1.0; - SHOW_WARN(ERROR::TOO_MANY_OMEGA_ITERATIONS); // show warning + // we have a solution from the root finder - it may not be an acceptable solution + // so we check if it is within our preferred tolerance + if (fabs(func(root.first + (root.second - root.first) / 2.0)) <= ROOT_ABS_TOLERANCE) { // solution within tolerance? + done = true; // yes - we're done + } + else { // no - try again + // we don't have an acceptable solution - reduce search step size and try again + factorFrac /= 2.0; // reduce fractional part of factor + factor = 1.0 + factorFrac; // new search step size + tries++; // increment number of tries + if (tries > TIDES_OMEGA_MAX_TRIES || fabs(factor - 1.0) <= ROOT_ABS_TOLERANCE) { // too many tries, or step size 0.0? + // we've tried as much as we can - fail here with -ve return value + root.first = -1.0; // yes - set error return + root.second = -1.0; + SHOW_WARN(ERROR::TOO_MANY_OMEGA_TRIES); // show warning + done = true; // we're done + } + } } return root.first + (root.second - root.first) / 2.0; // midway between brackets (could return brackets...) diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index 21308dc3c..c1b8fe2ac 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; @@ -377,6 +377,7 @@ COMPAS_VARIABLE BaseStar::StellarPropertyValue(const T_ANY_PROPERTY p_Property) case ANY_STAR_PROPERTY::TIMESCALE_MS: value = Timescale(TIMESCALE::tMS); break; case ANY_STAR_PROPERTY::TOTAL_MASS_AT_COMPACT_OBJECT_FORMATION: value = SN_TotalMassAtCOFormation(); break; case ANY_STAR_PROPERTY::TRUE_ANOMALY: value = SN_TrueAnomaly(); break; + case ANY_STAR_PROPERTY::TZAMS: value = TZAMS()*TSOL; break; case ANY_STAR_PROPERTY::ZETA_HURLEY: value = CalculateZetaAdiabaticHurley2002(m_CoreMass); break; case ANY_STAR_PROPERTY::ZETA_HURLEY_HE: value = CalculateZetaAdiabaticHurley2002(m_HeCoreMass); break; case ANY_STAR_PROPERTY::ZETA_SOBERMAN: value = CalculateZetaAdiabaticSPH(m_CoreMass); break; @@ -2587,7 +2588,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 +2606,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: 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; } @@ -2621,18 +2622,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; @@ -2684,11 +2683,11 @@ double BaseStar::CalculateMassLossRate() { default: // unknown mass-loss prescription SHOW_WARN(ERROR::UNKNOWN_MASS_LOSS_PRESCRIPTION, "Using HURLEY"); // 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 +2697,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 +2726,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 +2740,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 +2755,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,21 +2778,28 @@ 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); - - UpdateInitialMass(); // update effective initial mass (MS, HG & HeMS) - UpdateAgeAfterMassLoss(); // update age (MS, HG & HeMS) - ApplyMassTransferRejuvenationFactor(); // apply age rejuvenation factor + double mass = CalculateMassLossValues(true, p_UpdateMDt); // calculate new values assuming mass loss applied + + // JR: this is here 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()) + // Don't resolve envelope loss here (JR: we're not going to switch anyway... need to revisit this) + STELLAR_TYPE st = UpdateAttributesAndAgeOneTimestep(mass - m_Mass, 0.0, 0.0, false, false); // recalculate stellar attributes + if (st != m_StellarType) { // should switch? + SHOW_WARN(ERROR::SWITCH_NOT_TAKEN); // show warning if we think we should switch again... + } + + UpdateInitialMass(); // update effective initial mass (MS, HG & HeMS) + UpdateAgeAfterMassLoss(); // update age (MS, HG & HeMS) + ApplyMassTransferRejuvenationFactor(); // apply age rejuvenation factor } } @@ -3843,6 +3853,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 +3880,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 +3956,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) @@ -3970,46 +3982,49 @@ void BaseStar::UpdateAttributesAndAgeOneTimestepPreamble(const double p_DeltaMas * @param [IN] p_DeltaTime The timestep to evolve in Myr * @param [IN] p_ForceRecalculate Specifies whether the star's attributes should be recalculated even if the three deltas are 0.0 * (optional, default = false) + * @param [IN] p_ResolveEnvelopeLoss Specifies whether envelope loss should be resolved here + * (optional, default = true) + * JR: this is a bit of a kludge to resolve problems introduced by modifying stellar attributes in + * anticipation of switching stellar type, but using those attributes before the actual switch + * to the new stellar type - we need to resolve those situations in the code. + * The bottom line is that this parameter is a temporary fix while I develop the permannet fix. * @return Stellar type to which star should evolve */ STELLAR_TYPE BaseStar::UpdateAttributesAndAgeOneTimestep(const double p_DeltaMass, const double p_DeltaMass0, const double p_DeltaTime, - const bool p_ForceRecalculate) { - STELLAR_TYPE stellarType = m_StellarType; // default is no change + const bool p_ForceRecalculate, + const bool p_ResolveEnvelopeLoss) { + STELLAR_TYPE stellarType = m_StellarType; // default is no change - if (ShouldBeMasslessRemnant()) { // Do not update the star if it lost all of its mass + if (ShouldBeMasslessRemnant()) { // Do not update the star if it lost all of its mass stellarType = STELLAR_TYPE::MASSLESS_REMNANT; } else { - stellarType = ResolveSupernova(); // handle supernova - - if (stellarType == m_StellarType) { // still on phase? + 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 - - 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? - // yes - update attributes - AgeOneTimestepPreamble(p_DeltaTime); // advance dt, age, simulation time - - if (ShouldSkipPhase()) { // skip phase? - stellarType = ResolveSkippedPhase(); // phase skipped - do what's 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 - } + 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... + p_DeltaTime > 0) { // age/time advance? (don't use utils::Compare() here) + // yes - update attributes + 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 + stellarType = EvolveOnPhase(p_DeltaTime); // evolve on phase + if (stellarType == m_StellarType) { // need to switch to new stellar type? + stellarType = ResolveEndOfPhase(p_ResolveEnvelopeLoss); // no - check for need to move off phase + } } } } } - return stellarType; // stellar type to which star should evolve + return stellarType; // stellar type to which star should evolve } @@ -4028,7 +4043,6 @@ void BaseStar::AgeOneTimestepPreamble(const double p_DeltaTime) { m_Age += p_DeltaTime; // advance age of star m_Dt = p_DeltaTime; // set timestep } - EvolveOneTimestepPreamble(); } @@ -4081,16 +4095,16 @@ void BaseStar::UpdateMassTransferDonorHistory() { * Evolve the star on it's current phase - take one timestep on the current phase * * - * STELLAR_TYPE EvolveOnPhase() + * STELLAR_TYPE EvolveOnPhase(const double p_DeltaTime) * + * @param [IN] p_DeltaTime Timestep in Myr * @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(); @@ -4107,9 +4121,11 @@ STELLAR_TYPE BaseStar::EvolveOnPhase() { m_Temperature = CalculateTemperatureOnPhase(); - 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; + 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; + } } } @@ -4121,17 +4137,23 @@ STELLAR_TYPE BaseStar::EvolveOnPhase() { * Evolve the star onto the next phase if necessary - take one timestep at the end of the current phase * * - * STELLAR_TYPE ResolveEndOfPhase() + * STELLAR_TYPE ResolveEndOfPhase(const bool p_ResolveEnvelopeLoss) * + * @param [IN] p_ResolveEnvelopeLoss Specifies whether envelope loss should be resolved here + * (optional, default = true) + * JR: this is a bit of a kludge to resolve problems introduced by modifying stellar attributes in + * anticipation of switching stellar type, but using those attributes before the actual switch + * to the new stellar type - we need to resolve those situations in the code. * @return Stellar Type to which star should evolve - unchanged if not moving off current phase */ -STELLAR_TYPE BaseStar::ResolveEndOfPhase() { +STELLAR_TYPE BaseStar::ResolveEndOfPhase(const bool p_ResolveEnvelopeLoss) { STELLAR_TYPE stellarType = m_StellarType; if (IsEndOfPhase()) { // end of phase - stellarType = ResolveEnvelopeLoss(); // resolve envelope loss if it occurs + if (p_ResolveEnvelopeLoss) stellarType = ResolveEnvelopeLoss(); // if required, resolve envelope loss if it occurs + if (stellarType == m_StellarType) { // staying on phase? m_Tau = CalculateTauAtPhaseEnd(); diff --git a/src/BaseStar.h b/src/BaseStar.h index 665462be8..b5a5c7f1c 100644 --- a/src/BaseStar.h +++ b/src/BaseStar.h @@ -137,6 +137,7 @@ class BaseStar { double Temperature() const { return m_Temperature; } double Time() const { return m_Time; } double Timescale(TIMESCALE p_Timescale) const { return m_Timescales[static_cast(p_Timescale)]; } + double TZAMS() const { return m_TZAMS; } virtual ACCRETION_REGIME WhiteDwarfAccretionRegime() const { return ACCRETION_REGIME::NONE; } double XExponent() const { return m_XExponent; } @@ -236,7 +237,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; } @@ -250,7 +251,8 @@ class BaseStar { STELLAR_TYPE UpdateAttributesAndAgeOneTimestep(const double p_DeltaMass, const double p_DeltaMass0, const double p_DeltaTime, - const bool p_ForceRecalculate); + const bool p_ForceRecalculate = false, + const bool p_ResolveEnvelopeLoss = true); virtual void UpdateInitialMass() { } // Default is NO-OP @@ -264,8 +266,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 } @@ -583,7 +584,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; } @@ -615,14 +616,14 @@ class BaseStar { virtual void PerturbLuminosityAndRadius() { } // NO-OP virtual void PerturbLuminosityAndRadiusAtPhaseEnd() { PerturbLuminosityAndRadiusOnPhase(); } // Same as on phase virtual void PerturbLuminosityAndRadiusOnPhase() { PerturbLuminosityAndRadius(); } - STELLAR_TYPE ResolveEndOfPhase(); + STELLAR_TYPE ResolveEndOfPhase(const bool p_ResolveEnvelopeLoss = true); virtual void ResolveHeliumFlash() { } virtual STELLAR_TYPE ResolveSkippedPhase() { return EvolveToNextPhase(); } // Default is evolve to next phase virtual STELLAR_TYPE ResolveSupernova() { return m_StellarType; } // Default is NO-OP 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..0d8d0c227 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); // He mass capped at total mass (should become HeMS star) } @@ -1358,6 +1358,7 @@ double CHeB::CalculateLifetimeOnBluePhase(const double p_Mass) { #undef b } + /* * Determine whether star should continue to evolve on phase * @@ -1366,15 +1367,15 @@ double CHeB::CalculateLifetimeOnBluePhase(const double p_Mass) { * * @return true if evolution should continue on phase, false otherwise */ - bool CHeB::ShouldEvolveOnPhase() const { - bool afterHeIgnition = (m_Age >= m_Timescales[static_cast(TIMESCALE::tHeI)]); + bool afterHeIgnition = (m_Age >= m_Timescales[static_cast(TIMESCALE::tHeI)]); bool beforeEndOfHeBurning = (m_Age < (m_Timescales[static_cast(TIMESCALE::tHeI)] + m_Timescales[static_cast(TIMESCALE::tHe)])); - bool coreIsNotTooMassive = (m_HeCoreMass < m_Mass); + bool coreIsNotTooMassive = (m_HeCoreMass < m_Mass); // Evolve on CHeB phase if age after He Ign and while He Burning and He core mass does not exceed total mass (could happen due to mass loss) return (afterHeIgnition && beforeEndOfHeBurning && coreIsNotTooMassive && !ShouldEnvelopeBeExpelledByPulsations()); } + /////////////////////////////////////////////////////////////////////////////////////// // // // MISCELLANEOUS FUNCTIONS / CONTROL FUNCTIONS // @@ -1432,11 +1433,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 +1467,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 +1485,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..b09a95beb 100755 --- a/src/EAGB.cpp +++ b/src/EAGB.cpp @@ -772,6 +772,12 @@ double EAGB::CalculateRadiusOnPhase_Static(const double p_Mass, const DBL_VECTOR &p_BnCoefficients) { #define b p_BnCoefficients // for convenience and readability - undefined at end of function + // sanity check for mass and luminosity - just return 0.0 if mass or luminosity <= 0 + // doing this will save some compute cycles - but it is not strictly what the equation in Hurley says + // (Hurley et al. 2000, eq 74 and immediately following) - there if mass is 0.0 but luminosity is + // non-zero we get a non-zero value for radius (shouldn't happen, but we need to code for all possibilities) + if (utils::Compare(p_Mass, 0.0) <= 0 || utils::Compare(p_Luminosity, 0.0) <= 0) return 0.0; + // calculate radius constant A (Hurley et al. 2000, eq 74) // and coefficient b(50) double A; @@ -1000,17 +1006,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 +1026,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 +1038,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..936d76f4d 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; 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..d631284f2 100755 --- a/src/HG.h +++ b/src/HG.h @@ -5,6 +5,7 @@ #include "typedefs.h" #include "profiling.h" #include "utils.h" + #include #include "GiantBranch.h" @@ -39,21 +40,27 @@ class HG: virtual public BaseStar, public GiantBranch { CalculateTimescales(); // Initialise timescales m_Age = m_Timescales[static_cast(TIMESCALE::tMS)]; // Set age appropriately - //Update stellar properties at start of HG phase (since core definition changes) + // update stellar properties at start of HG phase (since core definition changes) CalculateGBParams(); - //update effective "initial" mass m_Mass0 so that the core mass is at least equal to the minimum core mass (only relevant if RetainCoreMassDuringCaseAMassTransfer() ) but no more than total mass + // update effective "initial" mass (m_Mass0) so that the core mass is at least equal to the minimum core mass but no more than total mass + // (only relevant if RetainCoreMassDuringCaseAMassTransfer()) if(utils::Compare(CalculateCoreMassOnPhase(m_Mass0, m_Age), std::min(m_Mass, MinimumCoreMass())) < 0) { - m_Mass0 = Mass0ToMatchDesiredCoreMass(this, std::min(m_Mass,MinimumCoreMass())); + double desiredCoreMass = std::min(m_Mass, MinimumCoreMass()); // desired core mass + m_Mass0 = Mass0ToMatchDesiredCoreMass(this, desiredCoreMass); // use root finder to find new core mass estimate + if (m_Mass0 <= 0.0) { // no root found - no solution for estimated core mass + // if no root found we keep m_Mass0 equal to the total mass + m_Mass0 = m_Mass; + } CalculateTimescales(); 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 } @@ -118,66 +125,123 @@ class HG: virtual public BaseStar, public GiantBranch { void UpdateInitialMass(); // Per Hurley et al. 2000, section 7.1 - - //Functor for the boost root finder to determine the "initial mass" m_Mass0 based on desired core mass + + /* + * Functor for Mass0ToMatchDesiredCoreMass() + * + * + * Constructor: initialise the class + * template Mass0YieldsDesiredCoreMassFunctor(HG *p_Star, double p_DesiredCoreMass) + * + * @param [IN] p_Star (Pointer to) The star under examination + * @param [IN] p_DesiredCoreMass The desired core mass + * + * Function: core mass difference after mass loss + * T RadiusEqualsRocheLobeFunctor(double const& p_dM) + * + * @param [IN] p_GuessMass0 Guess for Mass0 + * @return Difference between estimated core mass (based on p_GuessMass0) and desired core mass (p_DesiredCoreMass) + */ template - struct Mass0YieldsDesiredCoreMassFunctor - { - Mass0YieldsDesiredCoreMassFunctor(HG *p_Star, double p_DesiredCoreMass, ERROR *p_Error) - { + struct Mass0YieldsDesiredCoreMassFunctor { + Mass0YieldsDesiredCoreMassFunctor(HG *p_Star, double p_DesiredCoreMass) { m_Star = p_Star; m_DesiredCoreMass = p_DesiredCoreMass; - m_Error = p_Error; } - T operator()(double const& guessMass0) - { + T operator()(double const& p_GuessMass0) { HG *copy = new HG(*m_Star, false); - copy->UpdateAttributesAndAgeOneTimestep(0.0, guessMass0 - copy->Mass0(), 0.0, true); - double coreMassEstimate = copy->CalculateCoreMassOnPhase(guessMass0, copy->Age()); + copy->UpdateAttributesAndAgeOneTimestep(0.0, p_GuessMass0 - copy->Mass0(), 0.0, true); + double coreMassEstimate = copy->CalculateCoreMassOnPhase(p_GuessMass0, copy->Age()); delete copy; copy = nullptr; return (coreMassEstimate - m_DesiredCoreMass); } private: HG *m_Star; double m_DesiredCoreMass; - ERROR *m_Error; }; - //Root solver to determine "initial mass" m_Mass0 based on desired core mass - double Mass0ToMatchDesiredCoreMass(HG * p_Star, double p_DesiredCoreMass) - { - using namespace std; // Help ADL of std functions. - using namespace boost::math::tools; // For bracket_and_solve_root. - - double guess = p_Star->Mass(); // Rough guess at solution - double factor = ADAPTIVE_MASS0_SEARCH_FACTOR; // Size of search steps - - const boost::uintmax_t maxit = ADAPTIVE_MASS0_MAX_ITERATIONS; // Limit to maximum iterations. - boost::uintmax_t it = maxit; // Initially our chosen max iterations, but updated with actual. - bool is_rising = true; // So if result with guess is too low, then try increasing guess. - int digits = std::numeric_limits::digits; // Maximum possible binary digits accuracy for type T. + /* + * Root solver to determine "initial" mass (m_Mass0) based on desired core mass + * + * Uses boost::math::tools::bracket_and_solve_root() + * + * + * double Mass0ToMatchDesiredCoreMass(HG * p_Star, double p_DesiredCoreMass) + * + * @param [IN] p_Star (Pointer to) The star under examination + * @param [IN] p_DesiredCoreMass The desired core mass + * @return Root found: will be -1.0 if no acceptable real root found + */ + double Mass0ToMatchDesiredCoreMass(HG * p_Star, double p_DesiredCoreMass) { + + const boost::uintmax_t maxit = ADAPTIVE_MASS0_MAX_ITERATIONS; // Limit to maximum iterations. + boost::uintmax_t it = maxit; // Initially our chosen max iterations, but updated with actual. + + // find root + // we use an iterative algorithm to find the root here: + // - if the root finder throws an exception, we stop and return a negative value for the root (indicating no root found) + // - if the root finder reaches the maximum number of (internal) iterations, we stop and return a negative value for the root (indicating no root found) + // - if the root finder returns a solution, we check that func(solution) = 0.0 +/ ROOT_ABS_TOLERANCE + // - if the solution is acceptable, we stop and return the solution + // - if the solution is not acceptable, we reduce the search step size and try again + // - if we reach the maximum number of search step reduction iterations, or the search step factor reduces to 1.0 (so search step size = 0.0), + // we stop and return a negative value for the root (indicating no root found) + + double guess = p_Star->Mass(); // Rough guess at solution - // Some fraction of digits is used to control how accurate to try to make the result. - int get_digits = digits - 5; // We have to have a non-zero interval at each step, so - - // maximum accuracy is digits - 1. But we also have to - // allow for inaccuracy in f(x), otherwise the last few - // iterations just thrash around. - eps_tolerance tol(get_digits); // Set the tolerance. - - std::pair root; - try { - ERROR error = ERROR::NONE; - root = bracket_and_solve_root(Mass0YieldsDesiredCoreMassFunctor(p_Star, p_DesiredCoreMass, &error), guess, factor, is_rising, tol, it); - if (error != ERROR::NONE) SHOW_WARN(error); - } - catch(exception& e) { - SHOW_ERROR(ERROR::TOO_MANY_MASS0_ITERATIONS, e.what()); // Catch generic boost root finding error + double factorFrac = ADAPTIVE_MASS0_SEARCH_FACTOR_FRAC; // search step size factor fractional part + double factor = 1.0 + factorFrac; // factor to determine search step size (size = guess * factor) + + std::pair root(p_DesiredCoreMass, 0.0); // initialise root - default return + std::size_t tries = 0; // number of tries + bool done = false; // finished (found root or exceed maximum tries)? + Mass0YieldsDesiredCoreMassFunctor func = Mass0YieldsDesiredCoreMassFunctor(p_Star, p_DesiredCoreMass); + while (!done) { // while no acceptable root found + bool isRising = func((const double)guess) >= func((const double)guess * factor) ? false : true; // gradient direction from guess to upper search increment + + // run the root finder + // regardless of any exceptions or errors, display any problems as a warning, then + // check if the root returned is within tolerance - so even if the root finder + // bumped up against the maximum iterations, or couldn't bracket the root, use + // whatever value it ended with and check if it's good enough for us - not finding + // an acceptable root should be the exception rather than the rule, so this strategy + // shouldn't cause undue performance issues. + try { + root = boost::math::tools::bracket_and_solve_root(func, guess, factor, isRising, utils::BracketTolerance, it); // find root + // root finder returned without raising an exception + if (it >= maxit) { SHOW_WARN(ERROR::TOO_MANY_MASS0_ITERATIONS); } // too many root finder iterations + } + catch(std::exception& e) { // catch generic boost root finding error + // root finder exception + // could be too many iterations, or unable to bracket root - it may not + // be a hard error - so no matter what the reason is that we are here, + // we'll just emit a warning and keep trying + if (it >= maxit) { SHOW_WARN(ERROR::TOO_MANY_MASS0_ITERATIONS); } // too many root finder iterations + else { SHOW_WARN(ERROR::ROOT_FINDER_FAILED, e.what()); } // some other problem - show it as a warning + } + + // we have a solution from the root finder - it may not be an acceptable solution + // so we check if it is within our preferred tolerance + if (fabs(func(root.first + (root.second - root.first) / 2.0)) <= ROOT_ABS_TOLERANCE) { // solution within tolerance? + done = true; // yes - we're done + } + else { // no - try again + // we don't have an acceptable solution - reduce search step size and try again + factorFrac /= 2.0; // reduce fractional part of factor + factor = 1.0 + factorFrac; // new search step size + tries++; // increment number of tries + if (tries > ADAPTIVE_MASS0_MAX_TRIES || fabs(factor - 1.0) <= ROOT_ABS_TOLERANCE) { // too many tries, or step size 0.0? + // we've tried as much as we can - fail here with -ve return value + root.first = -1.0; // yes - set error return + root.second = -1.0; + SHOW_WARN(ERROR::TOO_MANY_MASS0_TRIES); // show warning + done = true; // we're done + } + } } - SHOW_WARN_IF(it>=maxit, ERROR::TOO_MANY_MASS0_ITERATIONS); - return root.first + (root.second - root.first) / 2.0; // Midway between brackets is our result, if necessary we could return the result as an interval here. + return root.first + (root.second - root.first) / 2.0; // Midway between brackets is our result, if necessary we could return the result as an interval here. } }; diff --git a/src/HeGB.cpp b/src/HeGB.cpp index e4ca36bef..1a38953b4 100755 --- a/src/HeGB.cpp +++ b/src/HeGB.cpp @@ -41,6 +41,9 @@ double HeGB::CalculateLuminosityOnPhase_Static(const double p_CoreMass, const do */ std::tuple HeGB::CalculateRadiusOnPhase_Static(const double p_Mass, const double p_Luminosity) { + // sanity check for mass and luminosity - just return 0.0 if mass or luminosity <= 0 + if (utils::Compare(p_Mass, 0.0) <= 0 || utils::Compare(p_Luminosity, 0.0) <= 0) return std::make_tuple(0.0, 0.0); + double RZHe = HeMS::CalculateRadiusAtZAMS_Static(p_Mass); double LTHe = CalculateLuminosityAtPhaseEnd_Static(p_Mass); 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/HeHG.h b/src/HeHG.h index a198e8983..cc711e14f 100755 --- a/src/HeHG.h +++ b/src/HeHG.h @@ -77,6 +77,8 @@ class HeHG: virtual public BaseStar, public HeMS { double CalculateMassTransferRejuvenationFactor() const; + double CalculateMomentOfInertia() const { return GiantBranch::CalculateMomentOfInertia(); } + double CalculatePerturbationMu() const; double CalculatePerturbationMuAtPhaseEnd() const { return m_Mu; } // NO-OP diff --git a/src/HeMS.cpp b/src/HeMS.cpp index 8febed40a..da53c2c01 100755 --- a/src/HeMS.cpp +++ b/src/HeMS.cpp @@ -144,6 +144,9 @@ double HeMS::CalculateRadiusAtZAMS_Static(const double p_Mass) { */ double HeMS::CalculateRadiusOnPhase_Static(const double p_Mass, const double p_Tau) { + // sanity check for mass - just return 0.0 if mass <= 0 + if (utils::Compare(p_Mass, 0.0) <= 0) return 0.0; + double tau_6 = p_Tau * p_Tau * p_Tau * p_Tau * p_Tau * p_Tau; // pow() is slow - use multiplication double beta = std::max(0.0, 0.4 - 0.22 * log10(p_Mass)); @@ -464,8 +467,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/HeMS.h b/src/HeMS.h index 6fb67fb4f..806fdeb20 100644 --- a/src/HeMS.h +++ b/src/HeMS.h @@ -84,6 +84,8 @@ class HeMS: virtual public BaseStar, public TPAGB { double CalculateMassTransferRejuvenationFactor() const; + double CalculateMomentOfInertia() const { return MainSequence::CalculateMomentOfInertia(); } + double CalculatePerturbationMu() const { return 5.0; } // Hurley et al. 2000, eqs 97 & 98 double CalculateRadialExtentConvectiveEnvelope() const { return BaseStar::CalculateRadialExtentConvectiveEnvelope(); } // HeMS stars don't have a convective envelope 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 7d0c78ff4..9b7d353d4 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; @@ -972,7 +974,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() ) @@ -1811,6 +1813,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 a5063ae70..fc0ba14f5 100755 --- a/src/Options.h +++ b/src/Options.h @@ -532,6 +532,8 @@ class Options { "store-input-files", "switch-log", + "timesteps-filename", + "use-mass-loss", "VMS-mass-loss", @@ -611,6 +613,8 @@ class Options { "store-input-files", "switch-log", + "timesteps-filename", + "version", "v", "yaml-template" @@ -682,12 +686,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 + unsigned long 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) @@ -1374,7 +1380,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); } + unsigned long 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); } @@ -1480,6 +1486,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..4bf414010 100644 --- a/src/Star.cpp +++ b/src/Star.cpp @@ -297,7 +297,10 @@ STELLAR_TYPE Star::UpdateAttributesAndAgeOneTimestep(const double p_DeltaMass, // (this could get recursive, but shouldn't...) if (stellarTypePrev == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS && stellarType == STELLAR_TYPE::NAKED_HELIUM_STAR_MS) { // discontinuous transition? - UpdateAttributes(0.0, 0.0, true); // yes - recalculate stellar attributes + if (UpdateAttributes(0.0, 0.0, true) != stellarType) { // yes - recalculate stellar attributes + // JR: need to revisit this - should we (queue) switch here? + SHOW_WARN(ERROR::SWITCH_NOT_TAKEN); // show warning if we think we should switch again... + } } } @@ -330,8 +333,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) <--- 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 @@ -397,12 +400,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; @@ -425,7 +435,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 +462,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 + + (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 +489,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; + unsigned long int stepNum = 0; // initialise step number 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 BSE 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/WhiteDwarfs.cpp b/src/WhiteDwarfs.cpp index 28c33bf58..a2b1790ce 100644 --- a/src/WhiteDwarfs.cpp +++ b/src/WhiteDwarfs.cpp @@ -151,8 +151,13 @@ double WhiteDwarfs::CalculateLuminosityOnPhase_Static(const double p_Mass, const * @return Radius of a White Dwarf in Rsol (since WD is ~ Earth sized, expect answer around 0.009) */ double WhiteDwarfs::CalculateRadiusOnPhase_Static(const double p_Mass) { + + // sanity check for mass - just return 0.0 if mass <= 0 + if (utils::Compare(p_Mass, 0.0) <= 0) return 0.0; + double MCH_Mass_one_third = std::cbrt(MCH / p_Mass); double MCH_Mass_two_thirds = MCH_Mass_one_third * MCH_Mass_one_third; + return std::max(NEUTRON_STAR_RADIUS, 0.0115 * std::sqrt((MCH_Mass_two_thirds - 1.0 / MCH_Mass_two_thirds))); } @@ -210,6 +215,12 @@ STELLAR_TYPE WhiteDwarfs::ResolveAIC() { if (!IsSupernova()) return m_StellarType; // shouldn't be here if no SN + m_SupernovaDetails.totalMassAtCOFormation = m_Mass; + m_SupernovaDetails.HeCoreMassAtCOFormation = m_HeCoreMass; + m_SupernovaDetails.COCoreMassAtCOFormation = m_COCoreMass; + m_SupernovaDetails.coreMassAtCOFormation = m_CoreMass; + SetSNHydrogenContent(); // SN to be H-poor. + m_Mass = MECS_REM; // defined in constants.h m_SupernovaDetails.drawnKickMagnitude = 0.0; @@ -235,11 +246,17 @@ STELLAR_TYPE WhiteDwarfs::ResolveSNIa() { if (!IsSupernova()) return m_StellarType; // shouldn't be here if no SN + m_SupernovaDetails.totalMassAtCOFormation = m_Mass; + m_SupernovaDetails.HeCoreMassAtCOFormation = m_HeCoreMass; + m_SupernovaDetails.COCoreMassAtCOFormation = m_COCoreMass; + m_SupernovaDetails.coreMassAtCOFormation = m_CoreMass; + SetSNHydrogenContent(); // SN to be H-poor. + m_Mass = 0.0; m_Radius = 0.0; m_Luminosity = 0.0; m_Age = 0.0; - + m_SupernovaDetails.drawnKickMagnitude = 0.0; m_SupernovaDetails.kickMagnitude = 0.0; diff --git a/src/changelog.h b/src/changelog.h index 378633373..262acc8bf 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1090,8 +1090,30 @@ // - Fix for issue #1048 // 02.41.05 YS - Jan 31, 2024 - Bug fix: // - Fix for issue #1058: fixing calculation of pulsar spin period +// 02.41.06 JR - Feb 10, 2024 - Defect repair: +// - Fix for issue #1057: +// HeMS::CalculateMomentOfInertia() falls back to MainSequence::CalculateMomentOfInertia() +// HeHG::CalculateMomentOfInertia() falls back to GiantBranch::CalculateMomentOfInertia() +// - Added sanity checks for mass and luminosity where necessary in variants of CalculateRadiusOnPhase_Static() +// 02.42.00 JR - Jan 08, 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 +// - added warning for stellar type switch not taken - just a diagnostic for now +// 02.42.01 JR - Jan 21, 2024 - Defect repair +// - fix for issue 1066 - see issue/PR for explanation +// - cleaned up root solvers OmegaAfterSynchronisation(), MassLossToFitInsideRocheLobe(), and Mass0ToMatchDesiredCoreMass(), and their respective functors +// - MassLossToFitInsideRocheLobe(), and Mass0ToMatchDesiredCoreMass() now return -1.0 if no acceptable root found +// - calling code for MassLossToFitInsideRocheLobe() and Mass0ToMatchDesiredCoreMass() now handles -ve return: +// - if MassLossToFitInsideRocheLobe() returns -ve value (i.e. no root found), the binary immediately enters a CE phase +// - if Mass0ToMatchDesiredCoreMass() returns -ve value (i.e. no root found), an arbitrary value is used for core mass (see code for value) +// 02.42.02 RTW - Mar 21, 2024 - Minor edits: +// - Defect repair : Added explicit definition `bool isUnstable = false` to avoid confusion in BaseBinaryStar.cpp +// - Defect repair : Fixed erroneous core mass values in ResolveSNIa in WhiteDwarfs.cpp. Was previously 0 for all core masses. +// - Enhancement: Added output parameter TZAMS for internal variable m_TZAMS -const std::string VERSION_STRING = "02.41.05"; - +const std::string VERSION_STRING = "02.42.02"; # endif // __changelog_h__ diff --git a/src/constants.h b/src/constants.h index 898c00ffa..165f8456f 100755 --- a/src/constants.h +++ b/src/constants.h @@ -180,8 +180,11 @@ extern OBJECT_ID globalObjectId; #undef COMPARE_GLOBAL_TOLERANCE // define/undef this to compare floats with/without tolerance (see FLOAT_TOLERANCE_ABSOLUTE, FLOAT_TOLERANCE_RELATIVE and Compare() function) -constexpr double FLOAT_TOLERANCE_ABSOLUTE = 0.0000005; // Absolute tolerance for floating-point comparisons if COMPARE_GLOBAL_TOLERANCE is defined -constexpr double FLOAT_TOLERANCE_RELATIVE = 0.0000005; // Relative tolerance for floating-point comparisons if COMPARE_GLOBAL_TOLERANCE is defined +constexpr double FLOAT_TOLERANCE_ABSOLUTE = 0.0000005; // absolute tolerance for floating-point comparisons if COMPARE_GLOBAL_TOLERANCE is defined +constexpr double FLOAT_TOLERANCE_RELATIVE = 0.0000005; // relative tolerance for floating-point comparisons if COMPARE_GLOBAL_TOLERANCE is defined + +constexpr double ROOT_ABS_TOLERANCE = 1.0E-10; // absolute tolerance for root finder +constexpr double ROOT_REL_TOLERANCE = 1.0E-6; // relative tolerance for root finder // initialisation constants @@ -216,12 +219,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 @@ -260,7 +265,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) @@ -302,8 +307,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 @@ -330,18 +338,23 @@ constexpr double EPSILON_PULSAR = 1.0; constexpr double MIN_HMXRB_STAR_TO_ROCHE_LOBE_RADIUS_RATIO = 0.8; // Minimum value of stellar radius | Roche Lobe radius for visible HMXRBs constexpr double ADAPTIVE_RLOF_FRACTION_DONOR_GUESS = 0.001; // Fraction of donor mass to use as guess in BaseBinaryStar::MassLossToFitInsideRocheLobe() -constexpr int ADAPTIVE_RLOF_MAX_ITERATIONS = 50; // Maximum number of iterations in BaseBinaryStar::MassLossToFitInsideRocheLobe() -constexpr double ADAPTIVE_RLOF_SEARCH_FACTOR = 2.0; // Search factor in BaseBinaryStar::MassLossToFitInsideRocheLobe() +constexpr int ADAPTIVE_RLOF_MAX_TRIES = 30; // Maximum number of tries in BaseBinaryStar::MassLossToFitInsideRocheLobe() +constexpr int ADAPTIVE_RLOF_MAX_ITERATIONS = 50; // Maximum number of root finder iterations in BaseBinaryStar::MassLossToFitInsideRocheLobe() +constexpr double ADAPTIVE_RLOF_SEARCH_FACTOR_FRAC = 1.0; // Search size factor (fractional part) in BaseBinaryStar::MassLossToFitInsideRocheLobe() (added to 1.0) + +constexpr int ADAPTIVE_MASS0_MAX_TRIES = 30; // Maximum number of tries in HG::Mass0ToMatchDesiredCoreMass() constexpr int ADAPTIVE_MASS0_MAX_ITERATIONS = 50; // Maximum number of iterations in HG::Mass0ToMatchDesiredCoreMass() -constexpr double ADAPTIVE_MASS0_SEARCH_FACTOR = 2.0; // Search factor in HG::Mass0ToMatchDesiredCoreMass() +constexpr double ADAPTIVE_MASS0_SEARCH_FACTOR_FRAC = 1.0; // Search size factor (fractional part) in HG::Mass0ToMatchDesiredCoreMass() (added to 1.0) + +constexpr int TIDES_OMEGA_MAX_TRIES = 30; // Maximum number of tries in BaseBinaryStar::OmegaAfterCircularisation() +constexpr int TIDES_OMEGA_MAX_ITERATIONS = 50; // Maximum number of root finder iterations in BaseBinaryStar::OmegaAfterCircularisation() +constexpr double TIDES_OMEGA_SEARCH_FACTOR_FRAC = 1.0; // Search size factor (fractional part) in BaseBinaryStar::OmegaAfterCircularisation() (added to 1.0) + constexpr double FARMER_PPISN_UPP_LIM_LIN_REGIME = 38.0; // Maximum CO core mass to result in the linear remnant mass regime of the FARMER PPISN prescription constexpr double FARMER_PPISN_UPP_LIM_QUAD_REGIME = 60.0; // Maximum CO core mass to result in the quadratic remnant mass regime of the FARMER PPISN prescription constexpr double FARMER_PPISN_UPP_LIM_INSTABILLITY = 140.0; // Maximum CO core mass to result in PI (upper edge of PISN gap) from FARMER PPISN prescription constexpr double STARTRACK_PPISN_HE_CORE_MASS = 45.0; // Helium core mass remaining following PPISN as assumed in StarTrack (Belczynski et al. 2017 https://arxiv.org/abs/1607.03116) -constexpr int TIDES_OMEGA_MAX_ITERATIONS = 1000; // Maximum number of iterations in BaseBinaryStar::OmegaAfterCircularisation() -constexpr double TIDES_OMEGA_SEARCH_FACTOR = 1.1; // Search factor in BaseBinaryStar::OmegaAfterCircularisation() - // logging constants @@ -400,7 +413,11 @@ enum class BSE_DETAILED_RECORD_TYPE: unsigned int { }; 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 @@ -557,6 +574,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, // file 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 @@ -571,7 +590,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 @@ -588,7 +606,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,10 +630,17 @@ enum class ERROR: int { STELLAR_EVOLUTION_STOPPED, // evolution of current star stopped STELLAR_SIMULATION_STOPPED, // stellar simulation stopped SUGGEST_HELP, // suggest using --help + SWITCH_NOT_TAKEN, // switch to new stellar type not performed 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_MASS0_TRIES, // too many tries in MASS0 root finder TOO_MANY_OMEGA_ITERATIONS, // too many iterations in OMEGA root finder + TOO_MANY_OMEGA_TRIES, // too many tries in OMEGA root finder TOO_MANY_RLOF_ITERATIONS, // too many iterations in RLOF root finder + TOO_MANY_RLOF_TRIES, // too many tries 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 @@ -700,6 +726,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" }}, @@ -714,7 +742,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" }}, @@ -732,6 +759,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" }}, @@ -755,10 +783,17 @@ const COMPASUnorderedMap> ERROR_CATA { ERROR::STELLAR_EVOLUTION_STOPPED, { ERROR_SCOPE::ALWAYS, "Evolution of current star stopped" }}, { 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::SWITCH_NOT_TAKEN, { ERROR_SCOPE::ALWAYS, "Switch to new stellar type not performed" }}, { 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_MASS0_TRIES, { ERROR_SCOPE::ALWAYS, "Reached maximum number of tries 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_OMEGA_TRIES, { ERROR_SCOPE::ALWAYS, "Reached maximum number of tries 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_RLOF_TRIES, { ERROR_SCOPE::ALWAYS, "Reached maximum number of tries 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" }}, @@ -816,6 +851,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, @@ -835,6 +873,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" }, @@ -1909,6 +1950,7 @@ const COMPASUnorderedMap PROPERTY_TYPE_LABEL = { TIMESCALE_MS, \ TOTAL_MASS_AT_COMPACT_OBJECT_FORMATION, \ TRUE_ANOMALY, \ + TZAMS, \ ZETA_HURLEY, \ ZETA_HURLEY_HE, \ ZETA_SOBERMAN, \ @@ -2064,6 +2106,7 @@ const COMPASUnorderedMap STAR_PROPERTY_LABEL = { { STAR_PROPERTY::TIMESCALE_MS, "TIMESCALE_MS" }, { STAR_PROPERTY::TOTAL_MASS_AT_COMPACT_OBJECT_FORMATION, "TOTAL_MASS_AT_COMPACT_OBJECT_FORMATION" }, { STAR_PROPERTY::TRUE_ANOMALY, "TRUE_ANOMALY" }, + { STAR_PROPERTY::TZAMS, "TZAMS" }, { STAR_PROPERTY::ZETA_HURLEY, "ZETA_HURLEY" }, { STAR_PROPERTY::ZETA_HURLEY_HE, "ZETA_HURLEY_HE" }, { STAR_PROPERTY::ZETA_SOBERMAN, "ZETA_SOBERMAN" }, @@ -2820,8 +2863,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; @@ -2834,7 +2877,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 }}, @@ -2852,10 +2895,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 }}, @@ -2946,13 +2989,14 @@ 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 }}, @@ -3042,8 +3086,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 }}, @@ -3188,10 +3232,10 @@ 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 }}, + { PROGRAM_OPTION::MAXIMUM_TIMESTEPS, { TYPENAME::ULONGINT, "Max_Timesteps", "Count", 10, 1 }}, { PROGRAM_OPTION::MCBUR1, { TYPENAME::DOUBLE, "MCBUR1", "Msol", 14, 6 }}, @@ -3213,6 +3257,7 @@ const std::map PROGRAM_OPTION_DETAIL = { { PROGRAM_OPTION::MT_CRIT_MR_GIANT_DEGENERATE_ACCRETOR, { TYPENAME::DOUBLE, "MT_Crit_MR_Giant_Deg_Acc", "-", 14, 6 }}, { PROGRAM_OPTION::MT_CRIT_MR_GIANT_NON_DEGENERATE_ACCRETOR, { TYPENAME::DOUBLE, "MT_Crit_MR_Giant_NonDeg_Acc", "-", 14, 6 }}, { PROGRAM_OPTION::MT_CRIT_MR_HG_DEGENERATE_ACCRETOR, { TYPENAME::DOUBLE, "MT_Crit_MR_HG_Deg_Acc", "-", 14, 6 }}, + { PROGRAM_OPTION::MT_CRIT_MR_HG_NON_DEGENERATE_ACCRETOR, { TYPENAME::DOUBLE, "MT_Crit_MR_HG_NonDeg_Acc", "-", 14, 6 }}, { PROGRAM_OPTION::MT_CRIT_MR_HE_GIANT_DEGENERATE_ACCRETOR, { TYPENAME::DOUBLE, "MT_Crit_MR_HE_Giant_Deg_Acc", "-", 14, 6 }}, { PROGRAM_OPTION::MT_CRIT_MR_HE_GIANT_NON_DEGENERATE_ACCRETOR, { TYPENAME::DOUBLE, "MT_Crit_MR_HE_Giant_NonDeg_Acc", "-", 14, 6 }}, diff --git a/src/utils.cpp b/src/utils.cpp index 9c6507d04..32f0141d4 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -1375,7 +1375,6 @@ namespace utils { * @param [IN] p_A Coefficient of x^2 * @param [IN] p_B Coefficient of x^1 * @param [IN] p_C Coefficient of x^0 (Constant) - * @return Root found (see above) * @return Tuple containing (in order): error value, root found (see above) * The error value returned will be: * ERROR::NONE if no error occurred @@ -1411,6 +1410,25 @@ namespace utils { } + /* + * Tolerance for Boost bracket_and_solve_root() + * + * Determines if the brackets around the root are within the COMPAS defined tolerance. + * + * + * bool BracketTolerance(const double p_Bracket1, const double p_Bracket2) + * + * @param [IN] p_Bracket1 Bracket bound 1 + * @param [IN] p_Bracket2 Bracket bound 2 + * @return Boolean indicating if the brackets bounds are within tolerance + */ + bool BracketTolerance(const double p_Bracket1, const double p_Bracket2) { + double diff = fabs(p_Bracket1 - p_Bracket2); // absolute value of difference + double min = std::min(p_Bracket1, p_Bracket2); // minimum bracket value - could straddle 0.0 + return diff <= ROOT_ABS_TOLERANCE || fabs(diff / min) <= ROOT_REL_TOLERANCE; + } + + /* * Announce COMPAS * @@ -1438,4 +1456,117 @@ 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 + } + + rec = trim(rec); // remove leading and trailing blanks + + if (!(rec.empty() || rec[0] == '#')) { // blank record or comment? + 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 + } + + 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? + 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..b6414c609 100755 --- a/src/utils.h +++ b/src/utils.h @@ -136,8 +136,12 @@ namespace utils { std::tuple SolveQuadratic(const double p_A, const double p_B, double p_C); + bool BracketTolerance(const double p_Bracket1, const double p_Bracket2); + 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 b6e9b3ed5..46987e340 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",