diff --git a/compas_python_utils/preprocessing/compasConfigDefault.yaml b/compas_python_utils/preprocessing/compasConfigDefault.yaml index 2264d93cc..382502ed4 100644 --- a/compas_python_utils/preprocessing/compasConfigDefault.yaml +++ b/compas_python_utils/preprocessing/compasConfigDefault.yaml @@ -14,6 +14,7 @@ booleanChoices: # --enable-warnings: False # Default: False # option to enable/disable warning messages # --errors-to-file: False # Default: False # --evolve-unbound-systems: True # Default: True +# --emit-gravitational-radiation: True # Default: True # --population-data-printing: False # Default: False # --print-bool-as-string: False # Default: False # --quiet: False # Default: False diff --git a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst index 95de10ff1..60c9bb4bb 100644 --- a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst +++ b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst @@ -349,6 +349,10 @@ Default = 0.0 Multiplication factor for Eddington accretion for NS & BH (i.e. > 1 is super-eddington and 0 is no accretion). |br| Default = 1.0 +**--emit-gravitational-radiation** |br| +Emit gravitational radiation at each timestep of binary evolution according to Peters 1964. |br| +Default = TRUE + **--enable-warnings** |br| Display warning messages to stdout. |br| Default = FALSE @@ -1370,7 +1374,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, --mass-change-fraction, --maximum-evolution-time, --maximum-number-timestep-iterations, +--mode, --number-of-systems, --emit-gravitational-radiation, --evolve-double-white-dwarfs, --evolve-pulsars, --evolve-unbound-systems, --mass-change-fraction, --maximum-evolution-time, --maximum-number-timestep-iterations, --radial-change-fraction, --random-seed, --timestep-multiplier, --timestep-filename --grid, --grid-start-line, --grid-lines-to-process diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index d43f5cada..155e29c8a 100644 --- a/online-docs/pages/whats-new.rst +++ b/online-docs/pages/whats-new.rst @@ -6,6 +6,11 @@ Following is a brief list of important updates to the COMPAS code. A complete r **LATEST RELEASE** |br| +**02.47.00 May 17, 2024** + +* New option to emit gravitational radiation at each timestep of binary evolution: ``--emit-gravitational-radiation``. The effects of radiation are approximated by the change in semimajor axis and eccentricity from Peters 1946 equations 5.6 and 5.7. +* If ``--emit-gravitational-radiation``, we update timesteps dynamically as a function of :math:`da/dt` if gravitational radiation is large enough. + **02.46.00 May 13, 2024** * added options ``--radial-change-fraction`` and ``--mass-change-fraction``, as approximate desired fractional changes in stellar radius and mass on phase when setting SSE and BSE timesteps diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index fd75fb696..7a646235a 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -2370,6 +2370,69 @@ void BaseBinaryStar::ResolveMassChanges() { } +/* + * Update the GW effects on the binary + * + * Use Peters 1964's approximations of the effects of emitting a GW by: + * - Updating the binary's semi-major axis (eq 5.6) + * - Updating the binary's eccentricity (eq 5.7) + * + * void CalculateGravitationalRadiation() + * + * @return The change in semi major axis and eccentricity per time due to GW emission +*/ +std::tuple BaseBinaryStar::CalculateGravitationalRadiation() { + + // Useful values + double eccentricitySquared = m_Eccentricity * m_Eccentricity; + double oneMinusESq = 1.0 - eccentricitySquared; + double oneMinusESq_5 = oneMinusESq * oneMinusESq * oneMinusESq * oneMinusESq * oneMinusESq; + double G_AU_Msol_yr_3 = G_AU_Msol_yr * G_AU_Msol_yr * G_AU_Msol_yr; + double C_AU_Yr_5 = C_AU_yr * C_AU_yr * C_AU_yr * C_AU_yr * C_AU_yr; + double m_SemiMajorAxis_3 = m_SemiMajorAxis * m_SemiMajorAxis * m_SemiMajorAxis; + double massAndGAndCTerm = G_AU_Msol_yr_3 * m_Star1->Mass() * m_Star2->Mass() * (m_Star1->Mass() + m_Star2->Mass()) / C_AU_Yr_5; // G^3 * m1 * m2(m1 + m2) / c^5 in units of Msol, AU and yr + + // Approximate rate of change in semimajor axis + double numeratorA = -64.0 * massAndGAndCTerm; + double denominatorA = 5.0 * m_SemiMajorAxis_3 * std::sqrt(oneMinusESq_5 * oneMinusESq * oneMinusESq); + double DaDtGW = (numeratorA / denominatorA) * (1.0 + (73.0 / 24.0) * eccentricitySquared + (37.0 / 96.0) * eccentricitySquared * eccentricitySquared) * MYR_TO_YEAR; // units of AU Myr^-1 + + // Approximate rate of change in eccentricity + double numeratorE = -304.0 * m_Eccentricity * massAndGAndCTerm; + double denominatorE = 15.0 * m_SemiMajorAxis_3 * m_SemiMajorAxis * std::sqrt(oneMinusESq_5); + double DeDtGW = (numeratorE / denominatorE) * (1.0 + (121.0 / 304.0) * eccentricitySquared) * YEAR_TO_MYR; // units of Myr^-1 + + return std::make_tuple(DaDtGW, DeDtGW); +} + + +/* + * Emit a GW based on the effects calculated by BaseBinaryStar::CalculateGravitationalRadiation(). + * + * This function updates the sem-major axis, eccentricity, and previous eccentricity values + * (m_SemiMajorAxis, m_Eccentricity, and m_EccentricityPrev) as a result of emitting the GW. + * + * void EmitGravitationalRadiation(const double p_Dt) + * + * @param [IN] p_Dt timestep in Myr + * @param [IN] DaDtGW change in semimajor axis per time in AU/Myr due to GW radiation + * @param [IN] DeDtGW change in eccentricity per time in Myr^-1 due to GW radiation +*/ +void BaseBinaryStar::EmitGravitationalWave(const double p_Dt, const double DaDtGW, const double DeDtGW) { + + // Update semimajor axis + double aNew = m_SemiMajorAxis + (DaDtGW * p_Dt); + m_SemiMajorAxis = aNew < 0.0 ? 1E-20 : (aNew); // if <0, set to arbitrarily small number + + // Update the eccentricity + m_Eccentricity += DeDtGW * p_Dt; + + // Save values as previous timestep + m_SemiMajorAxisPrev = m_SemiMajorAxis; + m_EccentricityPrev = m_Eccentricity; +} + + /* * Evaluate the binary system * @@ -2674,6 +2737,9 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() { dt = std::max(std::round(dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, NUCLEAR_MINIMUM_TIMESTEP); // quantised } + double DaDtGW; + double DeDtGW; + unsigned long int stepNum = 1; // initialise step number while (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // perform binary evolution - iterate over timesteps until told to stop @@ -2708,6 +2774,10 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() { if (OPTIONS->RLOFPrinting()) StashRLOFProperties(MASS_TRANSFER_TIMING::PRE_MT); // stash properties immediately pre-Mass Transfer + if (OPTIONS->EmitGravitationalRadiation()) { // emit gravitational wave if desired + EmitGravitationalWave(dt, DaDtGW, DeDtGW); // this may update m_SemiMajorAxis, m_Eccentricity, and m_EccentricityPrev + } + EvaluateBinary(dt); // evaluate the binary at this timestep (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_BINARY_TIMESTEP); // print (log) detailed output @@ -2793,6 +2863,11 @@ 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 (evolutionStatus == EVOLUTION_STATUS::STELLAR_MERGER || evolutionStatus == EVOLUTION_STATUS::STARS_TOUCHING) { // are the stars touching/merging? + m_MassTransferTrackerHistory = MT_TRACKING::MERGER; // set flag for merger + (void)PrintCommonEnvelope(); // print (log) common envelope details + } + 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 @@ -2814,6 +2889,13 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() { 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 + + if (OPTIONS->EmitGravitationalRadiation()) { // if our binary is emitting GWs + std::tie(DaDtGW, DeDtGW) = CalculateGravitationalRadiation(); // update the GW radiation effects da/dt and de/dt + if (!usingProvidedTimesteps) { + dt = std::min(dt, std::round(-1E-2 * m_SemiMajorAxis / DaDtGW / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM); // make dt smaller if strong GW effects + } + } } if ((m_Star1->IsOneOf({ STELLAR_TYPE::MASSLESS_REMNANT }) || m_Star2->IsOneOf({ STELLAR_TYPE::MASSLESS_REMNANT })) || dt < NUCLEAR_MINIMUM_TIMESTEP) { diff --git a/src/BaseBinaryStar.h b/src/BaseBinaryStar.h index 0c8f32420..dfda4787a 100644 --- a/src/BaseBinaryStar.h +++ b/src/BaseBinaryStar.h @@ -416,6 +416,9 @@ class BaseBinaryStar { double CalculateAngularMomentum() const { return CalculateAngularMomentum(m_SemiMajorAxis, m_Eccentricity, m_Star1->Mass(), m_Star2->Mass(), m_Star1->Omega(), m_Star2->Omega(), m_Star1->CalculateMomentOfInertiaAU(), m_Star2->CalculateMomentOfInertiaAU()); } + std::tuple CalculateGravitationalRadiation(); + void EmitGravitationalWave(const double p_Dt, const double DaDtGW, const double DeDtGW); + void CalculateEnergyAndAngularMomentum(); double CalculateDEccentricityTidalDt(const DBL_DBL_DBL_DBL p_ImKlm, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity); diff --git a/src/Options.cpp b/src/Options.cpp index 05e9c279a..f4bfe4075 100644 --- a/src/Options.cpp +++ b/src/Options.cpp @@ -161,7 +161,8 @@ void Options::OptionValues::Initialise() { m_EvolveDoubleWhiteDwarfs = false; m_EvolvePulsars = false; m_EvolveUnboundSystems = true; - + m_EmitGravitationalRadiation = true; + m_NatalKickForPPISN = false; m_DetailedOutput = false; @@ -807,6 +808,11 @@ bool Options::AddOptions(OptionValues *p_Options, po::options_description *p_Opt po::value(&p_Options->m_EvolveUnboundSystems)->default_value(p_Options->m_EvolveUnboundSystems)->implicit_value(true), ("Continue evolving stars even if the binary is disrupted (default = " + std::string(p_Options->m_EvolveUnboundSystems ? "TRUE" : "FALSE") + ")").c_str() ) + ( + "emit-gravitational-radiation", + po::value(&p_Options->m_EmitGravitationalRadiation)->default_value(p_Options->m_EmitGravitationalRadiation)->implicit_value(true), + ("Emit gravitational radiation at each timestep of binary evolution (default = " + std::string(p_Options->m_EmitGravitationalRadiation ? "TRUE" : "FALSE") + ")").c_str() + ) ( "expel-convective-envelope-above-luminosity-threshold", po::value(&p_Options->m_ExpelConvectiveEnvelopeAboveLuminosityThreshold)->default_value(p_Options->m_ExpelConvectiveEnvelopeAboveLuminosityThreshold)->implicit_value(true), diff --git a/src/Options.h b/src/Options.h index 43496bf4d..342acbc3d 100755 --- a/src/Options.h +++ b/src/Options.h @@ -336,6 +336,7 @@ class Options { "evolve-double-white-dwarfs", "evolve-pulsars", "evolve-unbound-systems", + "emit-gravitational-radiation", "initial-mass-1", "initial-mass-2", @@ -450,6 +451,7 @@ class Options { "evolve-double-white-dwarfs" "evolve-pulsars", "evolve-unbound-systems", + "emit-gravitational-radiation", "fryer-supernova-engine", @@ -674,6 +676,7 @@ class Options { bool m_EvolvePulsars; // Whether to evolve pulsars or not bool m_NatalKickForPPISN; // Flag if PPISN remnant should receive a non-zero natal kick bool m_EvolveUnboundSystems; // Option to chose if unbound systems are evolved until death or the evolution stops after the system is unbound during a SN. + bool m_EmitGravitationalRadiation; // Option to emit gravitational radiation for each timestep of binary evolution bool m_DetailedOutput; // Print detailed output details to file (default = false) bool m_PopulationDataPrinting; // Print certain data for small populations, but not for larger one @@ -1265,6 +1268,7 @@ class Options { bool EvolveDoubleWhiteDwarfs() const { return OPT_VALUE("evolve-double-white-dwarfs", m_EvolveDoubleWhiteDwarfs, true); } bool EvolvePulsars() const { return OPT_VALUE("evolve-pulsars", m_EvolvePulsars, true); } bool EvolveUnboundSystems() const { return OPT_VALUE("evolve-unbound-systems", m_EvolveUnboundSystems, true); } + bool EmitGravitationalRadiation() const { return OPT_VALUE("emit-gravitational-radiation", m_EmitGravitationalRadiation, true); } bool ExpelConvectiveEnvelopeAboveLuminosityThreshold() const { return OPT_VALUE("expel-convective-envelope-above-luminosity-threshold", m_ExpelConvectiveEnvelopeAboveLuminosityThreshold, true); } bool FixedRandomSeedCmdLine() const { return m_CmdLine.optionValues.m_FixedRandomSeed; } diff --git a/src/changelog.h b/src/changelog.h index 117ef8fdf..484b02078 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1181,8 +1181,11 @@ // - fix for issue #744 - GB parameters `p` and `q` calculated differently for naked helium stars (see issue for details) // - changed name of `ResolveEnvelopeLoss()` parameter `p_NoCheck` to `p_Force` (it is required, and now we understand why... see issue #873) // - some code cleanup +// 02.47.00 APB - May 17, 2024 = Enhancement: +// - Implemented gravitational radiation at each timestep of binary evolution. Available with new '--emit-gravitational-radiation' option. +// - If --emit-gravitational-radiation, timestep dynamically as function of da/dt. -const std::string VERSION_STRING = "02.46.05"; +const std::string VERSION_STRING = "02.47.00"; # endif // __changelog_h__ diff --git a/src/constants.h b/src/constants.h index 7ccfdfbe6..6611bc53c 100755 --- a/src/constants.h +++ b/src/constants.h @@ -273,6 +273,7 @@ constexpr double LSOLW = 3.844E26; constexpr double AU = 149597870700.0; // 1 AU (Astronomical Unit) in metres constexpr double KM = 1000.0; // 1 km (Kilometre) in metres constexpr double C = 3.0E8; // Speed of light in m s^-1 +constexpr double C_AU_yr = C * 0.001 * KM_TO_AU * SECONDS_IN_YEAR; // Speed of light in AU yr^-1 constexpr double MU_0 = 4.0 * M_PI * 1.0E-7; // Vacuum permeability in m kg s-2 A-2 diff --git a/src/yaml.h b/src/yaml.h index b8b25ee5a..ee57e5c2b 100644 --- a/src/yaml.h +++ b/src/yaml.h @@ -80,6 +80,7 @@ namespace yaml { " --enable-warnings # option to enable/disable warning messages", " --errors-to-file", " --evolve-unbound-systems", + " --emit-gravitational-radiation", " --population-data-printing", " --print-bool-as-string", " --quiet",