Skip to content

Commit

Permalink
Add gravitational radiation emission to binary evolution.
Browse files Browse the repository at this point in the history
  • Loading branch information
@Adam-Boesky
Adam-Boesky committed Aug 20, 2024
1 parent 370a56f commit d7b90b0
Show file tree
Hide file tree
Showing 11 changed files with 442 additions and 38 deletions.
1 change: 1 addition & 0 deletions compas_python_utils/preprocessing/compasConfigDefault.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ booleanChoices:
# --errors-to-file: False # Default: False
# --evolve-main-sequence-mergers: False # Default: False
# --evolve-unbound-systems: True # Default: True
# --emit-gravitational-radiation: False # Default: False
# --population-data-printing: False # Default: False
# --print-bool-as-string: False # Default: False
# --quiet: False # Default: False
Expand Down
6 changes: 5 additions & 1 deletion online-docs/pages/User guide/Program options/program-options-list-defaults.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -379,6 +379,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 = FALSE

**--enable-warnings** |br|
Display warning messages to stdout. |br|
Default = FALSE
Expand Down Expand Up @@ -1440,7 +1444,7 @@ Go to :ref:`the top of this page <options-props-top>` for the full alphabetical

**Administrative**

--mode, --number-of-systems, --evolve-double-white-dwarfs, --evolve-main-sequence-mergers, --evolve-pulsars, --evolve-unbound-systems, --mass-change-fraction,
--mode, --number-of-systems, --emit-gravitational-radiation, --evolve-double-white-dwarfs, --evolve-main-sequence-mergers, --evolve-pulsars, --evolve-unbound-systems, --mass-change-fraction,
--maximum-evolution-time, --maximum-number-timestep-iterations,
--radial-change-fraction, --random-seed, --timestep-multiplier, --timestep-filename

Expand Down
5 changes: 5 additions & 0 deletions online-docs/pages/whats-new.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,11 @@ Following is a brief list of important updates to the COMPAS code. A complete r

**LATEST RELEASE** |br|

**02.50.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.

**03.00.00 Jul 26, 2024**

This is a major release of COMPAS. There are some significant changes in COMPAS operation and functionality in this release. The major change, and the impetus for
Expand Down
151 changes: 116 additions & 35 deletions src/BaseBinaryStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2490,6 +2490,94 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) {
}


/*
* Calculate and emit gravitational radiation.
*
* This function uses Peters 1964 to approximate the effects of GW emission with two steps:
* - Calculate the change in semi-major axis (m_SemiMajorAxis) per time given by eq 5.6.
* - Calculate the change in eccentricity (m_Eccentricity) per time given by eq 5.7.
*
* m_DaDtGW and m_DeDtGW are updated so that they can be used to calculate the timestep dynamically.
*
* void CalculateGravitationalRadiation()
*/
void 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);
m_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);
m_DeDtGW = (numeratorE / denominatorE) * (1.0 + (121.0 / 304.0) * eccentricitySquared) * YEAR_TO_MYR; // units of Myr^-1
}


/*
* Emit a GW based on the effects calculated by BaseBinaryStar::CalculateGravitationalRadiation().
*
* This function updates the semi-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
*/
void BaseBinaryStar::EmitGravitationalWave(const double p_Dt) {

// Update semimajor axis
double aNew = m_SemiMajorAxis + (m_DaDtGW * p_Dt);
m_SemiMajorAxis = utils::Compare(aNew, 0.0) > 0 ? aNew : 1E-20; // if <0, set to arbitrarily small number

// Update the eccentricity
m_Eccentricity += m_DeDtGW * p_Dt;

// Save values as previous timestep
m_SemiMajorAxisPrev = m_SemiMajorAxis;
m_EccentricityPrev = m_Eccentricity;

// std::cout << "a = " << m_SemiMajorAxis << " da = "<< (m_DaDtGW * p_Dt) << " dt = " << p_Dt << " dadt = " << m_DaDtGW << std::endl;

}


/*
* Choose a timestep based on the parameters of the binary.
*
* This function will return the minimum of (i) a timestep based on the
* orbital timescale of the binary and (ii) (if configured to emit GWs)
* a timestep based on the magnitude of gravitational radiation.
*
* double ChooseTimestep(double p_Dt)
*
* @param [IN] p_Dt previous timestep in Myr
* @return new timestep in Myr
*/
double BaseBinaryStar::ChooseTimestep(double p_Dt) {

p_Dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep()) * OPTIONS->TimestepMultiplier(); // calculate new timestep

if (OPTIONS->EmitGravitationalRadiation()) { // emitting GWs?
p_Dt = std::min(p_Dt, -1.0E-2 * m_SemiMajorAxis / m_DaDtGW); // make dt smaller for strong GWs
}

p_Dt = std::max(std::round(p_Dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, NUCLEAR_MINIMUM_TIMESTEP); // quantised and not less than minimum

return p_Dt;
}


/*
* Evaluate the binary system
*
Expand Down Expand Up @@ -2692,23 +2780,38 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
// evolve the current binary up to the maximum evolution time (and number of steps)

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
unsigned long int stepNum = 0; // initialise step number

while (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // perform binary evolution - iterate over timesteps until told to stop

stepNum++; // increment stepNum

// if user selects to emit GWs, calculate the effects of radiation
// - note that this is placed before ChooseTimestep() is called because
// the timestep is a function of graviational radiation
if (OPTIONS->EmitGravitationalRadiation()) {
CalculateGravitationalRadiation();
}

m_Star2->UpdatePreviousTimestepDuration();
m_Star1->UpdatePreviousTimestepDuration();
if (usingProvidedTimesteps) { // user-provided timesteps?
// select a timestep
// - don't quantise
// - don't apply timestep multiplier
// (we assume user wants the timesteps in the file)
dt = timesteps[stepNum - 1];
}
else { // no - not using user-provided timesteps
dt = ChooseTimestep(dt);
}

error = EvolveOneTimestep(dt); // evolve the binary system one timestep

if (OPTIONS->EmitGravitationalRadiation()) {
EmitGravitationalWave(dt); // emit a graviataional wave
}

if (error != ERROR::NONE) { // SSE error for either constituent star?
evolutionStatus = EVOLUTION_STATUS::SSE_ERROR; // yes - stop evolution
}
Expand Down Expand Up @@ -2795,28 +2898,6 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
SHOW_WARN(ERROR::TIMESTEPS_EXHAUSTED); // show warning
}

if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // continue evolution?

m_Star1->UpdatePreviousTimestepDuration();
m_Star2->UpdatePreviousTimestepDuration();

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
}

if (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?
Expand Down
8 changes: 8 additions & 0 deletions src/BaseBinaryStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -382,6 +382,9 @@ class BaseBinaryStar {
double m_OrbitalAngularMomentumPrev;
double m_OrbitalAngularMomentum;

double m_DaDtGW; // Change in semi-major axis per time due to gravitational radiation
double m_DeDtGW; // Change in eccentricity per time due to gravitational radiation

double m_OrbitalEnergyPrev;
double m_OrbitalEnergy;

Expand Down Expand Up @@ -413,6 +416,11 @@ 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()); }

void CalculateGravitationalRadiation();
void EmitGravitationalWave(const double p_Dt);

double ChooseTimestep(double p_Dt);

void CalculateEnergyAndAngularMomentum();

double CalculateDEccentricityTidalDt(const DBL_DBL_DBL_DBL p_ImKlm, const BinaryConstituentStar* p_Star);
Expand Down
8 changes: 7 additions & 1 deletion src/Options.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -170,7 +170,8 @@ void Options::OptionValues::Initialise() {
m_EvolveMainSequenceMergers = false;
m_EvolvePulsars = false;
m_EvolveUnboundSystems = true;

m_EmitGravitationalRadiation = false;

m_NatalKickForPPISN = false;

m_DetailedOutput = false;
Expand Down Expand Up @@ -814,6 +815,11 @@ bool Options::AddOptions(OptionValues *p_Options, po::options_description *p_Opt
po::value<bool>(&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<bool>(&p_Options->m_EmitGravitationalRadiation)->default_value(p_Options->m_EmitGravitationalRadiation)->implicit_value(false),
("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<bool>(&p_Options->m_ExpelConvectiveEnvelopeAboveLuminosityThreshold)->default_value(p_Options->m_ExpelConvectiveEnvelopeAboveLuminosityThreshold)->implicit_value(true),
Expand Down
4 changes: 4 additions & 0 deletions src/Options.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -338,6 +338,7 @@ class Options {
"evolve-double-white-dwarfs",
"evolve-pulsars",
"evolve-unbound-systems",
"emit-gravitational-radiation",

"initial-mass-1",
"initial-mass-2",
Expand Down Expand Up @@ -457,6 +458,7 @@ class Options {
"evolve-double-white-dwarfs"
"evolve-pulsars",
"evolve-unbound-systems",
"emit-gravitational-radiation",

"fp-error-mode",
"fryer-supernova-engine",
Expand Down Expand Up @@ -696,6 +698,7 @@ class Options {
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_EvolveMainSequenceMergers; // Option to evolve binaries in which two stars merged on the main sequence
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
Expand Down Expand Up @@ -1286,6 +1289,7 @@ class Options {
bool EvolveMainSequenceMergers() const { return OPT_VALUE("evolve-main-sequence-mergers", m_EvolveMainSequenceMergers, 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; }
Expand Down
6 changes: 5 additions & 1 deletion src/changelog.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1255,9 +1255,13 @@
// - Continue evolution of main sequence merger products beyond the main sequence
// - Remove spurious print statement
// - Typo fixes
// 03.01.00 APB - Aug 18, 2024 = Enhancement:
// - Implemented gravitational radiation at each timestep of binary evolution. Available with new '--emit-gravitational-radiation' option.
// - If --emit-gravitational-radiation is true, timestep dynamically as function of da/dt.


const std::string VERSION_STRING = "03.00.02";

const std::string VERSION_STRING = "03.01.00";


# endif // __changelog_h__
Loading

0 comments on commit d7b90b0

Please sign in to comment.