GYRE Floating point exception and grid refinement

[mathren90]

Hi everyone,

I am trying to verify the stability of some of my MESA models by looking at the imaginary parts of their eigen-frequencies, calculated with GYRE. This is the first time I try to use this tool and I am encountering some problems, described in the following.

First of all I am using:

Code: Select all

gfortran --version
GNU Fortran (GCC) 4.10.0 20140710 (experimental)
Copyright (C) 2014 Free Software Foundation, Inc.

and I tried both the GYRE version shipped with MESA 6794 and GYRE 4.0 (installed and tested following the tutorial here: https://bitbucket.org/rhdtownsend/gyre/ ... %20Started ). I encountered similar issues with both. The input file is a MESA 6794 model at oxygen depletion. I produced the GYRE profile by setting in the MESA controls namelist:

Code: Select all

write_pulse_info_with_profile = .true.
pulse_info_format = 'GYRE'

You can find the input profile I am using for this tests in a tarball here: http://www.tapir.caltech.edu/~mathren90 ... YRE.tar.gz

With the attached parameter file, I am able to evaluate only the modes with l=0, while modes with l=1,2 (no matter what m is set) give a different problem depending on the ivp_solver used:

* if

Code: Select all

ivp_solver = 'MAGNUS_GL2'

I am able to obtain results, but the n_pg values are not ordered in a monotonic-increasing fashion. To fix this, it is my understanding that I need a finer grid. However, if I try to set

Code: Select all

op_type = 'CREATE_UNIFORM' 

and a very large n in the &shoot_grid and clone this grid in the &recon_grid, I get (this also for l=0) the following output:

Code: Select all

     ASSERT 'i > 0 .AND. i < this%n' failed at line 433 <gyre_spline:f_1_>:
     Out-of-bounds interpolation

* if

Code: Select all

ivp_solver = 'MAGNUS_GL4'

or

Code: Select all

'MAGNUS_GL6' 

the l=0 modes give regular output cloning the grid of the input model, but modes with higher l(=1,2) give

Code: Select all

Mode Search
===========

Mode parameters
   l : 2
   m : 0

Building omega grid
  omega points : 1000
  omega range  :   0.1500000000000000E+00 ->   0.1500000000000000E+03

Building x grid
  x points : 4077 [after CREATE_CLONE op]
  x range  :  0.0000000000000000E+00 ->   0.9999999999986838E+00

Root bracketing

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

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

Backtrace for this error:

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

Backtrace for this error:

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

Backtrace for this error:

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

Backtrace for this error:

Backtrace for this error:

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

Backtrace for this error:
#0  0x2B2100520517
#1  0x2B2100520B10
#2  0x2B210136E02F
#3  0x4BABBA in __gyre_r_ext_MOD_exp_
#4  0x4E08B7 in __gyre_r_magnus_ivp_MOD_shoot_
#5  0x4D1ABF in __gyre_r_bvp_MOD_build_._omp_fn.2 at gyre_r_bvp.f90:0
#6  0x2B2100AC2CDD
#7  0x2B210112883C
#8  0x2B2101412FCC
#0  0x2B2100520517
#1  0x2B2100520B10
#2  0x2B210136E02F
#3  0x4BABBA in __gyre_r_ext_MOD_exp_
#4  0x4E08B7 in __gyre_r_magnus_ivp_MOD_shoot_
#5  0x4D1ABF in __gyre_r_bvp_MOD_build_._omp_fn.2 at gyre_r_bvp.f90:0
#6  0x2B2100AC2CDD
#7  0x2B210112883C
#8  0x2B2101412FCC
#0  0x2B2100520517
#1  0x2B2100520B10
#2  0x2B210136E02F
#3  0x4BABBA in __gyre_r_ext_MOD_exp_
#4  0x4E08B7 in __gyre_r_magnus_ivp_MOD_shoot_
#5  0x4D1ABF in __gyre_r_bvp_MOD_build_._omp_fn.2 at gyre_r_bvp.f90:0
#0  0x2B2100520517
#1  0x2B2100520B10
#2  0x2B210136E02F
#3  0x4BABBA in __gyre_r_ext_MOD_exp_#0  0x2B2100520517

#1  0x2B2100520B10
#2  0x2B210136E02F
#4  0x4E08B7 in __gyre_r_magnus_ivp_MOD_shoot_
#3  0x4BABBA in __gyre_r_ext_MOD_exp_
#6  0x#5  0x4D1ABF in __gyre_r_bvp_MOD_build_._omp_fn.2 at gyre_r_bvp.f90:0
2B2100AC2CDD
#6  0x2B2100AC2CDD#4  0x4E08B7 in __gyre_r_magnus_ivp_MOD_shoot_

#7  0x2B210112883C
#5  0x4D1ABF#7  0x2B210112883C
#8  0x2B2101412FCC
#8  0x2B2101412FCC

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

Backtrace for this error:
 in __gyre_r_bvp_MOD_build_._omp_fn.2 at gyre_r_bvp.f90:0

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

Backtrace for this error:
#6  0x2B2100AC2CDD
#7  0x2B210112883C
#8  0x2B2101412FCC
#0  0x2B2100520517
#1  0x2B2100520B10
#2  0x2B210136E02F

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

Backtrace for this error:
#0  0x2B2100520517
#1  0x2B2100520B10
#2  0x2B210136E02F
#3  #3  0x4BABBA in __gyre_r_ext_MOD_exp_
0x4BABBA in #4  0x4E08B7 in __gyre_r_magnus_ivp_MOD_shoot_
__gyre_r_ext_MOD_exp_
Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:

#5  0x4D1ABF in __gyre_r_bvp_MOD_build_._omp_fn.2 at gyre_r_bvp.f90:0
#6  0x2B2100AC2CDD
#7  0x2B210112883C
#8  0x2B2101412FCC

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

Backtrace for this error:
#0  0x2B2100520517
#0  0x2B2100520517
#1  0x2B2100520B10
#2  #0  0x2B21005205170x2B210136E02F
#4  0x4E08B7 in __gyre_r_magnus_ivp_MOD_shoot_


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

Backtrace for this error:
#1  0x2B2100520B10
#2  0x#3  0x4BABBA in 2B210136E02F
__gyre_r_ext_MOD_exp_
#1  0x2B2100520B10#4  0x4E08B7 in __gyre_r_magnus_ivp_MOD_shoot_
#3  0x4BABBA in __gyre_r_ext_MOD_exp_

#5  0x4D1ABF in __gyre_r_bvp_MOD_build_._omp_fn.2 at gyre_r_bvp.f90:0
#6  0x2B2100AC2CDD
#7  0x2B210112883C
#8  0x2B2101412FCC
#5  0x4D1ABF in __gyre_r_bvp_MOD_build_._omp_fn.2 at gyre_r_bvp.f90:0
#2  #4  0x4E08B7 in __gyre_r_magnus_ivp_MOD_shoot_
0x#5  0x4D1ABF in __gyre_r_bvp_MOD_build_._omp_fn.2 at gyre_r_bvp.f90:0
2B210136E02F#6  0x2B2100AC2CDD

#7  0x2B210112883C
#8  0x2B2101412FCC
#3  0x4BABBA in __gyre_r_ext_MOD_exp_
#6  0x2B2100AC2CDD
#7  0x2B210112883C
#8  0x2B2101412FCC
#4  0x4E08B7 in __gyre_r_magnus_ivp_MOD_shoot_
#5  0x4D1ABF in __gyre_r_bvp_MOD_build_._omp_fn.2 at gyre_r_bvp.f90:0
#6  0x2B2100AC2CDD
#7  0x2B210112883C
#8  0x2B2101412FCC
#0  0x2B2100520517
#1  0x2B2100520B10
#2  0x2B210136E02F
#3  0x4BABBA in __gyre_r_ext_MOD_exp_
#4  0x4E08B7 in __gyre_r_magnus_ivp_MOD_shoot_
#0  0x2B2100520517
#1  0x2B2100520B10
#2  0x2B210136E02F
#5  0x4D1ABF in __gyre_r_bvp_MOD_build_._omp_fn.2 at gyre_r_bvp.f90:0
#6  0x2B2100ABFBAE
#3  #7  0x4BABBA in __gyre_r_ext_MOD_exp_0x
4D1F1D in __gyre_r_bvp_MOD_build_
#4  0x4E08B7 in __gyre_r_magnus_ivp_MOD_shoot_
#8  0x4CFF2A in __gyre_r_bvp_MOD_discrim_
#5  0x4D1ABF in __gyre_r_bvp_MOD_build_._omp_fn.2 at gyre_r_bvp.f90:0
#6  0x2B2100AC2CDD
#7  0x2B210112883C
#8  0x2B2101412FCC
#9  0x44EEFE in __gyre_r_discfunc_MOD_eval_
#10  0x41A21B in __gyre_r_search_MOD_scan_search
#11  0x40404E in MAIN__ at gyre_ad.f90:0
Floating point exception

This happens with both GYRE version I tried, and with different outer boundary conditions, frequency interval and number of frequencies (n_freq) . Since this is my first attempt with GYRE, it is very likely that I am doing naive errors, and I would be glad if you could help me tracking them down. My aim is to find the imaginary part of the eigenfrequencies of my models (once I can run GYRE I would do this for a large grid of MESA models at different stages of their evolution), I would be glad to hear any suggestion on this problem too if you have any.

Mathieu Renzo

[rhtownsend]

Hi Mathieu --

Thanks for your post. I'm pretty sure I can help you fix this -- but first, could you re-upload the model you're using? The tarball you linked to contains garbage.

cheers,

Rich

[mathren90]

Oh, I am so sorry. I don't know how it happened. Now the tarball content should be right.

[rhtownsend]

OK, thanks, the tarball now works.

Looking through your namelist input file, there are a couple of potential issues I can spot:

1) You're using a CREATE_CLONE shooting grid, and nothing else. You'll certainly want to use a higher-resolution grid, but CREATE_UNIFORM isn't the way to do that. Have you read the tutorial on grids? If so, you'll see that the correct approach is to put in an additional namelist looking like this:

Code: Select all

&shoot_grid
	op_type = 'RESAMP_DISPERSION'	! Resample the grid based on the local dispersion relation
	alpha_osc = 10			! At least 5 points per oscillatory wavelength
	alpha_exp = 2			! At least 1 point per exponential 'wavelength'
/

2) Your frequency scan range is way too wide -- you'll end up getting very high-order modes, which are almost certainly not what your looking for (and will cause the code to crash due to too little RAM). Do you understand what the difference between 'INVERSE' and 'LINEAR' is for grid_type? And do you understand what freq_units = 'NONE' corresponds to? More generally, what sort of mode frequencies are you hoping to find? Answering this latter question will help guide a better choice of the frequency scan.

[mathren90]

Thank you very much! This is very helpful.

1) You're using a CREATE_CLONE shooting grid, and nothing else. You'll certainly want to use a higher-resolution grid, but CREATE_UNIFORM isn't the way to do that. Have you read the tutorial on grids? If so, you'll see that the correct approach is to put in an additional namelist looking like this:

Code: Select all

&shoot_grid
	op_type = 'RESAMP_DISPERSION'	! Resample the grid based on the local dispersion relation
	alpha_osc = 10			! At least 5 points per oscillatory wavelength
	alpha_exp = 2			! At least 1 point per exponential 'wavelength'
/

I must admit I have missed that tutorial. According to it, I would not need to re-grid for the reconstruction grid, since a rough characterization of the spectrum is exactly what I am looking for (at least for the moment).

2) Your frequency scan range is way too wide -- you'll end up getting very high-order modes, which are almost certainly not what your looking for (and will cause the code to crash due to too little RAM). Do you understand what the difference between 'INVERSE' and 'LINEAR' is for grid_type? And do you understand what freq_units = 'NONE' corresponds to? More generally, what sort of mode frequencies are you hoping to find? Answering this latter question will help guide a better choice of the frequency scan.

I tried to vary the frequency range while playing with the code to understand how it works. If I am correct, grid_type determines if we work in the period domain ('INVERSE') or frequency domain ('LINEAR'); and freq_units = 'NONE' means the freq_max and freq_mix are interpreted as being dimensionless as if they where normalized with the inverse of the free-fall time (eq. A4 of the gyre paper ). I was trying to sample from 1/100 of the dynamical frequency to 100 times it (but I agree this is a large interval, and I probably only need a much smaller one). How much is this wrong?

My aim is to investigate the dynamical stability of my models, so in principle I want to find the modes whose frequencies have the largest imaginary parts (i.e. the largest growth rate for the associated instability). Hopefully I will not find any, especially not in the test case in the tarball, which is for a "nice" 15 solar mass model that I (intuitively) expect to be stable. Probably I can limit myself to radial and l=1 modes for the time being, and I would expect only the firsts few order modes to be significant for my present purposes.

I added the grid refinement based on the local dispersion relation (using alpha_osc between 1 and 10, and the same for alpha_exp), I tried narrowing the frequency interval (freq_min>0.1, freq_max<2, freq_units = 'NONE'), changed the n_freq between a few tenths and a few hundreds, and I tried different solvers too. So far I haven't been able to get results, however I never run into "floating point exceptions" anymore. What happens is that both the grid construction and the root bracketing (for l=1) get extremely long (for the moment I am using a single node with 12 cores on a cluster). Also, depending on the values I am trying, I sometimes don't find any radial mode (but I suspect this is because in these cases I am looking to a too narrow frequency interval probably). I was expecting that using a small n_freq(~10) the root bracketing should not be so slow... I am certainly misunderstanding something in the program behavior here, and I would be glad if you could give me some more pointers in the right direction, especially for the choice of the frequency interval and its sampling.

Thank you again!

Mathieu