Issue 1034 - angular momentum not conserved in tides code

[jeffriley]

Fix for issue #1034

• This fix changes the functions
  . BaseBinaryStar::CalculateAngularMomentum(),
  . BaseBinaryStar::CalculateTotalEnergy(), and
  . BaseStar::AngularMomentum()
  to use moment of inertia rather than gyration radius. This has wider implications than just issue --enable-tides Angular Momentum is not conserved #1034 and may change DCO yields slightly.
• Removed some unused functions.
• Change to functionality (noted above) noted in 'What's New' online documentation page

@ilyamandel ready for review
@veome22 can you run your tests and check this does what you expect?

[ilyamandel]

I think there may be a problem with units.

Line 2074 of BaseBinaryStar.cpp seems to be adding quantities in which distances are measured in RSOL to quantities in which distances are measured in AU. There is a function CalculateMomentOfInertiaAU which would solve the problem if this was used on lines 275-276.

[Aside: might not be a bad idea to note the expected units of inputs and outputs in complex functions like CalculateTotalEnergy() ]

[jeffriley]

There is a function CalculateMomentOfInertiaAU

Ahhh - that's what I meant to use! Thanks @ilyamandel - I'll change that.

[ilyamandel]

Same unit issue in BaseStar::OrbitalAngularMomentum() will be addressed by the change above.

[ilyamandel]

Same change requested in

BaseStar.h,
doubleAngularMomentum() const { return CalculateMomentOfInertia() * m_Omega; }
-- should use CalculateMomentOfInertiaAU()

[ilyamandel]

I suggest we remove CalculateGyrationRadius() functions altogether -- they are currently inconsistent with the MomentOfInertia() calculations, and will only confuse the unwary user.

[jeffriley]

I suggest we remove CalculateGyrationRadius() functions altogether ...

It is used in BinaryConstituentStar::CalculateSynchronisationTimescale() - can we rewrite that?

[ilyamandel]

It is used in BinaryConstituentStar::CalculateSynchronisationTimescale() - can we rewrite that?

Yes. The gyration radius [squared] is just CalculateMomentOfInertia()/Mass()/Radius()/Radius() -- see text between Eq. 26 and 27 of Hurley, Pols, Tout (2002).

[jeffriley]

@ilyamandel

I suggest we remove CalculateGyrationRadius() functions altogether ...

Done.

CalculateGyrationRadius() was also used in MainSequence::CalculateMomentOfInertia(), so I moved GiantBranch::MomentOfInertia to the BaseStar class and removed both MainSequence::CalculateMomentOfInertia() and GiantBranch::MomentOfInertia (and fixed calls as appropriate). Hurely et al. 2000 seems to suggest this is correct (see paragraph following eq 109) - let me know if I've misinterpreted that.

Please look at what I did in BinaryConstituentStar::CalculateSynchronisationTimescale() - I think, based on the name of the variable, that there was an existing defect there that has gone unnoticed. Or, the variable could have been badly named... I've changed the code assuming the variable name was meant to indicated what should happen in the code - let me know if I need to change that.

[jeffriley]

@veome22 I'm not convinced these changes have completely addressed the original problem - let me know if, after running your tests, there's more work to do.

[ilyamandel]

BaseStar.h, line 190:

virtual double CalculateMomentOfInertiaAU(const double p_RemnantRadius = 0.0) const { return CalculateMomentOfInertia(p_RemnantRadius * RSOL_TO_AU) * RSOL_TO_AU * RSOL_TO_AU; }

Why are you converting p_RemnantRadius to AU? You need to convert the output to AU (which is already done), not the input...

[ilyamandel]

I like the general approach to calculating the moment of inertia -- but, at the moment, it seems to be called without an argument when calculating the angular momentum, whereas it should either have the core radius passed in (as written, via p_RemnantRadius argument) or should determine its own core radius (OOP). Otherwise, the (0.21 * m_CoreMass * p_RemnantRadius * p_RemnantRadius) term never contributes anything. Did I miss something?

[jeffriley]

Why are you converting p_RemnantRadius to AU?
Otherwise, the (0.21 * m_CoreMass * p_RemnantRadius * p_RemnantRadius) term never contributes anything. Did I miss something?

I don't think you did :-) I think the p_RemnantRadius is a hangover from the original refactor - it looks like we no longer use it. I'll take another look (I did look, but didn't want to deal with it now :-)) As far as why it is converted to AU - I can only think that in the original code it was passed in Rsol...

I'll look at it later today.

[jeffriley]

Fixed a few places that I had missed changing CalculateMomentOfInertia() to CalculateMomentOfInertiaAU().

I believe this now resolves the issue of angular momentum not being conserved when tides are enabled - @veome22 could you confirm?

[jeffriley]

@ilyamandel I have changed the implementation of CalculateMomentOfInertia() to properly implement Hurley et al., 2000 eq 109 - please check my changes. This version still resolves the issue of angular momentum not being conserved when tides are enabled.

EDIT: I have not checked if/how the changes made in this PR affect DCO yields.

EDIT: I probably should pull the implementation of GiantBranch::CalculateMomentOfInertia() out of the header and into the cpp file. If you're otherwise ok with the implementation I'll push a fix for that.

[ilyamandel]

BinaryConstituentStar.h, line 387:

double gyrationRadiusSquared = CalculateMomentOfInertiaAU() / Mass() / Radius() / Radius();

Isn't Radius() returning results in Rsun? Then this should use CalculateMomentOfInertia(), not CalculateMomentOfInertiaAU()

[ilyamandel]

BH.h, line 58 [and similarly in MR.h, NS.h, MS.h]:

double CalculateMomentOfInertiaAU() const { return CalculateMomentOfInertia() * RSOL_TO_AU * RSOL_TO_AU; }

Why do we need this? Isn't CalculateMomentOfInertiaAU() always just
{ return CalculateMomentOfInertia() * RSOL_TO_AU * RSOL_TO_AU; }
?
If so, we don't need to over-write it -- just define that once at the top of the inheritance hierarchy, it should automatically pick up the relevant CalculateMomentOfInertia(), right?

[ilyamandel]

Other than the two comments above (the first needs to be fixed, the second is optional), looks good, happy to approve -- thank you, @jeffriley !

[veome22]

Fixed a few places that I had missed changing CalculateMomentOfInertia() to CalculateMomentOfInertiaAU().

I believe this now resolves the issue of angular momentum not being conserved when tides are enabled - @veome22 could you confirm?

@jeffriley I just ran the same binary as before, and angular momentum is indeed conserved now! I am running a larger test now and will keep you posted.

[jeffriley]

double gyrationRadiusSquared = CalculateMomentOfInertiaAU() / Mass() / Radius() / Radius();

Ah - yes, I went too far :-) I'll fix that.

If so, we don't need to over-write it -- just define that once at the top of the inheritance hierarchy, it should automatically pick up the relevant CalculateMomentOfInertia(), right?

Yes, it should - I'll take that out. Good catch.

[jeffriley]

@jeffriley I just ran the same binary as before, and angular momentum is indeed conserved now! I am running a larger test now and will keep you posted.

Excellent - thanks @veome22

[jeffriley]

@ilyamandel Requested changes made. Thanks for the review and picking up my missteps :-)

[ilyamandel]

Looks good -- thanks again, @jeffriley !

Will leave it for you to merge.

[veome22]

Just completed the final check. Angular momentum is now conserved in all binaries!
Thank you @jeffriley :)