Possible errors in rotation kernels [2.2]

[warrick]

Hi,

I've been using GYRE (adiabatically) to compute some rotational kernels but something is going wrong. I've attached an input FGONG model (produced by MESA's solar calibration test case), my GYRE inlist, and output for one mode (n=4, l=2). Here's also a plot of the kernel for that mode as produced by GYRE (red) and my own evaluation of the formula in Aerts, Christensen-Dalsgaard and Kurtz's book (equation 8.42), with the normalization factor fitted roughly by eye. The formula I use, to be explicit, is



with the radial profiles read from the file I've attached (mode-0061.txt).



As you can see, they're different. In particular, the differences seem to begin when the radial displacement eigenfunction changes sign. i.e. the region where the two curves differ is bounded on the left by a point where changes sign. I hope that's helpful.

As usual, I provide the usual disclaimer that I might be doing something silly, but the kernel looks a bit fishy anyway...

Cheers,
Warrick

PS: I forgot to mention, this is GYRE 2.2, compiled with MESA SDK (8 April 2013 version), running on Scientific Linux 6.4 on 4 OpenMP cores.

[warrick]

Okay, so I dug into the source code* and found the probable error. Line 826 of src/ad/gyre_mode.f90 is this:

Code: Select all

K = (ABS(xi_r)**2 + (l*(l+1)-1)*ABS(xi_h)**2 - 2._WP*ABS(xi_r*xi_h))*U*x**2/c_1

I'm pretty sure the product "xi_r*xi_h" (the last term in the brackets) shouldn't be wrapped in an absolute value. If I modify the formula I used in the previous plot, it matches the GYRE curve pretty much perfectly. (Well, except for the normalization, which was by eye anyway.)

I'll presume that's the error and do a temporary fix myself.

Cheers,
Warrick

*...which I must say, is quite beautiful. Three greps of increasing intelligence, and the function was dead easy to find!

[rhtownsend]

Okay, so I dug into the source code* and found the probable error. Line 826 of src/ad/gyre_mode.f90 is this:

Code: Select all

K = (ABS(xi_r)**2 + (l*(l+1)-1)*ABS(xi_h)**2 - 2._WP*ABS(xi_r*xi_h))*U*x**2/c_1

I'm pretty sure the product "xi_r*xi_h" (the last term in the brackets) shouldn't be wrapped in an absolute value. If I modify the formula I used in the previous plot, it matches the GYRE curve pretty much perfectly. (Well, except for the normalization, which was by eye anyway.)

I'll presume that's the error and do a temporary fix myself.

Cheers,
Warrick

Sounds legit; I'll roll the fix into the next update. Stuff like this comes from mapping adiabatic expressions (which always have real eigenfunctions) over to non-adiabatic ones (which generally don't); I should probably be using real() instead of abs().

*...which I must say, is quite beautiful. Three greps of increasing intelligence, and the function was dead easy to find!

Why, thank you! It's not often that someone opens up the bonnet and looks at the engine...

[rhtownsend]

OK, this issue is now fixed in the development version. I'll be making a new release soon.