Floating point exception while searching for l=1 g modes

[earlbellinger]

I am searching for g-modes in a model with a helium burning core and a hydrogen burning shell.

It can do the very low frequency modes (n<<-1000) but once it reaches a radial order of -189 it crashes with the following FPE:

...
At line 2164 of file gyre_wave.f90
Fortran runtime warning: An array temporary was created
   1    0    -191      0    191  0.13474501E+01  0.00000000E+00 0.6408E-11      7
At line 2164 of file gyre_wave.f90
Fortran runtime warning: An array temporary was created
   1    0    -190      0    190  0.13533290E+01  0.00000000E+00 0.6627E-11      7
At line 2164 of file gyre_wave.f90
Fortran runtime warning: An array temporary was created
   1    0    -189      0    189  0.13622276E+01  0.00000000E+00 0.4887E-11      8

Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:
#0  0x7fa89405a4ff in ???
        at /tmp/S/glibc-2.37/signal/../sysdeps/unix/sysv/linux/x86_64/libc_sigaction.c:0
#1  0x6bbf19 in __gyre_r_block_sysmtx_MOD_scale_rows_
        at /afs/mpa/data/earl/gyre/gyre-6.0.1/src/build/gyre_r_block_sysmtx.f90:462
#2  0x6bc973 in __gyre_r_block_sysmtx_MOD_factor
        at /afs/mpa/data/earl/gyre/gyre-6.0.1/src/build/gyre_r_block_sysmtx.f90:339
#3  0x7fba77 in __gyre_r_bvp_MOD_factor
        at /afs/mpa/data/earl/gyre/gyre-6.0.1/src/build/gyre_r_bvp.f90:232
#4  0x568fdc in __gyre_r_discrim_func_MOD_eval_
        at /afs/mpa/data/earl/gyre/gyre-6.0.1/src/build/gyre_r_discrim_func.f90:140
#5  0x646f3d in narrow_brent_
        at /afs/mpa/data/earl/gyre/gyre-6.0.1/src/build/gyre_r_root.f90:369
#6  0x647081 in __gyre_r_root_MOD_narrow_
        at /afs/mpa/data/earl/gyre/gyre-6.0.1/src/build/gyre_r_root.f90:162
#7  0x647271 in __gyre_r_root_MOD_solve_
        at /afs/mpa/data/earl/gyre/gyre-6.0.1/src/build/gyre_r_root.f90:132
#8  0x80362b in __gyre_bracket_search_MOD_bracket_search
        at /afs/mpa/data/earl/gyre/gyre-6.0.1/src/build/gyre_bracket_search.f90:153
#9  0x40a0d3 in gyre
        at /afs/mpa/data/earl/gyre/gyre-6.0.1/src/build/gyre.f90:288
#10  0x40bd4a in main
        at /afs/mpa/data/earl/gyre/gyre-6.0.1/src/build/gyre.f90:54
Floating point exception

Using the GYRE version that is bundled with the latest MESA.

Here is the .GYRE file and the gyre.in file: gyre-fpe-hecore-gmodes.zip

[warrickball]

The error itself seems to be cropping up because the boundary condition (BC) matrix is all zero (you can see this being set in src/ad/gyre_ad_bound.fpp). GYRE tries to rescale this matrix by the largest absolute value but that causes an FPE error (0/0, I guess). I'd need a bit more time to decipher why all the elements of the matrix are becoming small. One suggestion is to print out various quantities used to construct B in a place like this: https://github.com/rhdtownsend/gyre/blob/b0b5bc94bb12d21edaa23cfee2e10e1c909cfda2/src/ad/gyre_ad_bound.fpp#L703-L766 There (JCD BCs), only B(1,1), B(1,2) and B(1,3) vary.

I was a bit surprised to see V dropping dramatically near the surface of the star, to ~1, but I quickly evolved my own 17 Msun star to the same Teff (with an Eddington atmosphere) and saw the same behaviour. I always thought the limiting behaviour at the surface should be V → ∞ but I guess I'll have to revisit that for such hot stars. The sound speed inversion near the surface is presumably related. Either way, I'm not sure if this violates some assumption that goes into the BCs

[warrickball]

I also just wanted to flag that I've reproduced this with the 7.0 and the latest commit (b0b5bc94bb12d21edaa23cfee2e10e1c909cfda2).