Page 1 of 1

Possible errors in rotation kernels [2.2]

Posted: Tue Oct 15, 2013 9:36 am
by 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

Image

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

Image

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 Image 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.

Re: Possible errors in rotation kernels [2.2]

Posted: Tue Oct 15, 2013 9:56 am
by 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!

Re: Possible errors in rotation kernels [2.2]

Posted: Tue Oct 15, 2013 10:00 am
by rhtownsend
warrick wrote: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...

Re: Possible errors in rotation kernels [2.2]

Posted: Thu Oct 17, 2013 3:45 am
by rhtownsend
OK, this issue is now fixed in the development version. I'll be making a new release soon.