Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Issue 1034 - angular momentum not conserved in tides code #1041

Merged
merged 19 commits into from
Dec 28, 2023
Merged

Issue 1034 - angular momentum not conserved in tides code #1041

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
19 commits
Select commit Hold shift + click to select a range
a3b1fac
Defect repair:
jeffriley Dec 27, 2023
31b53ed
Updated changelog to indicate some unused functions removed
jeffriley Dec 27, 2023
0ff9d12
Merge pull request #1040 from jeffriley/issue-1034
jeffriley Dec 27, 2023
1073ba6
Changed CalculateMomentOfInertia() to CalculateMomentOfInertiaAU() in…
jeffriley Dec 27, 2023
6def95c
Changed CalculateMomentOfInertia() to CalculateMomentOfInertiaAU() in…
jeffriley Dec 27, 2023
5167070
Merge pull request #1042 from jeffriley/issue-1034
jeffriley Dec 27, 2023
962a65c
Fixed typo
jeffriley Dec 27, 2023
d64b8d5
Merge pull request #1043 from jeffriley/issue-1034
jeffriley Dec 27, 2023
5154502
Removed function CalculateGyrationRadius() and all calls to it
jeffriley Dec 28, 2023
bc96fe5
Removed diagnostics...
jeffriley Dec 28, 2023
ac9285b
Merge pull request #1044 from jeffriley/issue-1034
jeffriley Dec 28, 2023
6dac310
Fixed units
jeffriley Dec 28, 2023
133e4be
Merge pull request #1045 from jeffriley/issue-1034
jeffriley Dec 28, 2023
3a94267
Changed CalculateMomentOfInertia() to properly implement Hurley et al…
jeffriley Dec 28, 2023
1354bbf
Update What's New
jeffriley Dec 28, 2023
18217b2
Update What's New
jeffriley Dec 28, 2023
2624083
Merge pull request #1046 from jeffriley/issue-1034
jeffriley Dec 28, 2023
175fccb
Requested changes
jeffriley Dec 28, 2023
12f58ea
Merge pull request #1047 from jeffriley/issue-1034
jeffriley Dec 28, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions online-docs/pages/whats-new.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,12 @@ Following is a brief list of important updates to the COMPAS code. A complete r

**LATEST RELEASE** |br|

**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.
* Changed CalculateMomentOfInertia() to properly implement Hurley et al., 2000 eq 109.
* This change may change DCO yields slightly when compared to previous versions of the code.

**02.41.00 Nov 02, 2023**

* Added a naive tides implementation.
Expand Down
5 changes: 1 addition & 4 deletions src/BH.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,14 +50,11 @@ class BH: virtual public BaseStar, public Remnants {

double CalculateEddingtonCriticalRate() const { return 2.6E-8 * m_Mass * MYR_TO_YEAR; } // E.g., Marchant+, 2017, Eq. 3, assuming accretion efficiency of 10%

double CalculateGyrationRadius() const { return 0.0; } // No tidal coupling to a BH

double CalculateLuminosityOnPhase() const { return CalculateLuminosityOnPhase_Static(); }

double CalculateMassLossRate() { return 0.0; } // Ensure that BHs don't lose mass in winds

double CalculateMomentOfInertia(const double p_RemnantRadius = 0.0) const { return (2.0 / 5.0) * m_Mass * m_Radius * m_Radius; }
double CalculateMomentOfInertiaAU(const double p_RemnantRadius = 0.0) const { return CalculateMomentOfInertia(p_RemnantRadius * RSOL_TO_AU) * RSOL_TO_AU * RSOL_TO_AU; }
double CalculateMomentOfInertia() const { return (2.0 / 5.0) * m_Mass * m_Radius * m_Radius; }

double CalculateRadiusOnPhase() const { return CalculateRadiusOnPhase_Static(m_Mass); } // Use class member variables - returns radius in Rsol
};
Expand Down
112 changes: 36 additions & 76 deletions src/BaseBinaryStar.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -272,12 +272,12 @@ void BaseBinaryStar::SetRemainingValues() {
m_SemiMajorAxisAtDCOFormation = DEFAULT_INITIAL_DOUBLE_VALUE;
m_EccentricityAtDCOFormation = DEFAULT_INITIAL_DOUBLE_VALUE;

double gyrationRadius1 = m_Star1->CalculateGyrationRadius();
double gyrationRadius2 = m_Star2->CalculateGyrationRadius();
double momentOfInertia1 = m_Star1->CalculateMomentOfInertiaAU();
double momentOfInertia2 = m_Star2->CalculateMomentOfInertiaAU();

m_TotalEnergy = CalculateTotalEnergy(m_SemiMajorAxis, m_Star1->Mass(), m_Star2->Mass(), m_Star1->RZAMS(), m_Star2->RZAMS(), m_Star1->Omega(), m_Star2->Omega(), gyrationRadius1, gyrationRadius2);
m_TotalEnergy = CalculateTotalEnergy(m_SemiMajorAxis, m_Star1->Mass(), m_Star2->Mass(), m_Star1->Omega(), m_Star2->Omega(), momentOfInertia1, momentOfInertia2);

m_TotalAngularMomentum = CalculateAngularMomentum(m_SemiMajorAxis, m_Eccentricity, m_Star1->Mass(), m_Star2->Mass(), m_Star1->RZAMS(), m_Star2->RZAMS(), m_Star1->Omega(), m_Star2->Omega(), gyrationRadius1, gyrationRadius2);
m_TotalAngularMomentum = CalculateAngularMomentum(m_SemiMajorAxis, m_Eccentricity, m_Star1->Mass(), m_Star2->Mass(), m_Star1->Omega(), m_Star2->Omega(), momentOfInertia1, momentOfInertia2);
m_TotalAngularMomentumPrev = m_TotalAngularMomentum;

m_Omega = 0.0;
Expand Down Expand Up @@ -2046,51 +2046,32 @@ void BaseBinaryStar::InitialiseMassTransfer() {
* double CalculateTotalEnergy(const double p_SemiMajorAxis,
* const double p_Star1Mass,
* const double p_Star2Mass,
* const double p_Star1Radius,
* const double p_Star2Radius,
* const double p_Star1_SpinAngularVelocity,
* const double p_Star1_SpinAngularVelocity,
* const double p_Star1_GyrationRadius,
* const double p_Star2_GyrationRadius)
* const double p_Star1SpinAngularVelocity,
* const double p_Star2SpinAngularVelocity,
* const double p_Star1MomentOfInertia,
* const double p_Star2MomentOfInertia)
*
* @param [IN] p_SemiMajorAxis Semi-major axis of the binary
* @param [IN] p_Star1Mass Mass of star 1
* @param [IN] p_Star2Mass Mass of star 2
* @param [IN] p_Star1Radius Radius of star 1
* @param [IN] p_Star2Radius Radius of star 2
* @param [IN] p_Star1_SpinAngularVelocity Spin angular velocity of star 1
* @param [IN] p_Star1_SpinAngularVelocity Spin angular velocity of star 1
* @param [IN] p_Star1_GyrationRadius Gyration radius of star 1
* @param [IN] p_Star2_GyrationRadius Gyration radius of star 2
* @param [IN] p_Star1SpinAngularVelocity Spin angular velocity of star 1
* @param [IN] p_Star2SpinAngularVelocity Spin angular velocity of star 2
* @param [IN] p_Star1MomentOfInertia Moment of inertia of star 1
* @param [IN] p_Star2MomentOfInertia Moment of inertia of star 2
* @return Total energy of the binary
*/
double BaseBinaryStar::CalculateTotalEnergy(const double p_SemiMajorAxis,
const double p_Star1Mass,
const double p_Star2Mass,
const double p_Star1Radius,
const double p_Star2Radius,
const double p_Star1_SpinAngularVelocity,
const double p_Star2_SpinAngularVelocity,
const double p_Star1_GyrationRadius,
const double p_Star2_GyrationRadius) const {
double m1 = p_Star1Mass;
double m2 = p_Star2Mass;
const double p_Star1SpinAngularVelocity,
const double p_Star2SpinAngularVelocity,
const double p_Star1MomentOfInertia,
const double p_Star2MomentOfInertia) const {

double R1 = p_Star1Radius;
double R2 = p_Star2Radius;
double w1_2 = p_Star1SpinAngularVelocity * p_Star1SpinAngularVelocity;
double w2_2 = p_Star2SpinAngularVelocity * p_Star2SpinAngularVelocity;

double w1 = p_Star1_SpinAngularVelocity;
double w2 = p_Star2_SpinAngularVelocity;

double ks1 = p_Star1_GyrationRadius;
double ks2 = p_Star2_GyrationRadius;

constexpr double RSOL_TO_AU_2 = RSOL_TO_AU * RSOL_TO_AU;

double Is1 = ks1 * m1 * R1 * R1 * RSOL_TO_AU_2;
double Is2 = ks2 * m2 * R2 * R2 * RSOL_TO_AU_2;

return (0.5 * Is1 * w1 * w1) + (0.5 * Is2 * w2 * w2) - (0.5 * G_AU_Msol_yr * m1 * m2 / p_SemiMajorAxis);
return 0.5 * ((p_Star1MomentOfInertia * w1_2) + (p_Star2MomentOfInertia * w2_2) - (G_AU_Msol_yr * p_Star1Mass * p_Star2Mass / p_SemiMajorAxis));
}


Expand All @@ -2101,57 +2082,36 @@ double BaseBinaryStar::CalculateTotalEnergy(const double p_SemiMajorAxis,
*
*
* double CalculateAngularMomentum(const double p_SemiMajorAxis,
* const double p_Eccentricity,
* const double p_Star1Mass,
* const double p_Star2Mass,
* const double p_Star1Radius,
* const double p_Star2Radius,
* const double p_Star1_SpinAngularVelocity,
* const double p_Star1_SpinAngularVelocity,
* const double p_Star1_GyrationRadius,
* const double p_Star2_GyrationRadius)
* const double p_Eccentricity,
* const double p_Star1Mass,
* const double p_Star2Mass,
* const double p_Star1SpinAngularVelocity,
* const double p_Star1SpinAngularVelocity,
* const double p_Star1MomentOfInertia,
* const double p_Star2MomentOfInertia)
*
* @param [IN] p_SemiMajorAxis Semi-major axis of the binary
* @param [IN] p_Eccentricity Eccentricity of the binary
* @param [IN] p_Star1Mass Mass of the primary
* @param [IN] p_Star2Mass Mass of the secondary
* @param [IN] p_Star1Radius Radius of the primary
* @param [IN] p_Star2Radius Radius of the secondary
* @param [IN] p_Star1_SpinAngularVelocity Orbital frequency of the primary
* @param [IN] p_Star1_SpinAngularVelocity Orbital frequency of the secondary
* @param [IN] p_Star1_GyrationRadius Gyration radius of the primary
* @param [IN] p_Star2_GyrationRadius Gyration radius of the secondary
* @param [IN] p_Star1SpinAngularVelocity Orbital frequency of the primary
* @param [IN] p_Star1SpinAngularVelocity Orbital frequency of the secondary
* @param [IN] p_Star1MomentOfInertia Moment of inertia of the primary
* @param [IN] p_Star2MomentOfInertia Moment of inertia of the secondary
* @return Angular momentum of the binary
*/
double BaseBinaryStar::CalculateAngularMomentum(const double p_SemiMajorAxis,
const double p_Eccentricity,
const double p_Star1Mass,
const double p_Star2Mass,
const double p_Star1Radius,
const double p_Star2Radius,
const double p_Star1_SpinAngularVelocity,
const double p_Star2_SpinAngularVelocity,
const double p_Star1_GyrationRadius,
const double p_Star2_GyrationRadius) const {

double m1 = p_Star1Mass;
double m2 = p_Star2Mass;

double R1 = p_Star1Radius * RSOL_TO_AU;
double R2 = p_Star2Radius * RSOL_TO_AU;

double w1 = p_Star1_SpinAngularVelocity;
double w2 = p_Star2_SpinAngularVelocity;

double ks1 = p_Star1_GyrationRadius;
double ks2 = p_Star2_GyrationRadius;

double Is1 = ks1 * m1 * R1 * R1;
double Is2 = ks2 * m2 * R2 * R2;
const double p_Star1SpinAngularVelocity,
const double p_Star2SpinAngularVelocity,
const double p_Star1MomentOfInertia,
const double p_Star2MomentOfInertia) const {

double Jorb = CalculateOrbitalAngularMomentum(m1, m2, p_SemiMajorAxis, p_Eccentricity);
double Jorb = CalculateOrbitalAngularMomentum(p_Star1Mass, p_Star2Mass, p_SemiMajorAxis, p_Eccentricity);

return (Is1 * w1) + (Is2 * w2) + Jorb;
return (p_Star1MomentOfInertia * p_Star1SpinAngularVelocity) + (p_Star2MomentOfInertia * p_Star2SpinAngularVelocity) + Jorb;
}


Expand Down
54 changes: 13 additions & 41 deletions src/BaseBinaryStar.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -407,35 +407,16 @@ class BaseBinaryStar {

// member functions - alphabetical in groups (sort of...)

// CalculateAngularMomentum - the actual function takes 10 parameters because of the various calling permutations
// - various signatures are defined here - they just assemble the parameters as required
// and call the actual function
// JR: todo: note in the original code the binary orbital velicity was passed in as a parameter but never used - I removed it

void SetInitialValues(const unsigned long int p_Seed, const long int p_Id);
void SetRemainingValues();

double CalculateAngularMomentum(const double p_SemiMajorAxis,
const double p_Eccentricity,
const double p_Star1Mass,
const double p_Star2Mass,
const double p_Star1Radius,
const double p_Star2Radius,
const double p_Star1_SpinAngularVelocity,
const double p_Star2_SpinAngularVelocity,
const double p_Star1_GyrationRadius,
const double p_Star2_GyrationRadius) const;

double CalculateAngularMomentum() const { return CalculateAngularMomentum(m_SemiMajorAxis, m_Eccentricity, m_Star1->Mass(), m_Star2->Mass(), m_Star1->Radius(), m_Star2->Radius(), m_Star1->Omega(), m_Star2->Omega(), m_Star1->CalculateGyrationRadius(), m_Star2->CalculateGyrationRadius()); }

double CalculateAngularMomentum(const double p_SemiMajorAxis,
const double p_Eccentricity,
const double p_Star1_SpinAngularVelocity,
const double p_Star2_SpinAngularVelocity,
const double p_Star1_GyrationRadius,
const double p_Star2_GyrationRadius) const { return CalculateAngularMomentum(p_SemiMajorAxis, p_Eccentricity, m_Star1->Mass(), m_Star2->Mass(), m_Star1->Radius(), m_Star2->Radius(), p_Star1_SpinAngularVelocity, p_Star2_SpinAngularVelocity, p_Star1_GyrationRadius, p_Star2_GyrationRadius); }
const double p_Star1SpinAngularVelocity,
const double p_Star2SpinAngularVelocity,
const double p_Star1MomentOfInertia,
const double p_Star2MomentOfInertia) const;

double CalculateAngularMomentumPrev() const { return CalculateAngularMomentum(m_SemiMajorAxisPrev, m_EccentricityPrev, m_Star1->MassPrev(), m_Star2->MassPrev(), m_Star1->RadiusPrev(), m_Star2->RadiusPrev(), m_Star1->OmegaPrev(), m_Star2->OmegaPrev(), m_Star1->CalculateGyrationRadius(), m_Star2->CalculateGyrationRadius()); }
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 CalculateEnergyAndAngularMomentum();

Expand Down Expand Up @@ -466,27 +447,15 @@ class BaseBinaryStar {

double CalculateTimeToCoalescence(double a0, double e0, double m1, double m2) const;

// CalculateTotalEnergy - the actual function takes 9 parameters because of the various calling permutations
// - various signatures are defined here - they just assemble the parameters as required
// and call the actual function
double CalculateTotalEnergy(const double p_SemiMajorAxis,
const double p_Star1Mass,
const double p_Star2Mass,
const double p_Star1Radius,
const double p_Star2Radius,
const double p_Star1_SpinAngularVelocity,
const double p_Star2_SpinAngularVelocity,
const double p_Star1GyrationRadius,
const double p_Star2GyrationRadius) const;

double CalculateTotalEnergy() const { return CalculateTotalEnergy(m_SemiMajorAxis, m_Star1->Mass(), m_Star2->Mass(), m_Star1->Radius(), m_Star2->Radius(), m_Star1->Omega(), m_Star2->Omega(), m_Star1->CalculateGyrationRadius(), m_Star2->CalculateGyrationRadius()); }

double CalculateTotalEnergy(const double p_SemiMajorAxis,
const double p_Star1_SpinAngularVelocity,
const double p_Star2_SpinAngularVelocity,
const double p_Star1_GyrationRadius,
const double p_Star2_GyrationRadius) const { return CalculateTotalEnergy(p_SemiMajorAxis, m_Star1->Mass(), m_Star2->Mass(), m_Star1->Radius(), m_Star2->Radius(), p_Star1_SpinAngularVelocity, p_Star2_SpinAngularVelocity, p_Star1_GyrationRadius, p_Star2_GyrationRadius); }
const double p_Star1SpinAngularVelocity,
const double p_Star2SpinAngularVelocity,
const double p_Star1MomentOfInertia,
const double p_Star2MomentOfInertia) const;

double CalculateTotalEnergy() const { return CalculateTotalEnergy(m_SemiMajorAxis, m_Star1->Mass(), m_Star2->Mass(), m_Star1->Omega(), m_Star2->Omega(), m_Star1->CalculateMomentOfInertiaAU(), m_Star2->CalculateMomentOfInertiaAU()); }

void EvaluateBinary(const double p_Dt);

Expand All @@ -500,6 +469,9 @@ class BaseBinaryStar {
void ResolveMassChanges();
bool ResolveSupernova();

void SetInitialValues(const unsigned long int p_Seed, const long int p_Id);
void SetRemainingValues();

void SetPostCEEValues(const double p_SemiMajorAxis,
const double p_Eccentricity,
const double p_RocheLobe1to2,
Expand Down
Loading