Naive tides implementation

[jeffriley]

Enhancement, a little cleanup:

• Added naive tides implementation. Functionality enabled with new option --enable-tides. Default is no tides.
• Fixed CalculateOrbitalAngularMomentum() (now uses eccentricity)
• Added links to online documentation to splash string

[avivajpeyi]

@jeffriley i've added the changes to the pytest to use the fiducial-bbh config (and not use the method-paper config)

[avivajpeyi]

Err, i think i messed the config up.

Commandline Options error: Unknown Mass Loss Prescription

https://github.com/TeamCOMPAS/COMPAS/actions/runs/6702099754/job/18210562003?pr=1012#step:6:354

[avivajpeyi]

CI/Pytest changes

• Previously, we were using the 'Methods paper detailed evolution diagram' config (a COMPAS config that generated a DCO that merged within hubble time)
• Now, the new config creates a DCO but this DCO does not merge within hubble time. We needed to make some changes to the config to allow for this new tides option (right, @jeffriley ?)

[jeffriley]

BaseBinaryStar.cpp, dealing with CHE stars (lines 285-292, also 1613-1614):
I don’t think we should decide what to do with stellar rotation based on whether CHE options are on — I think it should
be based on whether the stars are actually in CHE. If stars are not in CHE, tidal options determine whether stellar rotations
should be updated.

Lines 285-292 are in code that is called in the constructor - at that stage we don't know if the stars are rotating fast enough to evolve as CH stars, we just know that CHE is enabled (in either mode) and that we need to check whether the stars are rotating fast enough to be evolving as CH stars. Is that not right?

BaseBinaryStar.cpp, lines 294—324: why is the treatment asymmetric for star1 vs. star2?

Do you mean this:

if (utils::Compare(m_Star1->Omega(), m_Star2->OmegaCHE()) >= 0) {

That would be a typo... and I think it has probably been there since CHE was implemented (unless I accidentally changed it in this PR, but it looks like I copy-pasted the code above it and missed changing that digit). But it shouldn't make any difference - m_Star1->Omega() and m_Star2->Omega() are the same at that point - really the same variable could be used in each conditional, but using the variable that pertains to the star makes it more apparent what the condition is checking. I will fix it though :-)

BaseBinaryStar.cpp, lines 1219-1223: I’d prefer using OrbitalAngularVelocity() and converting the output to appropriate
units rather than introducing another formula for calculating omega… oh, never mind, this is using previous values, OK.
Regardless, G_SN should really be renamed into a variable that describes the units of G rather than the function in which
it was first used. :-)

Yes, I believe that was code @reinhold-willcox added (including G_SN). Maybe @reinhold-willcox can suggest a better name, but if not I' can do it with this PR.

[jeffriley]

In Grid_demo.txt, why do we create the same binary twice? I think one binary should be sufficient for the detailed evolution
plot testing.

I don't know the answer to that - that's the way it was before this PR. EDIT: see Avi's answer above.

I see that the example bbh yaml has the mass loss prescription set to VINK rather than the new default FLEXIBLE2023

I thought I changed that to FLEXIBLE2023. (Just checked - I did in 9cd502d)

the WR multiplier set to 0.1, rather than the default of 1, and no settings at all for some of the other mass-loss choices.

I didn't change those - do we need to use the defaults?

It would be better to pick another set of initial conditions (masses, separations) rather than changing evolution settings defaults.

Ok - all I cared about here was getting the pytests to run - isn't the goal just to test that the functionality works, rather than produce a plot that is in some way useful? If that's not so, then sure, if you want to suggest what initial conditions and prescriptions should be used, I'll happily change it.

[jeffriley]

In conftest.py, we have
TEST_CONFIG_FNAME = os.path.join(TEST_CONFIG_DIR, "fiducial_bbh_config.yaml") but that's not the same as
example_bbh_compas_config.yaml . Is that intentional? If so, why do we need two different example bbh config
yaml's -- that seems like a recipe for future confusion? [Note: I haven't checked fiducial_bbh_config.yaml ]

I think we added one for the pytests - just to make sure we could test functionality. They don't need to be the same do they. But yes, maybe it is better if we combine them and use a single yaml file - but I don't think the fiducial yaml file allowed the pytests to run (I'm struggling to remember what we did...).

BTW, basic request: it would be nice to change one thing at a time (e.g., just tides in this case), makes it much harder to
review if many different changes appear at once.

Do you mean the inclusion of the code cleanups (if so, I addressed that above), or the inclusion of the changes to tests etc? If tests, we need to change some things to make the tests run... Or am I missing the point?

[jeffriley]

BaseBinaryStar.cpp, dealing with CHE stars (lines 285-292, also 1613-1614):

I forgot lines 1613-1614. Isn't that how we decided to implement CHE initially? (And partly why I asked about the interaction of CHE with the tides code). Prior to the implementation of CHE (I think - my memory may be failing me) we didn't set omega for the constituents of a binary - when we implemented CHE we chose to set omega for both stars to the same value (omega for the binary) - and then we decide if either (or both) is rotating fast enough to be evolving as a CH star. I haven't actually changed the functionality here - just rearranged it for performance).

[ilyamandel]

@jeffriley , @avivajpeyi , a few comments (and note my review is still not complete, a couple more items to come :) ) --

• Thank you for clarifying the distinction between updates to the old script used for the methods paper and the python testing script. Personally, I'd be in favour of minimising the number of such scripts and of keeping the test scripts as close to defaults as possible (we really want to make sure that default options produce sensible outcomes -- they are the recommended defaults, after all -- and because folks might be tempted to use these as examples). But I leave it to you to decide what's easiest to do in this PR.
• I take your point about take advantages of the opportunity to improve the code (and apologise for the use of "prettifying" -- it was not meant dismissively). I agree that this is more important than review convenience, so if the options are either to make those improvements simultaneously with new development or not at all, I definitely vote for making improvements when the opportunity arises.
• CHE and tides: my logic is that if someone turns on the CHE option, they are letting us know that they meant that tides should ensure spin-up and CHE if we actually satisfy the CHE criteria, but not necessarily if we don't (i.e., if the CHE option is turned on but stars don't satisfy the CHE threshold, and tides are off, there's no co-rotation). Setting omega initially: I was thinking of the possibility that, in the future, we may wish to assign some meaningful initial spins to stars and, if tides are off and CHE is not happening, we shouldn't over-ride these. BTW, see also lines 2229-2230.

[ilyamandel]

why is the treatment asymmetric for star1 vs. star2?

I meant that star1 has three if-else clauses on lines 304-313 of BaseBinaryStar.cpp, but star2 has only one clause on lines 322-324.

G_SN:

constexpr double G = 6.67E-11; // Gravitational constant in m^3 kg^-1 s^-2 (more accurately known as G M_sol)
constexpr double G_CGS = 6.6743E-8; // Gravitational constant in cm^3 g^-1 s^-2
constexpr double G1 = 4.0 * M_PI * M_PI; // Gravitational constant in AU^3 Msol^-1 yr^-2
constexpr double G_SN = G * 1.0E-9 / KG_TO_MSOL; // Gravitational constant in km^3 Msol^-1 s^-2, for use in the ResolveSupernova() function
constexpr double G_SOLAR_YEAR = 3.14E7; // Gravitational constant in Lsol Rsol yr Msol^-2 for calculating photon tiring limit

I'd suggest renaming G1 and G_SN in a way that describes their unit systems (like G_CGS and G_SOLAR_YEAR), e.g., G1 -> G_AUMSOLYR, G_SN->G_KMMSOLS, and removing the comment "for use in the ResolveSupernova() function" since it can be (and now is) used elsewhere.

[ilyamandel]

Continuing review (still not complete):

5. BaseBinaryStar.cpp, lines 2327: As mentioned previously, I think it's dangerous to do nothing if a root is not found. This should generally mean a Darwin-instability driven common envelope event.

6. BaseBinaryStar.cpp, line 2332: Since you don’t like non-integer powers (yes, I understand why :) ), could simplify this expression to replace three with one:
   m_SemiMajorAxis = std::cbrt(G1 * (m_Star1->Mass() + m_Star2->Mass()) / m_Omega / m_Omega);

[ilyamandel]

7. BaseBinaryStar.h, line 681:
   if (it < maxit) { // too many iterations?
   I don’t think that’s the right comment here (unlike on line 686)

8. It looks like m_AngularMomentum is now defined in BaseStar and set to the default value in BaseStar.cpp, line 144, but does not have a meaningful (physical) value until the star becomes an NS. Do we really need this? I think it’s risky to keep around a variable that’s incorrect for most stellar types: some day someone may output it and will get meaningless values…

9. utils.cpp::Compare() [lines 193—217]
   I think if(p_Tolerance > 0.0) can be taken outside the check for #ifdef COMPARE_GLOBAL_TOLERANCE — no need to repeat that code block twice. BTW, in expressions such as std::max(std::abs(p_X), fabs(p_Y)), why is there an asymmetry in the treatment of p_X and p_Y?

10. utils.cpp, lines 235–239:
    This has to be one of the most round-about ways to do the calculation. :)
    How about just the following one-liner?
    return std::cbrt((p_Mass1 + m_Mass2) * p_Period * p_Period / 365.25 / 365.25)
    [If you don’t like 365.25, you can define DAYS_IN_YEAR as 365.25, and then change the definition of
    SECONDS_IN_DAY = SECONDS_IN_YEAR / DAYS_IN_YEAR in constants.h]

11. constants.h, line 237:
    constexpr double _2_PI = M_PI * 2; // 2PI
    I am surprised you didn’t change that to 2.0. :)
    BTW, what’s the logic behind using M in constant names? E.g., why _2_PI but SQRT_M_2_PI?

12. I see STAR_PROPERTY::MOMENT_OF_INERTIA now listed in constants.h, but I don’t see it defined in BaseStar.cpp — is GitHub not showing that change to me?

[ilyamandel]

OK, I think this is the last batch. :)

13. NS.cpp, line 288:
    double r_cm = m_Radius * KM_TO_CM;
    I think this is a mistake — if you want to convert the radius from km into cm, this should be
    double r_cm = radius * KM_TO_CM;

14. NS.cpp, line 308:

• @param [IN] p_Radius Radius of the Neutron Star in metres
  I presume that should say in km, given that we later convert it from km into cm; otherwise, the conversion is incorrect!

15. NS.cpp, line 359:
    m_AngularMomentum = _2_PI * m_MomentOfInertia_CGS / (m_PulsarDetails.spinPeriod * SECONDS_IN_MS) * factor;
    We’ve already computed the spin frequency, so can simplify this to
    m_AngularMomentum = m_MomentOfInertia_CGS * factor * m_PulsarDetails.spinFrequency
    (But see below!)

16. NS.cpp, line 417:
    m_AngularMomentum = m_PulsarDetails.spinFrequency * m_MomentOfInertia_CGS; // angular momentum of star
    This will clearly have different units from the previous calculation of the angular momentum (cf. previous comment) — we should be consistent in our units! [I think we assume elsewhere that the angular momentum is in CGS units — if so, we don’t need the factor in item 15]

17. NS.cpp, line 460:
    // his part of the code
    ->
    // This part of the code
    :)

[reinhold-willcox]

Point taken about G_SN. I'll fix that, and G1, now.

[jeffriley]

@ilyamandel I think I have addressed all your review points - I may not have fixed them all yet, but I have addressed them. Notes below.

1. You're right - it does make reviewing difficult. I think what I'll do in future is note/record changes I want to make, and push them as a separate PR.

2. See my response in the PR - re timings. Thoughts?

3. Per our discussion on slack, this is withdrawn?

4. G1 and G_SN renamed by @Reinhold (see his commits).

5. Done - but you should check what I did :-)

6. Very nice - should have seen that myself. Done.

7. Well, it kind-of is in a negative way :-) (especially if you read the comment on the next line). I changed the comment.

8. My fault. I moved it into BaseStar (for some reason I thought we'd need it). I have now moved it back to NS and renamed it m_AngularMomentum_CGS (see point 16 below).

9a. > I think if(p_Tolerance > 0.0) can be taken outside the check ...

I suppose it could, but then we'd need to add more #ifdef and #else statement to cover the else part of the if (p_Tolerance > 0.0) conditionals. I think it's cleaner/more straightfoward the way it is.

EDIT: Oh wait - I see what you mean. Yes, I'll take another look.

9b. > BTW, in expressions such as std::max(std::abs(p_X) ...

I think probably distraction and forgetfulness :-) fabs (an overload introduced in C++11) is probably marginally safer, so I converted the remaining std::abs calls to fabs (in utils::Compare())

10. Fixed as suggested (I don't know who wrote that function - could have been me :-))

11. I missed it... now fixed. M_2_PI is a constant defined in the cmath header. I added _2_PI and SQRT_M_2_PI (see constants.h lines 169..178).

12. Missed it. Now fixed (added to BaseStar::StellarPropertyValue()).

13. Fixed (slightly different solution).

14. Fixed (km is correct - called from NS::CalculateAndSetPulsarParameters() with radius in km).

15. Fixed as suggested.

16. Moved class variable m_AngularMomentum from BaseStar to NS (see point 8 above) and renamed it m_AngularMomentum_CGS. Made allcalculation in CGS.

17. Fixed (also fixed typo on the same line).

(I had to play around a bit with git to get Reinhold's commit for the G constants - I don't think I broke anything or lost any changes - at least I hope not).

[jeffriley]

@ilyamandel Last change pushed - ready for you to re-review (at your convenience).

[ilyamandel]

Hi @jeffriley ,
I am very sorry for taking so long to review the latest changes.
I am happy to approve this PR. I leave it to you whether to address the following (none are technically incorrect or mission-critical, though the first item may merit returning to):

---

BaseBinaryStar.cpp, dealing with CHE stars (lines 285-292, also 1613-1614):
I don’t think we should decide what to do with stellar rotation based on whether CHE options are on — I think it should be based on whether the stars are actually in CHE. If stars are not in CHE, tidal options determine whether stellar rotations should be updated.
(And, from another comment on this:
CHE and tides: my logic is that if someone turns on the CHE option, they are letting us know that they meant that tides should ensure spin-up and CHE if we actually satisfy the CHE criteria, but not necessarily if we don't (i.e., if the CHE option is turned on but stars don't satisfy the CHE threshold, and tides are off, there's no co-rotation). Setting omega initially: I was thinking of the possibility that, in the future, we may wish to assign some meaningful initial spins to stars and, if tides are off and CHE is not happening, we shouldn't over-ride these. BTW, see also lines 2229-2230.)

---

Calling OrbitalAngularVelocity() instead of introducing m_Omega -- OK, let's leave things as you've implemented them. :-)

---

utils.cpp::Compare() [lines 193—217]
I think if(p_Tolerance > 0.0) can be taken outside the check for #ifdef COMPARE_GLOBAL_TOLERANCE — no need to repeat that code block twice.

I am not sure if you addressed this one, but it's not a big deal either way.

[jeffriley]

Thanks @ilyamandel

BaseBinaryStar.cpp, dealing with CHE stars (lines 285-292, also 1613-1614):

We may need to discuss so I better understand what you want here - also see my comment here:

#1012 (comment)

utils.cpp::Compare() [lines 193—217]
I am not sure if you addressed this one, but it's not a big deal either way.

I did. At least I think I did - I pulled that check out of the #ifdef as you suggested,

[ilyamandel]

Ahh, right -- sorry! Ok, I suggest pulling the PR at this stage. Thank you once again!

…

On Mon, 27 Nov 2023, 11:56 pm jeffriley, ***@***.***> wrote: Thanks @ilyamandel <https://github.com/ilyamandel> BaseBinaryStar.cpp, dealing with CHE stars (lines 285-292, also 1613-1614): We may need to discuss so I better understand what you want here - also see my comment here: #1012 (comment) <#1012 (comment)> utils.cpp::Compare() [lines 193—217] I am not sure if you addressed this one, but it's not a big deal either way. I did. At least I think I did - I pulled that check out of the #ifdef as you suggested, — Reply to this email directly, view it on GitHub <#1012 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAD5BVJ7MKSUZ64TODZT7ZLYGULDVAVCNFSM6AAAAAA6WNX7ACVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQMRYG43DKMRVGM> . You are receiving this because you were mentioned.Message ID: ***@***.***>