ehsan wrote:

Thanks Rich for your thorough reply.

I admit that this was a pure exploratory experiment, to exploit the challenges that come up. Indeed, the limited floating point arithmetic is the limitation here. That is the reason we did not trust in what we saw in the output.

Thanks Rich for your thorough reply.

I've had a look through the inlist you sent me, and I've run it through GYRE with the corresponding model. The curious cases you have found unfortunately demonstrate what happens when you run GYRE outside its limits of validity! For the region of frequency space you are examining, the radial orders of the *adiabatic* modes are in the thousands. For such high-order modes, it is very difficult to obtain accurate numerical approximations to the corresponding *non-adiabatic* modes, due to the finite precision of computer floating-point arithmetic. As a result the non-adiabatic modes returned by GYRE, for the cases you consider, are in effect numerical noise. Increasing the grid resolution won't change this behavior; you've run into a limitation tied to the computer hardware.

I admit that this was a pure exploratory experiment, to exploit the challenges that come up. Indeed, the limited floating point arithmetic is the limitation here. That is the reason we did not trust in what we saw in the output.

It's possible to avoid these sorts of unnecessary calculations by choosing freq_min and freq_max to match the expected frequencies of modes with radial orders in some 'reasonable' range. This is what the ACOUSTC_DELTA and GRAVITY_DELTA frequency units are designed to do. So, for instance, if you're only interested in g modes with radial orders 1 <= n <= 100, you can do this (in GYRE 5.0):

Code:

`&scan`

grid_type = 'INVERSE'

freq_min = 0.01 ! 1/100

freq_max = 1 ! 1/1

freq_min_units = 'GRAVITY_DELTA'

freq_max_units = 'GRAVITY_DELTA'

n_freq = 500

/

This should serve as a cautionary tale: you can't treat GYRE like a black box! As a specific piece of advise for deciding whether you can believe GYRE's results, you should get into the habit of looking at the chi (convergence) parameters printed to the screen during a GYRE run. These should be in the typical range 1E-15 -> 1E-8. Much larger, or much smaller, suggests that GYRE either isn't converging, or is converging to a bogus solution.

Admittedly, the chi variable comes very handy in these cases, where we need to decide whether or not to rely on the GYRE output. The slim point here is that we need to run our GYRE computations on a computing cluster, and we typically do not sit in front of the screen looking at the GYRE STDOUT (which we should). We mainly look into the summary and eigenfunction files. Perhaps, this motivates putting chi as an additional output in the summary/eigenfunction files, so the user would not miss this valuable piece of info from the STDOUT.

Thanks again for enlightening us.

Ehsan.

I'll see if I can add chi to the variables which are written out to the files.

cheers,

Rich

Statistics: Posted by rhtownsend — Tue Apr 18, 2017 4:19 pm

]]>

I've had a look through the inlist you sent me, and I've run it through GYRE with the corresponding model. The curious cases you have found unfortunately demonstrate what happens when you run GYRE outside its limits of validity! For the region of frequency space you are examining, the radial orders of the *adiabatic* modes are in the thousands. For such high-order modes, it is very difficult to obtain accurate numerical approximations to the corresponding *non-adiabatic* modes, due to the finite precision of computer floating-point arithmetic. As a result the non-adiabatic modes returned by GYRE, for the cases you consider, are in effect numerical noise. Increasing the grid resolution won't change this behavior; you've run into a limitation tied to the computer hardware.

I admit that this was a pure exploratory experiment, to exploit the challenges that come up. Indeed, the limited floating point arithmetic is the limitation here. That is the reason we did not trust in what we saw in the output.

This should serve as a cautionary tale: you can't treat GYRE like a black box! As a specific piece of advise for deciding whether you can believe GYRE's results, you should get into the habit of looking at the chi (convergence) parameters printed to the screen during a GYRE run. These should be in the typical range 1E-15 -> 1E-8. Much larger, or much smaller, suggests that GYRE either isn't converging, or is converging to a bogus solution.

Admittedly, the chi variable comes very handy in these cases, where we need to decide whether or not to rely on the GYRE output. The slim point here is that we need to run our GYRE computations on a computing cluster, and we typically do not sit in front of the screen looking at the GYRE STDOUT (which we should). We mainly look into the summary and eigenfunction files. Perhaps, this motivates putting chi as an additional output in the summary/eigenfunction files, so the user would not miss this valuable piece of info from the STDOUT.

Thanks again for enlightening us.

Ehsan.

Statistics: Posted by ehsan — Tue Apr 18, 2017 3:09 am

]]>

ehsan wrote:

Dear Rich,

Mathias (and I) are making some tests with MESA+GYRE to look at the mode stability in core helium burning blue supergiant stars, and a possible trace of excitation by the epsilon mechanism. We have come across two curious cases:

I have hard time understanding two points here:

I also compared the work integral for epsilon-mechanism from Unno et al. (Eq. 25.9), and found it different from its implementation in line 1701 of the src/build/gyre_mode.f90 file. How do I track down that formula in Unno or other literature?

I will send you the model and inlist separately, as the file is large.

Kind regards,

Ehsan & Mathias

Dear Rich,

Mathias (and I) are making some tests with MESA+GYRE to look at the mode stability in core helium burning blue supergiant stars, and a possible trace of excitation by the epsilon mechanism. We have come across two curious cases:

- a g-mode penetrating the fully convective core, and having negligible amplitude elsewhere. this has a positive excitation
- a g-mode with positive dW_eps_dx in the convective core, followed by a negative dW_eps_dx in the radiative zone above the core

I have hard time understanding two points here:

- Why a g-mode would penetrate a convective core?
- Why dW_eps_dx be negative (while Unno et al. mentions that it only contributes towards destabilisation of a mode?

I also compared the work integral for epsilon-mechanism from Unno et al. (Eq. 25.9), and found it different from its implementation in line 1701 of the src/build/gyre_mode.f90 file. How do I track down that formula in Unno or other literature?

I will send you the model and inlist separately, as the file is large.

Kind regards,

Ehsan & Mathias

Hi Ehsan & Mathias --

I've had a look through the inlist you sent me, and I've run it through GYRE with the corresponding model. The curious cases you have found unfortunately demonstrate what happens when you run GYRE outside its limits of validity! For the region of frequency space you are examining, the radial orders of the *adiabatic* modes are in the thousands. For such high-order modes, it is very difficult to obtain accurate numerical approximations to the corresponding *non-adiabatic* modes, due to the finite precision of computer floating-point arithmetic. As a result the non-adiabatic modes returned by GYRE, for the cases you consider, are in effect numerical noise. Increasing the grid resolution won't change this behavior; you've run into a limitation tied to the computer hardware.

This should serve as a cautionary tale: you can't treat GYRE like a black box! As a specific piece of advise for deciding whether you can believe GYRE's results, you should get into the habit of looking at the chi (convergence) parameters printed to the screen during a GYRE run. These should be in the typical range 1E-15 -> 1E-8. Much larger, or much smaller, suggests that GYRE either isn't converging, or is converging to a bogus solution.

Regarding the work integral for the epsilon mechanism: Note that it's only the first (e_N) term inside the parentheses of Unno's eqn. 25.9 which is used in evaluating the nuclear contribution to the work integral.

cheers,

Rich

Statistics: Posted by rhtownsend — Sun Apr 16, 2017 1:50 pm

]]>

jingluan wrote:

Hi Rich,

Hi Rich,

Hi, Jing!

1, Which dependent variables and which set of ODE's do you use for non-adiabatic modes please? For example, (24.1) to (24.6) in Unno's are one choice of dependent variables, and equations (24.7) to (24.12) in Unno's are one choice of ODE's. The reason I like to ask this question is explained in the following part.

For the list of dependent variables and ODEs, see the file gyre/doc/equations.pdf included in the GYRE distribution. The ODEs are similar to those given in Unno et al, but instead of adopting equation (21.6) for the perturbed radiative heating term, I use (21.7).

2, Does the radiation diffusion make your ODE's a stiff system, especially when it is only a tiny little non-adiabatic please? I did a simple plane parallel model trying to check something, but I found the solution blows up (grows to very large amplitude and also oscillate crazily) due to the following equation,

i*omega*rho*kB*T*z0/F *\delta s = d(\delta F/F)/dz,

where z is the depth, and z0 is the length which I use to scale the depth, \delta s is the Lagrangian entropy perturbation, \delta F/F is the fractional radiative flux perturbation, I have converted partial derivative over time to i*omega. Because the coefficient, omega*rho*kB*T*z0/F, is very big, it brings in a branch of solutions which grows very fast with depth. The other branch declines very fast with depth, which is the physical branch I want. However, the fast growing branch blows up my solution. I tried some integration method for stiff system, but still failed.

So I am guessing maybe my choice of dependent variables and set of equations is ill-conditioned themselves??? That is why I ask about your choice.

You are correct -- when the solutions are only a little bit non-adiabatic, the ODEs are extremely stiff. This can make them very difficult to solve, especially when using shooting approaches which integrate from one boundary to the other. GYRE attempts to get around the stiffness issue by using multiple shooting (see Ascher, Mattheij & Russell, "Numerical Solution of Boundary Value Problems for Ordinary Differential Equations") and adopting implicit integration schemes.

I don't think this problem can be avoided by choosing a different set of dependent variables; the stiffness is intrinsic to the physics underlying the equations.

Thank you so much for reading this long message. Any suggestions, comments, guesses are highly appreciated please. Thank you

Happy to help, and also happy to talk a little about some of the challenges GYRE (and other oscillation codes) must overcome.

cheers,

Rich

Statistics: Posted by rhtownsend — Thu Apr 13, 2017 10:13 pm

]]>

Mathias (and I) are making some tests with MESA+GYRE to look at the mode stability in core helium burning blue supergiant stars, and a possible trace of excitation by the epsilon mechanism. We have come across two curious cases:

- a g-mode penetrating the fully convective core, and having negligible amplitude elsewhere. this has a positive excitation
- a g-mode with positive dW_eps_dx in the convective core, followed by a negative dW_eps_dx in the radiative zone above the core

I have hard time understanding two points here:

- Why a g-mode would penetrate a convective core?
- Why dW_eps_dx be negative (while Unno et al. mentions that it only contributes towards destabilisation of a mode?

I also compared the work integral for epsilon-mechanism from Unno et al. (Eq. 25.9), and found it different from its implementation in line 1701 of the src/build/gyre_mode.f90 file. How do I track down that formula in Unno or other literature?

I will send you the model and inlist separately, as the file is large.

Kind regards,

Ehsan & Mathias

Statistics: Posted by ehsan — Thu Apr 13, 2017 5:19 am

]]>

1, Which dependent variables and which set of ODE's do you use for non-adiabatic modes please? For example, (24.1) to (24.6) in Unno's are one choice of dependent variables, and equations (24.7) to (24.12) in Unno's are one choice of ODE's. The reason I like to ask this question is explained in the following part.

2, Does the radiation diffusion make your ODE's a stiff system, especially when it is only a tiny little non-adiabatic please? I did a simple plane parallel model trying to check something, but I found the solution blows up (grows to very large amplitude and also oscillate crazily) due to the following equation,

i*omega*rho*kB*T*z0/F *\delta s = d(\delta F/F)/dz,

where z is the depth, and z0 is the length which I use to scale the depth, \delta s is the Lagrangian entropy perturbation, \delta F/F is the fractional radiative flux perturbation, I have converted partial derivative over time to i*omega. Because the coefficient, omega*rho*kB*T*z0/F, is very big, it brings in a branch of solutions which grows very fast with depth. The other branch declines very fast with depth, which is the physical branch I want. However, the fast growing branch blows up my solution. I tried some integration method for stiff system, but still failed.

So I am guessing maybe my choice of dependent variables and set of equations is ill-conditioned themselves??? That is why I ask about your choice.

Thank you so much for reading this long message. Any suggestions, comments, guesses are highly appreciated please. Thank you

Sincerely,

Jing

Statistics: Posted by jingluan — Fri Mar 24, 2017 11:25 am

]]>

rhtownsend wrote:

Thanks for sharing your thoughts, Chris. From my math, the definition of dC_dx in gyre_mode.fpp means that

K_nl = dC/dx / C

...where C is the integral of dC_dx over 0 < x < 1. Likewise,

beta_nl = C

(but there's no factor of 4*pi, can you point out where that creeps in?).

I probably shouldn't have used C as my symbol, since the Ledoux splitting coefficient is

C_nl = 1 - beta_nl = 1 - C

(cf. Aerts 2010, 3.361). In fact, your suggestion of dbeta_dx is much better -- I'll make the change once we agree on the 4pi.

cheers,

Rich

mankovich wrote:I've been looking into this today; to recover a normalized rotation kernel K from dC_dx, you can divide dC_dx by (C * Rstar) where C is the integral of dC_dx from x=0 to x=1.

edit: I posted a confusing paragraph here earlier; suffice it to say that evidently C / (4*pi) gives beta_nl. Should that normalization be put in and these variables be renamed dbeta_dx and beta?

chris

Thanks for sharing your thoughts, Chris. From my math, the definition of dC_dx in gyre_mode.fpp means that

K_nl = dC/dx / C

...where C is the integral of dC_dx over 0 < x < 1. Likewise,

beta_nl = C

(but there's no factor of 4*pi, can you point out where that creeps in?).

I probably shouldn't have used C as my symbol, since the Ledoux splitting coefficient is

C_nl = 1 - beta_nl = 1 - C

(cf. Aerts 2010, 3.361). In fact, your suggestion of dbeta_dx is much better -- I'll make the change once we agree on the 4pi.

cheers,

Rich

OK, the 4pi comes in because of the way Aerts et al define the inertia in terms of \tilde{xi} = xi/sqrt(4pi), not xi. My expressions are correct... BUT I've decided to follow the Aerts convention, so I'll be making changes accordingly.

cheers,

Rich

Statistics: Posted by rhtownsend — Sat Mar 11, 2017 10:53 pm

]]>

Thanks for your help in figuring out the problem, I guess I'll upgrade soon.

Just to let you know, by setting "add_atmosphere_to_pulse_data" to false, the models that previously resulted in errors now seem to run correctly.

The GYRE model and MESA profile have the same amount of grind points, and the GYRE error message is gone.

I didn't check every calculated model, but the ones I checked run fine.

Best regards

Mathias

Statistics: Posted by Mathias — Fri Mar 10, 2017 6:39 am

]]>

Mathias wrote:

Hi Rich

Apparently, add_atmosphere is set to true internally, which explains the extra points in the surface layers.

Maybe I should set it to false for the post-MS models.

Best regards

Mathias

Hi Rich

Apparently, add_atmosphere is set to true internally, which explains the extra points in the surface layers.

Maybe I should set it to false for the post-MS models.

Best regards

Mathias

Hi Mathias --

Aha, you're running into a couple of bugs that we recently fixed in MESA. One bug adds on an atmosphere to the GYRE models, even when it isn't requested by the user; and the other causes mismatches between the atmosphere and the interior.

Both of these bugs are fixed in the latest release (9575) of MESA -- so, I'm afraid you'll have to upgrade to get a model that works OK.

cheers,

Rich

Statistics: Posted by rhtownsend — Wed Mar 08, 2017 3:59 pm

]]>

Apparently, add_atmosphere is set to true internally, which explains the extra points in the surface layers.

Maybe I should set it to false for the post-MS models.

Best regards

Mathias

Statistics: Posted by Mathias — Wed Mar 08, 2017 9:40 am

]]>

I also noticed it myself just now. The inlists I send you should be the ones used to produce the MESA profile and GYRE model.

I think the difference might be due to some extra options in the "run_star_extras" , but there are some methods used which are defined in other files, which I don't have on this computer right now. I'll check it out tomorrow when I have acces to the system used for my MESA runs.

Best regards

Mathias

Statistics: Posted by Mathias — Tue Mar 07, 2017 2:18 pm

]]>

Mathias wrote:

Hi Rich

This is the MESA profile corresponding to the problematic model file.

Best regards

Mathias

Hi Rich

This is the MESA profile corresponding to the problematic model file.

Best regards

Mathias

Hi Mathias --

Comparing the profile and the GYRE model, it seems the GYRE model has an extra 50 points in the surface layers, which presumably represent the stellar atmosphere. However, looking through your inlists I can't see any flags which tell MESA to add an atmosphere onto the GYRE model.

Is it possible that the inlists you sent me aren't exactly the same ones that you use to calculate the profile & GYRE model?

cheers,

Rich

Statistics: Posted by rhtownsend — Tue Mar 07, 2017 1:33 pm

]]>

This is the MESA profile corresponding to the problematic model file.

Best regards

Mathias

- profile.zip

Statistics: Posted by Mathias — Tue Mar 07, 2017 12:32 pm

]]>

Mathias wrote:

Dear Rich

The inlists are attached, I'm using MESA version 8845.

Best regards

Mathias

Dear Rich

The inlists are attached, I'm using MESA version 8845.

Best regards

Mathias

Hi Mathias --

One more thing -- can you post a MESA profile for the problem model?

cheers,

RIch

Statistics: Posted by rhtownsend — Mon Mar 06, 2017 9:32 am

]]>

The inlists are attached, I'm using MESA version 8845.

Best regards

Mathias

- inlists.zip

Statistics: Posted by Mathias — Mon Mar 06, 2017 5:05 am

]]>