Floating-point exception in bracket search for higher l

[jj_zanazzi]

profile150.data.GYRE

MESA input model

(3.57 MiB) Downloaded 2490 times

Hi Rich,

I hope you are doing well, I have been having a blast learning how to use GYRE. I have been able to run GYRE fine for l=2 oscillations to calculate g and f modes for both stars and planets, but when I try to go to higher l, I run into a floating-point exception when I think the code is closing in on a particular oscillation frequency. Depending on what I take to be the inner/outer boundary condition, the cutoff frequencies, as well as the number of frequencies in the scan, the code breaks at different times. I looked through the other reported floating-point error posts, and they don't seem to describe this one which occurs during the bracket search.

Below is the output from a run I did, setting l=6:

Code: Select all

gyre [7.1]
----------

OpenMP Threads   : 1
Input filename   : gyre.in

Model Init
----------

Reading from MESA file
   File name LOGS_to_end_core_h_burn/profile150.data.GYRE
   File version 1.01
   Read 7589 points
   No need to add central point

Mode Search
-----------

Mode parameters
   l : 6
   m : 0

Building frequency grid (REAL axis)
   added scan interval :  0.1000E-01 ->  0.1000E+01 (2000 points, INVERSE)

Building spatial grid
   Scaffold grid from model
   Refined 2432 subinterval(s) in iteration 1
   Refined 3590 subinterval(s) in iteration 2
   Refined 2772 subinterval(s) in iteration 3
   Refined 2010 subinterval(s) in iteration 4
   Refined 2052 subinterval(s) in iteration 5
   Refined 2204 subinterval(s) in iteration 6
   Refined 1699 subinterval(s) in iteration 7
   Refined 516 subinterval(s) in iteration 8
   Refined 24 subinterval(s) in iteration 9
   Refined 20 subinterval(s) in iteration 10
   Refined 20 subinterval(s) in iteration 11
   Refined 19 subinterval(s) in iteration 12
   Refined 20 subinterval(s) in iteration 13
   Refined 20 subinterval(s) in iteration 14
   Refined 20 subinterval(s) in iteration 15
   Refined 0 subinterval(s) in iteration 16
   Final grid has 1 segment(s) and 25007 point(s):
      Segment 1 : x range 0.0000 -> 1.0000 (1 -> 25007)

Starting search (adiabatic)

Evaluating discriminant
  Time elapsed :   329.207 s

Root Solving
   l    m    n_pg    n_p    n_g       Re(omega)       Im(omega)        chi n_iter
   6    0   -2718      0   2718  0.10000853E-01  0.00000000E+00 0.2172E+00     27
   6    0   -2714      0   2714  0.10009915E-01  0.00000000E+00 0.1354E+00     19
   6    0   -2713      0   2713  0.10011690E-01  0.00000000E+00 0.4615E+00     34
   6    0   -2710      0   2710  0.10029174E-01  0.00000000E+00 0.1270E-01     26
   6    0   -2711      0   2711  0.10030527E-01  0.00000000E+00 0.2416E+00     36
   6    0   -2705      1   2706  0.10043512E-01  0.00000000E+00 0.3085E+00     41
   6    0   -2693      0   2693  0.10093597E-01  0.00000000E+00 0.1924E+00     30
   6    0   -2690      1   2691  0.10103939E-01  0.00000000E+00 0.1167E-01     37
   6    0   -2689      0   2689  0.10114367E-01  0.00000000E+00 0.2648E+00     43
   6    0   -2684      0   2684  0.10132367E-01  0.00000000E+00 0.3789E-01     40
   6    0   -2679      0   2679  0.10140582E-01  0.00000000E+00 0.2578E-02     19
   6    0   -2678      0   2678  0.10145564E-01  0.00000000E+00 0.4677E+00     38
   6    0   -2677      0   2677  0.10147450E-01  0.00000000E+00 0.1258E+01     40
   6    0   -2674      0   2674  0.10160653E-01  0.00000000E+00 0.1564E+00     31
   6    0   -2676      0   2676  0.10166059E-01  0.00000000E+00 0.1139E+02     22
   6    0   -2672      0   2672  0.10171165E-01  0.00000000E+00 0.5279E-01     47
   6    0   -2673      0   2673  0.10171540E-01  0.00000000E+00 0.8719E-01     36
   6    0   -2671      0   2671  0.10176910E-01  0.00000000E+00 0.1135E+00     30
   6    0   -2668      0   2668  0.10191284E-01  0.00000000E+00 0.1254E-01     26
   6    0   -2663      2   2665  0.10196019E-01  0.00000000E+00 0.1589E+00     29
   6    0   -2663      0   2663  0.10201544E-01  0.00000000E+00 0.8886E-01     32
   6    0   -2665      0   2665  0.10205699E-01  0.00000000E+00 0.1341E+00     40
   6    0   -2664      0   2664  0.10210190E-01  0.00000000E+00 0.4943E+00     42
   6    0   -2661      0   2661  0.10219730E-01  0.00000000E+00 0.1267E+00     30
   6    0   -2659      0   2659  0.10227907E-01  0.00000000E+00 0.2176E-01     39
   6    0   -2656      1   2657  0.10233125E-01  0.00000000E+00 0.6806E+00     31
   6    0   -2655      0   2655  0.10233174E-01  0.00000000E+00 0.0000E+00     10
   6    0   -2651      0   2651  0.10250699E-01  0.00000000E+00 0.8865E-01     29
   6    0   -2649      0   2649  0.10259081E-01  0.00000000E+00 0.9792E-01     26
   6    0   -2646      0   2646  0.10275843E-01  0.00000000E+00 0.1078E-01     40
   6    0   -2644      0   2644  0.10284984E-01  0.00000000E+00 0.1314E-01     33
   6    0   -2644      0   2644  0.10289130E-01  0.00000000E+00 0.8599E-02     27
   6    0   -2636      0   2636  0.10304942E-01  0.00000000E+00 0.1805E+00     41
   6    0   -2637      0   2637  0.10306252E-01  0.00000000E+00 0.2143E+00     41
   6    0   -2635      1   2636  0.10311927E-01  0.00000000E+00 0.8024E+00     32
   6    0   -2634      1   2635  0.10318300E-01  0.00000000E+00 0.1816E+00     30
   6    0   -2630      0   2630  0.10336438E-01  0.00000000E+00 0.1385E+00     38
   6    0   -2629      0   2629  0.10343723E-01  0.00000000E+00 0.2442E+00     37
   6    0   -2625      0   2625  0.10350636E-01  0.00000000E+00 0.1504E-07     41
   6    0   -2621      1   2622  0.10364265E-01  0.00000000E+00 0.2029E+00     27
   6    0   -2622      0   2622  0.10369710E-01  0.00000000E+00 0.3185E+00     23
   6    0   -2619      0   2619  0.10375308E-01  0.00000000E+00 0.6544E-01     32
   6    0   -2617      0   2617  0.10391103E-01  0.00000000E+00 0.1183E+00     31
   6    0   -2617      0   2617  0.10391517E-01  0.00000000E+00 0.1999E-01     32
   6    0   -2612      0   2612  0.10399543E-01  0.00000000E+00 0.1280E+00     31
   6    0   -2610      0   2610  0.10411598E-01  0.00000000E+00 0.2088E-01     35
   6    0   -2612      0   2612  0.10413265E-01  0.00000000E+00 0.1223E+00     29
   6    0   -2606      0   2606  0.10433869E-01  0.00000000E+00 0.4872E+00     38
   6    0   -2599      0   2599  0.10456350E-01  0.00000000E+00 0.2204E+00     38
   6    0   -2599      0   2599  0.10466110E-01  0.00000000E+00 0.1924E-01     40
   6    0   -2596      0   2596  0.10469126E-01  0.00000000E+00 0.3002E+03     19
   6    0   -2597      0   2597  0.10472311E-01  0.00000000E+00 0.8328E-01     37
   6    0   -2588      0   2588  0.10502142E-01  0.00000000E+00 0.2168E-07     28
   6    0   -2585      0   2585  0.10511876E-01  0.00000000E+00 0.1199E+01     35

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

Backtrace for this error:
#0  0x7febda91a08f in ???
#1  0x5226f0 in __gyre_wave_MOD_f_t
#2  0x53d7a0 in __gyre_summary_MOD_cache
#3  0x409556 in process_mode_ad.1
#4  0x4d1632 in __gyre_bracket_search_MOD_bracket_search
#5  0x4071a0 in MAIN__
#6  0x405bb6 in main
[4]+  Done                    gedit gyre.in
Floating point exception (core dumped)

Here is my gyre inlist:

Code: Select all

&constants
/

&model
  model_type = 'EVOL' 
  file = 'profile150.data.GYRE'
  file_format = 'MESA' 
  add_center = .True.
/

&mode
  l = 6 ! Harmonic degree
/


&osc
  inner_bound = 'REGULAR' 
  outer_bound = 'VACUUM'
  nonadiabatic = .TRUE.
/

&rot
/

&num
  diff_scheme = 'MAGNUS_GL4' !
  nad_search = 'AD'
/

&scan
  grid_type = 'INVERSE' 
  freq_min = 0.01      
  freq_max = 1.0      
  n_freq = 2000         
  axis = 'REAL'
/

&grid
  w_osc = 10 ! Oscillatory region weight parameter
  w_exp = 2  ! Exponential region weight parameter
  w_ctr = 10 ! Central region weight parameter
  w_str = 10 ! Structure parameter
/


&ad_output
  summary_file = 'summary/summary_ad.h5'                       
  summary_item_list = 'M_star,R_star,l,n_pg,omega,E_norm,f_T,psi_T,Delta_g' 
  detail_item_list = 'l,n_pg,omega,x,xi_r,xi_h,eul_Phi,eul_rho,lag_L,lag_T,rho' 	     
/

&nad_output
  summary_file = 'summary/summary_nad.h5'                         
  summary_item_list = 'M_star,R_star,l,n_pg,omega,E_norm,f_T,psi_T,Delta_g' 
  detail_template = 'detail/detail.l%l.n%n.h5'        	   
  detail_item_list = 'l,n_pg,omega,x,xi_r,xi_h,eul_Phi,eul_rho,lag_L,lag_T,rho' 	   
/

I have also attached the mesa model it is being run on.

Is there something silly with my GYRE inlist, or a problem with my MESA model? I am using the most recent GYRE (7.1), MESA (r24.03.1), and MESA SDK (23.7.3) versions.

Thanks again for all your work on GYRE,
J. J. Zanazzi.

[warrick]

Hi JJ,

You can get the backtrace to show you line numbers of most of the source code in it by building GYRE with `DEBUG=yes make` instead of just `make`. I did so and reproduced the error, now with the backtrace

Code: Select all

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

Backtrace for this error:
#0  0x7f35ede237f2 in ???
#1  0x7f35ede22985 in ???
#2  0x7f35edc5c99f in ???
#3  0x5cf95d in __gyre_wave_MOD_f_t
	at /home/wball/code/gyre/7.1/src/build/gyre_wave.f90:2262
#4  0x69ebab in cache_wave_
	at /home/wball/code/gyre/7.1/src/build/gyre_summary.f90:518
#5  0x69fc39 in __gyre_summary_MOD_cache
	at /home/wball/code/gyre/7.1/src/build/gyre_summary.f90:418
#6  0x40b1b4 in process_mode_ad
	at /home/wball/code/gyre/7.1/src/build/gyre.f90:426
#7  0x635aa9 in __gyre_bracket_search_MOD_bracket_search
	at /home/wball/code/gyre/7.1/src/build/gyre_bracket_search.f90:186
#8  0x40866a in gyre
	at /home/wball/code/gyre/7.1/src/build/gyre.f90:289
#9  0x40b27a in main
	at /home/wball/code/gyre/7.1/src/build/gyre.f90:54

The offending lines in gyre_wave.f90 are

Code: Select all

  function f_T (this)

    class(wave_t), intent(in) :: this
    real(WP)                  :: f_T

    complex(WP) :: C_T

    ! Evaluate the non-adiabatic f_T parameter. This is expression is                                                                                                                                                                                         
    ! based on eqn. (5) of [Dupret:2003]                                                                                                                                                                                                                      

    associate (j => this%j_ref)

      C_T = this%lag_T_eff()/this%xi_r(j) ! this is line 2262

      f_T = abs(C_T)

    end associate

    ! Finish                                                                                                                                                                                                                                                  

    return

  end function f_T

I dumped the values of the top and bottom of the quotient and that particular mode has zero for both, hence the floating-point error.

I'm not very familiar with non-adiabatic calculations. My hunch would be that C_T should be zero is lag_T_eff is but you might have a better idea. The stellar model looks fine (it's always the first thing I check) and the facts that other modes work and the failing mode depends on the input parameters makes me thing this isn't because you're trying to compute something unreasonable.

[rhtownsend]

Hi JJ --

Further to Warrick's comments, it looks like both lag_T_eff and xi_r are zero at the stellar surface, causing the error. This shouldn't usually happen, but may be a consequence of using the MAGNUS_GL4 diff scheme. Could you try using MAGNUS_GL2 instead and see if that works better?

(The higher-order Magnus schemes can sometimes result in rather crappy eigenfunctions).

cheers,

Rich

[jj_zanazzi]

Hey folks,

@warrick, thanks so much for finding the offending lines! I ran my code asking it to not compute f_T, and it works fine. I will make sure to always run it with debug=yes to find the errors myself from now on.

@rhtownsend, I did try switching to MAGNUS_GL2, and it gives better eigenfunctions, but still crashes when xi_r is zero at the surface. I will make sure to run my code with MAGNUS_GL2 from now on though.

Thanks so much for your help on this issue!

Best wishes,
J. J.

[rhtownsend]

Hi JJ --

I'll take a look at the crashing case in the coming days, to see if I can pinpoint why things are going awry.

cheers,

Rich