Normalized inertias: GYRE vs ADIPLS

[warrick]

Hi,

I've been mucking around with some results that use mode inertias. I found myself struggling to recover some previous results I computed with ADIPLS, and I seem to have narrowed the issue down to the mode inertias. Now, before I get too far ahead, I'm not sure this is really a bug rather than me not being careful about what is being computed, but the outputs take the same name and I'm not sure where the difference is.

Here's a plot of the normalized inertias for ADIPLS (blue) vs those of GYRE (green) for Model S from the ADIPLS distribution nested in MESA (mesa/adipls/adipls_tar_files/models.c.tar.gz/fgong.l5bi.d.15). The input files for GYRE and ADIPLS are attached. (Calculations done with GYRE 2.3.)

Plot of normalized mode intertias for Model S as computed by ADIPLS (blue) and GYRE (green)

inertia_adipls_gyre_modelS.png (51.58 KiB) Viewed 6064 times

I'm not worried about the overall difference but the shapes of the curves are obviously different, especially at high frequencies. The ADIPLS curve flattens out (and actually turns around) whereas the GYRE curve keeps decreasing.

I've dug around quite a lot into the definitions of the normalizations because I know the linear perturbation analysis only defines the solution up to a constant scaling. I've used ekinr=1 in ADIPLS (see eqn 4.3b in the ADIPLS notes), which should correspond to inertia_norm_type = 'BOTH' in GYRE. I've dug around in the (beautifully clear!) GYRE code and it looks like the definition of E there is fine, up to a constant.

If I look at the unnormalized inertias (i.e. 'E' in the output), I find that they are always basically one. So I'm inferring that the normalization condition is such that this is true. For ADIPLS, I'm not so sure what's happening but it's implied that the condition is y_1 = xi_r / R = 1. (BTW, can phpBB do maths?) This could easily be the cause of the problem, but I'm not actually sure if this is the condition he uses and, in that case, which is right...

I would check the eigenfunctions directly but I'm not yet well-versed enough in ADIPLS to extract that kind of information. I'll get cracking, though.

W

[rhtownsend]

Hi Warrick --

This does indeed look strange. I've reproduced your GYRE results using the latest version, and I'm now looking into what ADIPLS predicts. Thanks for bringing this to my attention.

cheers,

Rich

[rhtownsend]

Hi Warrick --

Can you post the 'agsm' input file being used for the frequency search in your adipls case?

cheers,

Rich

[warrick]

Here it is. Let me know if I can be of any help.

[rhtownsend]

Here it is. Let me know if I can be of any help.

Hmm, I can't see it -- can you repost, or mail it to me directly?

[warrick]

Ah, I missed the little warning that says "The extension agsm is not allowed." I'll email it directly to you.

[warrick]

Hmmm... I've been looking at various definitions and normalizations for the eigenvectors, but then I realised that they're irrelevant, since the definition of the normalized inertia is invariant with a rescaling of the eigenvectors. (Also, for the two that I checked, the results are indistinguishable.) But, for what it's worth, I think the missing factor is (4*pi)**2. In particular, when I expand the definition of dE_dx in gyre/src/ad/gyre_mode.f90, it looks like the leading 4*pi should be in the denominator, not the numerator. Here's a calculation that hopefully makes some sense:

CodeCogsEqn.gif (2.47 KiB) Viewed 6048 times

This times xi_r^2+l(l+1)xi_h^2 is what is integrated in GYRE, and then normalized by dividing by the surface value of xi_r^2+l(l+1)xi_h^2. So the M_* is in dE/dx to start off with.

So I think the equations are correctly implemented. It also doesn't seem to be related to the surface boundary condition. That introduces some differences, but not much and certainly not enough to close the gap. Could there be some other surface-related issue, like how much of the atmosphere is included?

[rhtownsend]

Hmmm... I've been looking at various definitions and normalizations for the eigenvectors, but then I realised that they're irrelevant, since the definition of the normalized inertia is invariant with a rescaling of the eigenvectors. (Also, for the two that I checked, the results are indistinguishable.) But, for what it's worth, I think the missing factor is (4*pi)**2. In particular, when I expand the definition of dE_dx in gyre/src/ad/gyre_mode.f90, it looks like the leading 4*pi should be in the denominator, not the numerator. Here's a calculation that hopefully makes some sense:

CodeCogsEqn.gif

This times xi_r^2+l(l+1)xi_h^2 is what is integrated in GYRE, and then normalized by dividing by the surface value of xi_r^2+l(l+1)xi_h^2. So the M_* is in dE/dx to start off with.

So I think the equations are correctly implemented. It also doesn't seem to be related to the surface boundary condition. That introduces some differences, but not much and certainly not enough to close the gap. Could there be some other surface-related issue, like how much of the atmosphere is included?

Thanks for spotting the 4pi issue. But it will only give a factor of 4pi -- rather than going from the numerator to the denominator in my expression, it should just be removed from the numerator I think.

Regarding the mismatch between the codes, I wonder whether things will get better if you switch both to a different normalization condition for the inertia?

[warrick]

Good idea. When I set both to the radial normalization (iekinr = 0 and RADIAL), the agreement is fantastic. In fact, the issue now is that I get the same normalized intertiae whether I use RADIAL or BOTH in GYRE. I've tried looking at the source to see if there's something up but it doesn't look like it (to my relatively-untrained eye, anyway).

While we're at it, though, I'll try patching to the latest version (I'm on 2.3) and see if anything changes.

[rhtownsend]

Good idea. When I set both to the radial normalization (iekinr = 0 and RADIAL), the agreement is fantastic. In fact, the issue now is that I get the same normalized intertiae whether I use RADIAL or BOTH in GYRE. I've tried looking at the source to see if there's something up but it doesn't look like it (to my relatively-untrained eye, anyway).

While we're at it, though, I'll try patching to the latest version (I'm on 2.3) and see if anything changes.

For l=0 modes the two inertia normalization methods are identical because xi_h = 0 always.

For l > 0 modes, the two methods aren't the same; but in this specific case, because we're looking at p-modes, xi_h is small and so the two methods give pretty close results (esp. at highest frequencies).

I've checked the GYRE data, and this seems to explain what we're seeing: the l=0 inertias are identical, but the l>0 inertias are similar but not identical.

Regarding the mismatch between ADIPLS and GYRE in the ieknir=1/BOTH case, I think this is because I set R_phot = R_s in GYRE. I'll look into whether there's a better way to get R_phot.

cheers,

Rich