Skip to content

Commit

Permalink
Making dt update earlier and fix some other comments
Browse files Browse the repository at this point in the history
  • Loading branch information
@Adam-Boesky
Adam-Boesky committed May 27, 2024
1 parent cd06201 commit 42fe84e
Show file tree
Hide file tree
Showing 2 changed files with 78 additions and 45 deletions.
117 changes: 72 additions & 45 deletions src/BaseBinaryStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2423,7 +2423,7 @@ void BaseBinaryStar::EmitGravitationalWave(const double p_Dt, const double DaDtG

// 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
m_SemiMajorAxis = utils::Compare(aNew, 0.0) > 0 ? aNew : 1E-20; // if <0, set to arbitrarily small number

// Update the eccentricity
m_Eccentricity += DeDtGW * p_Dt;
Expand All @@ -2434,6 +2434,59 @@ void BaseBinaryStar::EmitGravitationalWave(const double p_Dt, const double DaDtG
}


void BaseBinaryStar::EmitGW2(double p_Dt) {

// 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

// Update semimajor axis and eccentricity
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
m_Eccentricity += m_DeDtGW * p_Dt;

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


double BaseBinaryStar::ChooseTimestep(double p_Dt) {

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


p_Dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep()) * OPTIONS->TimestepMultiplier(); // update dt for orbital timescale
p_Dt = std::round(p_Dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM; // quantised

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

if ((m_Star1->IsOneOf({ STELLAR_TYPE::MASSLESS_REMNANT }) || m_Star2->IsOneOf({ STELLAR_TYPE::MASSLESS_REMNANT })) || p_Dt < NUCLEAR_MINIMUM_TIMESTEP) {
p_Dt = NUCLEAR_MINIMUM_TIMESTEP; // but not less than minimum
}

return p_Dt;
}


/*
* Evaluate the binary system
*
Expand Down Expand Up @@ -2724,7 +2777,7 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {

if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // continue evolution
// 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
Expand All @@ -2738,15 +2791,28 @@ 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
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 (stepNum > 1) { // do not update first timestep
if (usingProvidedTimesteps) { // user-provided timesteps?
dt = timesteps[stepNum - 1]; // yes - select from given dt values
}
else {
dt = ChooseTimestep(dt); // no - dynamically timestep
}
}

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

// emit gravitational radiation if selected by user
if (OPTIONS->EmitGravitationalRadiation()) {
EmitGW2(dt);
}

// check for problems
if (m_Error != ERROR::NONE) { // SSE error for either constituent star?
evolutionStatus = EVOLUTION_STATUS::SSE_ERROR; // yes - stop evolution
Expand Down Expand Up @@ -2775,10 +2841,6 @@ 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
Expand Down Expand Up @@ -2864,46 +2926,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
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 (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) {
dt = NUCLEAR_MINIMUM_TIMESTEP; // but not less than minimum
}
stepNum++; // increment stepNum
}
}

if (usingProvidedTimesteps && timesteps.size() > stepNum) { // all user-defined timesteps consumed?
Expand Down
6 changes: 6 additions & 0 deletions src/BaseBinaryStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -383,6 +383,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 @@ -418,6 +421,9 @@ class BaseBinaryStar {

std::tuple<double, double> CalculateGravitationalRadiation();
void EmitGravitationalWave(const double p_Dt, const double DaDtGW, const double DeDtGW);
void EmitGW2(double p_Dt);

double ChooseTimestep(double p_Dt);

void CalculateEnergyAndAngularMomentum();

Expand Down

0 comments on commit 42fe84e

Please sign in to comment.