NAD convergence problems of blue-looping stars

[afg]

I encounter and fail persistently to overcome convergence problems
when I try to compute nonadiabatic radial pulsations of
intermediate-mass stars in the classical instability strip while the
stars are on their blue loops.

I do manage to get nonadiabatic radial fundamental, first-, and
second-overtone results when the stars cross the strip for the first
time (i.e. as long as they are decently compact) but not anymore once
they are more puffed up.

Typically, what I get to see in a GYRE (V. 4.3 is running here) session is:

Code: Select all

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

Mode parameters
   l : 0
   m : 0

Building omega grid
   added scan interval :   0.1500000000000000E+01 ->   0.7500000000000000E+01 (60 points, LINEAR)

Building x grid
  x points : 1601 [after CREATE_MIDPOINT op]
  x range  :  0.0000000000000000E+00 ->   0.9999999985603535E+00

Root bracketing
  Time elapsed :     0.364 s

Root Solving
         l      n_pg       n_p       n_g                 Re(omega)                 Im(omega)                       chi  n_iter        n
         0         1         1         0    0.2655955222291173E+01    0.0000000000000000E+00    0.5018777496194273E-14       7     1601
         0         2         2         0    0.3889235307472945E+01    0.0000000000000000E+00    0.1058218887806180E-13       6     1601
         0         3         3         0    0.5344936627626918E+01    0.0000000000000000E+00    0.6300390309901204E-14       6     1601
  Time elapsed :      0.162 s
Root Solving
         l      n_pg       n_p       n_g                 Re(omega)                 Im(omega)                       chi  n_iter        n
         0         1         1         0    0.2661960010712412E+01    0.2627637059281112E-02    0.5084080344596849E-12       6     1601
Failed during deflate narrow : maximum iterations exceeded
   n_iter_def : 51
   n_iter     : 0
   omega_a    :   0.1583803223265424E+01 -0.1268621277039857E+01
   omega_b    :   0.1539169900795748E+01 -0.1271863390338167E+01
         0         3         3         0    0.5299325378482669E+01   -0.1708691639619444E-01    0.1246651790965083E-12       7     1601
  Time elapsed :      7.197 s

It is not always the 1st overtone, which goes haywire. It can be any of the 3 modes I am interested in.


Adiabatic computations are no problem at all, independent of the IVP
solver used (as suggested in the GYRE forum, I use usually one of the
COLLOC_GL<*> integrators. To me, the problem seems to be associated
with the (complex) root finding process. Actually, only the RIDDERS
method works reliably, if at all. To my surprise, even the SECANT
method fails most of the time.

To overcome the convergence problems I decided to try to set the
innermost grid point not at the star's center but resort to "envelope"
computations only (Just in case this all should be connected after all
to the full-star integration of the nonadiabatic problem). Therefore,
in the gyre.in file I request:

Code: Select all

&shoot_grid ! If left empty, adoption of defaults only
	op_type = 'CREATE_MIDPOINT'
	x_i     = 0.1
/

&recon_grid ! If left empty, adoption of defaults only
	op_type = 'CREATE_CLONE' 
	x_i     = 0.1
/
 ... and 

&osc
	inner_bound = 'ZERO'  

At startup, I get the stopping message:

Code: Select all

Root bracketing
 ASSERT 'this%x_i /= 0._WP' failed at line 210 <gyre_rad_bound:B_i_zero_>:
 Boundary condition invalid for x_i == 0

Checking the source of GYRE, I also got the impression that setting
x_i =/ 0 has no effect in nad/gyre_grid_par.f90. I admit though, that
I do not understand every single assignment in the routine.

Is it possible that x_i, x_o can be set via namelist but are not yet wired
for the intended use in the NAD case?

Any help is appreciated

Alfred

[rhtownsend]

Hi Alfred --

Please could you try upgrading to GYRE 4.4, and using the MAGNUS_GL2 scheme.

cheers,

Rich

[afg]

Rich

Thanks so much for Version 4.4. - and sorry for my sloppy realizing of the activities in the GYRE forum.

I just downloaded Version 4.4 and compiled it for standalone operation. This works very well. Now, I have convergence on the models which did not work before; i.e. I find continuous NAD mode orders with MAGNUS_GL2 with a convergence speed that is comparable with that on AD modes.

As of now, I just failed to embed the Version 4.4. in MESA so that I can call the correct GYRE version from within MESA. When I do a build_and_test from mesa/gyre/ (with the correct version plugged in in build_and_test) the process fails during the test stage. If Version 4.4. succeeds in all tests on your machines, I would be happy (also for the students next summer) if you could add V. 4.4. to the next MESA distribution!

I am extremely grateful for your efforts and your willingness to help. Now, I am very positive that I can successfully continue with my original plans for the summer school.

Alfred

[rhtownsend]

Hi Alfred --

Glad to hear things are working!

You'll be pleased to hear that GYRE 4.4 is already included in the development version of MESA, and we are planning a MESA release soon (in good time for the summer school).

cheers,

Rich

[afg]

Rich
I am just trying to scan the classical instability strip. Convergence with Version 4.4 is good, just as as tried out last night on two selected models.
However, I have the impression that the time dependence is back to exp(i*omega*t) as it was in the GYRE versions before 4.3 -- or somehow the instability strip has turned inside out . Is that choice a feature?

Alfred

[rhtownsend]

Rich
I am just trying to scan the classical instability strip. Convergence with Version 4.4 is good, just as as tried out last night on two selected models.
However, I have the impression that the time dependence is back to exp(i*omega*t) as it was in the GYRE versions before 4.3 -- or somehow the instability strip has turned inside out . Is that choice a feature?

Alfred

Hmmm, that's puzzling. The time dependence is still exp(-i*omega*t), as one can see from running the test cases in test/nad/mesa (the unstable modes show up as having Im(omega) > 0). So I'm not sure what's going on with your calculations.

cheers,

Rich

[afg]

If I do a GYRE standalone computation,I get the correct signs of omega_I (i.e. positive for unstable modes) so that indeed everything is fine, and in accordance to the post-V.4.3. definition of the time dependence. When I call GYRE from within MESA, I apparently mess things up. I tap the eigendata in gyre_support.f, there I have the following code segment (it is sufficient to look at bottom part starting with !>afg - experimental). The digits of omeg_r and omeg_i from gyre_support output (at the infamous write (*,*) checkpoint) agree exactly with those printed out in the standalone run of GYRE. The only exception is that omeg_i has the opposite sign to that from the standalone output (that's the origin of the question from the last posting).

Code: Select all

      
      subroutine gyre_call_back(md, ipar, rpar, ierr)
         use astero_data, only: store_new_oscillation_results
         
         type(mode_t), intent(in) :: md
         integer, intent(inout) :: ipar(:)
         real(dp), intent(inout) :: rpar(:)
         integer, intent(out) :: ierr

         integer :: new_el, new_order, new_em
         real(dp) :: new_inertia, new_linear_freq, omeg_r, omeg_i, &
                     dmp_freq, p_in_days, eta
         
         include 'formats'

         ierr = 0

         new_el          = md% mp% l
         new_order       = md% n_pg
         new_inertia     = md% E_norm()
         new_linear_freq = md% freq('UHZ')
         new_em          = 0
         
         call store_new_oscillation_results( &
            new_el, new_order, new_em, new_inertia, new_linear_freq, ierr)
         
        !>afg - experimental

         omeg_r   =  REALPART(md% omega)   ! dimensionless frequency
         omeg_i   =  IMAGPART(md% omega)

         write (*,*) omeg_i, omeg_r

         p_in_days = 1.15741d1/new_linear_freq
         eta       = omeg_i/omeg_r

Apparently I am doing something weird here, but what? If the digits of the omeg's would differ, I would presume that I tap the wrong variable; but the current situation just beats me.

Any advice of what I could do or check - or what not to touch?

Thanks, Alfred