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:

gyre/src/ad/gyre_ad_bound.fpp

Lines 703 to 766 in b0b5bc9

	subroutine build_jcd_o_ (this, st, B, scl)
	class(ad_bound_t), intent(in) :: this
	class(r_state_t), intent(in) :: st
	real(WP), intent(out) :: B(:,:)
	real(WP), intent(out) :: scl(:)
	real(WP) :: Omega_rot
	real(WP) :: omega_c
	real(WP) :: lambda
	real(WP) :: l_e
	real(WP) :: chi
	real(WP) :: b_11
	real(WP) :: b_12
	$CHECK_BOUNDS(SIZE(B, 1),this%n_o)
	$CHECK_BOUNDS(SIZE(B, 2),this%n_e)
	$CHECK_BOUNDS(SIZE(scl),this%n_o)
	! Evaluate the outer boundary conditions ([Chr2008] formulation)
	associate( &
	V => this%coeff(2,J_V), &
	As => this%coeff(2,J_AS), &
	c_1 => this%coeff(2,J_C_1), &
	Gamma_1 => this%coeff(2,J_GAMMA_1), &
	pt => this%pt(2), &
	alpha_grv => this%alpha_grv, &
	alpha_omg => this%alpha_omg)
	Omega_rot = this%cx%Omega_rot(pt)
	omega_c = this%cx%omega_c(Omega_rot, st)
	lambda = this%cx%lambda(Omega_rot, st)
	l_e = this%cx%l_e(Omega_rot, st)
	chi = atmos_chi(V, As, c_1, Gamma_1, omega_c, lambda)
	b_11 = V/Gamma_1 - 3._WP
	b_12 = lambda/(c_1*alpha_omg*omega_c**2) - V/Gamma_1
	! Set up the boundary conditions
	B(1,1) = chi - b_11
	B(1,2) = -b_12
	B(1,3) = alpha_grv*((lambda/(c_1*alpha_omg*omega_c**2) - l_e - 1._WP)*b_12/(V/Gamma_1 + As))
	B(1,4) = alpha_grv*(0._WP)
	B(2,1) = alpha_grv*(0._WP)
	B(2,2) = alpha_grv*(0._WP)
	B(2,3) = alpha_grv*(l_e + 1._WP) + (1._WP - alpha_grv)
	B(2,4) = alpha_grv*(1._WP)
	scl = 1._WP
	end associate
	! Finish
	return
	end subroutine build_jcd_o_

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 (b0b5bc9).

[rhdtownsend]

Addressing this old issue: for the given model and frequency, the atmospheric boundary condition is being evaluated way below the lower atmospheric cutoff frequency. Therefore, we're in a regime where the boundary condition is leaky, and no strictly-standing solutions are possible.

GYRE sort-of flags this earlier on in the calculation output:

INFO at line 93 of /Users/townsend/devel/gyre/src/common/atmos_m.fypp
Discarding imaginary part of atmospheric radial wavenumber

Basically, the boundary condition is complex, but the imaginary part gets thrown away (because adiabatic calculations are performed using real arithmetic).

So, this isn't so much a bug, as using GYRE to model a situation that's inherently unphysical. Still, it would be good to have GYRE do something more elegant than crash!

[warrickball]

... the atmospheric boundary condition is being evaluated way below the lower atmospheric cutoff frequency.

I don't understand the physics/mathematics of what's happening here but thought this argument was odd because the lower-frequency reportedly do work.

Still, it would be good to have GYRE do something more elegant than crash!

Is there some condition we can evaluate to print a warning like "this doesn't make sense, consider using some_control = 'OTHER_VALUE' instead"?