How to set i_rot when omega is zero?

[rjfarmer]

From Warwick due to changes made in commit 12599:

This is one of the errors I wanted to bring up in more detail. The problematic line (before this commit) was

 s% i_rot(k) = s% j_rot(k)/s% omega(k)

in hydro_rotation. If the model had zero rotation when we reach this point, we divide by zero.

I noticed, however, that the call to this subroutine on line 125 of read_model is preceded by the comment:

 ! older MESA versions stored only omega in saved models. However, when
 ! using rotation dependent moments of inertia one actually needs to store
 ! the angular momentum in order to initialize the model. This flag is here
 ! to account for the loading of old saved models.
 ! TODO: discuss inconsistencies

so I think we're being forced to discuss the inconsistencies.

Fundamentally, the issue is trying to recover the moment of inertia from the angular momentum and rotation rate when the rotation rate is zero. Rather than just set the moment of inertia to zero, I think we should either raise an error or use one of the other routines to recompute the moment of inertia from the structure. e.g. set_i_rot in hydro_rotation. But I defer to one of the rotation experts.

[orlox]

This should be an easy fix. when omega is zero the moment of inertia is just that of a shell. I'll try and get a fix tomorrow.

[aarondotter]

When s% rotation_flag = .false. I use the following to calculate s% i_rot:

s% i_rot = 0.0d0
do i=s% nz, 2, -1
    s% i_rot(k) = 0.4d0*(pow5(s% r(i)) - pow5(s% r(i-1)) )/(pow3(s% r(i)) - pow3(s% r(i-1)))   
enddo

[evbauer]

Just curious, why not just... s% i_rot(k) = (2d0/3d0)*pow2(s% r(k))

[orlox]

Evan, this is because the deformation causes a digression from being simply spherical shells (figure 33 in MESA V). In the limit of critical rotation we actually have i_rot = r_e^2.

What you compute Aaron is the moment of inertia of a finite sized shell. This only leads to very small differences compared to using the thin shell approximation, and we don't have the equivalent expression that accounts for the deformation from rotation which is a much larger effect.

On Sat, Jan 18, 2020, 1:39 AM Evan B. Bauer notifications@github.com wrote:

Just curious, why not just... s% i_rot(k) = (2d0/3d0)*pow2(s% r(k))

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/MESAHub/MESA/issues/1?email_source=notifications&email_token=AAEDFIBV5O5QYVAF62TVTA3Q6JFVDA5CNFSM4KHW3K52YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEJJLQFY#issuecomment-575846423, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAEDFIG3ODFONC2RYUBXSPLQ6JFVDANCNFSM4KHW3K5Q .

[aarondotter]

Pablo, I don't understand why you mention deformation due to rotation because we're discussing strictly the non-rotating case here, i.e., when s% rotation_flag = .false.. Please enlighten me!

I can imagine that the formula I use (for finite shell) could lead to some numerical issues when taking R^3 or R^5 and subtracting but I haven't experienced this yet. This would not be present in the thin shell approximation. I opted for the finite version over the thin shell because I thought it would be better not to assume dR << R in all cases.

In any case, we should be able to choose from thin-shell or finite-shell for the non-rotating case when we need to populate s% i_rot. Just need to decide which one.

[rjfarmer]

Technically we have rotation_flag true (this came up in relation to the check_redo test case) but the new omega was 0. When it tries to do a restart its loading from data with omega=0

[orlox]

As Rob says, we're discussing what to do in cases when the rotation flag is on. When rotation is not set, we do not even evaluate the moment of inertia as its an unnecessary overhead.

On Sat, Jan 18, 2020 at 3:47 PM Robert Farmer notifications@github.com wrote:

Technically we have rotation_flag true (this came up in relation to the check_redo test case) but the new omega was 0. When it tries to do a restart its loading from data with omega=0

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/MESAHub/MESA/issues/1?email_source=notifications&email_token=AAEDFICNFD66OXDFQ5AMLU3Q6MI7DA5CNFSM4KHW3K52YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEJJZ3BA#issuecomment-575905156, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAEDFIB7AW2EPEH32PEMTLLQ6MI7DANCNFSM4KHW3K5Q .

-- Pablo Marchant Campos M.Sc on Astrophysics, Universidad Católica de Chile PhD on Astrophysics, Argelander-Institut für Astronomie, Universität Bonn

[evbauer]

I've gotten lost on the goal here. I thought the only thing that needed correction was the moment of inertia in the case when s% omega(k) == 0. In that case, there isn't any deformation, and the moment of inertia is perfectly well defined (with some minor choice to be made regarding how precise to be about finite shell thickness). Isn't that all that needs to be implemented?

Of course the deformation should be accounted for with non-zero rotation, but I thought that case was already handled properly.

[orlox]

That's right Evan, sorry if I misunderstood. Thought you meant using 2/3r^2 always, while the single thing needed is to setup the exception for omega=0.

Need to make the fix, have not had the time but will put it up tomorrow. In any case, the issue tracker is a pretty good way to keep a record of things still not fixed despite earlier promises to do so shortly :).

On Mon, Jan 20, 2020, 10:34 PM Evan B. Bauer notifications@github.com wrote:

I've gotten lost on the goal here. I thought the only thing that needed correction was the moment of inertia in the case when s% omega(k) == 0. In that case, there isn't any deformation, and the moment of inertia is perfectly well defined (with some minor choice to be made regarding how precise to be about finite shell thickness). Isn't that all that needs to be implemented?

Of course the deformation should be accounted for with non-zero rotation, but I thought that case was already handled properly.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/MESAHub/MESA/issues/1?email_source=notifications&email_token=AAEDFIHU22MJDAX7EMASO73Q6YKGPA5CNFSM4KHW3K52YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEJN3USY#issuecomment-576436811, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAEDFIFVNOQTQBSIXASE4S3Q6YKGPANCNFSM4KHW3K5Q .

[evbauer]

Ah, thanks. That makes sense. I had really just meant to ask why we care about finite shell thickness for the moment of inertia in a simple case like this, but that was rather tangential anyway, so I should have been more clear what I was getting at.

[jschwab]

@orlox Has this been resolved?

[rjfarmer]

Anyone want to push this forward @orlox @evbauer @aarondotter before we release? Otherwise I'll revert the changes i did and we ship with something broken but at least not wrong

[aarondotter]

I'll do it later today unless someone else says they want to do it first.

My understanding is that when omega = 0 we will use the thin shell approximation to set i_rot.

Aaron

On Fri, Feb 21, 2020 at 8:10 AM Robert Farmer notifications@github.com wrote:

Anyone want to push this forward @orlox https://github.com/orlox @evbauer https://github.com/evbauer @aarondotter https://github.com/aarondotter before we release? Otherwise I'll revert the changes i did and we ship with something broken but at least not wrong

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/MESAHub/MESA/issues/1?email_source=notifications&email_token=AAYW4L4YEUWB44EIV5DSGMLRD7HDTA5CNFSM4KHW3K52YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEMSUWSA#issuecomment-589646664, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAYW4LYBTLAP674WDP53ZRDRD7HDTANCNFSM4KHW3K5Q .

[orlox]

I'll get this done today. Sorry for the delay.

On Fri, Feb 21, 2020 at 2:24 PM Aaron Dotter notifications@github.com wrote:

I'll do it later today unless someone else says they want to do it first.

My understanding is that when omega = 0 we will use the thin shell approximation to set i_rot.

Aaron

On Fri, Feb 21, 2020 at 8:10 AM Robert Farmer notifications@github.com wrote:

Anyone want to push this forward @orlox https://github.com/orlox @evbauer https://github.com/evbauer @aarondotter https://github.com/aarondotter before we release? Otherwise I'll revert the changes i did and we ship with something broken but at least not wrong

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub < https://github.com/MESAHub/MESA/issues/1?email_source=notifications&email_token=AAYW4L4YEUWB44EIV5DSGMLRD7HDTA5CNFSM4KHW3K52YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEMSUWSA#issuecomment-589646664 , or unsubscribe < https://github.com/notifications/unsubscribe-auth/AAYW4LYBTLAP674WDP53ZRDRD7HDTANCNFSM4KHW3K5Q

.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/MESAHub/MESA/issues/1?email_source=notifications&email_token=AAEDFIAP53VVXUTKSG7JRWDRD7IZDA5CNFSM4KHW3K52YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEMSV4AY#issuecomment-589651459, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAEDFIB3AU7D5IT2WMPOQLDRD7IZDANCNFSM4KHW3K5Q .

-- Pablo Marchant Campos M.Sc on Astrophysics, Universidad Católica de Chile PhD on Astrophysics, Argelander-Institut für Astronomie, Universität Bonn

[orlox]

Should be fixed now, made some cleanup of TODO's I left for myself as well. Most of those were deficient explanations of what was actually being done, but MESA V contains those details so we can refer to it.