Hi Rich,
I'm having some trouble with non adiabatic calculations of white dwarf models. Basically, I expect dW_dx to be small everywhere except in the outer layers, but the opposite seems to be the case. I'm attaching plots of example calculations I've done of a 0.15 Msun WD.
The bumps that I see at log (1-Mr/Mstar) = -13 are about what I'd expect, but everything below (1-Mr/Mstar) = -10 looks suspect. The three plots have theta_ad=0, 1.e-03, and 1.e-06. Clearly, they are all quite different. I'm also attaching a plot of delS, which looks very well behaved, so I'm wondering what's going on with the dW_dx calculation. Do you have any ideas?
Cheers,
Mike
Non-adiabatic driving
-
rhtownsend
- Site Admin
- Posts: 397
- Joined: Sun Mar 31, 2013 4:22 pm
Re: Non-adiabatic driving
Hi Mike --
As I mentioned in a separate email, this undesirable behavior stems from numerical noise in the entropy perturbation deep in the interior. This noise in turn stems from the fact that the non-adiabatic pulsation equations are quite ill-behaved in regions where the thermal timescale is long; although the equations are formally 6th order, they 'want' to behave like a 4th order (adiabatic) system.
There have been various attempts in the literature to get around the bad behavior. For instance, Unno et al. (1989) discuss a modification to the Henyey scheme proposed by Sugimoto (1970), which involves discretizing the thermal equations using asymmetric finite differences; my experience with older (pre-GYRE) codes has been that this works quite well in suppressing the noise. In his PhD thesis, Marc-Antoine Dupret discusses an alternative approach used in his MAD code, which also seems to work well.
For the multiple shooting scheme used by GYRE, the literature is unfortunately pretty sparse and I'm still looking for an efficient fix for the noise. In the meantime the noise can generally be eliminated by increasing the grid resolution where the thermal timescale is long. To do this, use the RESAMP_THERMAL operation when setting up the shooting grid:
&shoot_grid
op_type = 'CREATE_CLONE'
/
&shoot_grid
op_type = 'RESAMP_THERMAL'
alpha_thm = 1E-2
/
The key parameter here is alpha_thm (introduced in release 1.1 of the code); choosing larger values will increase the grid resolution, with the benefit of reduced entropy noise and the drawback of longer execution times.
Unfortunately, in the specific case of the input files you've sent me, it seems the number of grid points required for an accurate entropy eigenfunction is extremely large; here, the cure may be worse than the disease!
The upshot is: you can try playing around with alpha_thm (after upgrading to release 1.1), but I imagine a better fix will be necessary.
cheers,
Rich
As I mentioned in a separate email, this undesirable behavior stems from numerical noise in the entropy perturbation deep in the interior. This noise in turn stems from the fact that the non-adiabatic pulsation equations are quite ill-behaved in regions where the thermal timescale is long; although the equations are formally 6th order, they 'want' to behave like a 4th order (adiabatic) system.
There have been various attempts in the literature to get around the bad behavior. For instance, Unno et al. (1989) discuss a modification to the Henyey scheme proposed by Sugimoto (1970), which involves discretizing the thermal equations using asymmetric finite differences; my experience with older (pre-GYRE) codes has been that this works quite well in suppressing the noise. In his PhD thesis, Marc-Antoine Dupret discusses an alternative approach used in his MAD code, which also seems to work well.
For the multiple shooting scheme used by GYRE, the literature is unfortunately pretty sparse and I'm still looking for an efficient fix for the noise. In the meantime the noise can generally be eliminated by increasing the grid resolution where the thermal timescale is long. To do this, use the RESAMP_THERMAL operation when setting up the shooting grid:
&shoot_grid
op_type = 'CREATE_CLONE'
/
&shoot_grid
op_type = 'RESAMP_THERMAL'
alpha_thm = 1E-2
/
The key parameter here is alpha_thm (introduced in release 1.1 of the code); choosing larger values will increase the grid resolution, with the benefit of reduced entropy noise and the drawback of longer execution times.
Unfortunately, in the specific case of the input files you've sent me, it seems the number of grid points required for an accurate entropy eigenfunction is extremely large; here, the cure may be worse than the disease!
The upshot is: you can try playing around with alpha_thm (after upgrading to release 1.1), but I imagine a better fix will be necessary.
cheers,
Rich
Re: Non-adiabatic driving
Hi Rich,
First of all, I tried what you suggested. It helps with some of the noise in the deep interior but that's not really what I'm talking about. I think there must be a fundamental problem with the way dW_dx is being calculated, and also (possibly) with the non adiabatic results the code gives (i.e., the imaginary part of the frequency).
I attached some plots in a separate email to you but I don't know how to do that here. Unfortunately, they are quite central to what I'm talking about. I'm copying that email below so we'll have a record of this exchange.
I've gone and calculated the thermal timescale as a function of position for one of my models. This is plotted as the green curve in the lower panel of the plots I'm attaching, and the B-V frequency is plotted in blue. The red, dashed, vertical line denotes the point in the model where the thermal timescale equals the period of the mode (both in seconds).
The upper panels of all the plots are different, showing delta S/S, delta L/L, and dW_dx, respectively.
The first two plots, delta S/S and delta L/L, respectively, look just as I would expect: delta S/S is
large *only* above the adiabatic/nonadiabatic transition point, and delta L/L is essentially spatially
constant above this point (at least it is until it runs into the convection zone).
On the other hand dW_dx is large only below this point, in the adiabatic region. Deeper down it
certainly becomes noisy, but at the transition point and below it doesn't appear to be noisy at
all. So to me this doesn't look like a numerical problem. Rather, it looks like an error of some
sort in the computation of dW_dx. I also had the impression that by changing the theta_ad
parameter that I could change the value of the imaginary part of the frequency by factors of
several, so this problem may extend to the non adiabatic eigenfrequency calculation.
Of course, the error may somehow come from incorrect input from MESA via the GYRE output format.
I wouldn't be surprised at that at all. I've already found what I think is a small problem in "As". Using
the add_atmosphere option=True it seems that As=0 in the atmosphere when it should be non-zero.
This is why my B-V frequency, which is computed from As, dives down to zero at log dq = -15.
Anyway, please let me know if I haven't explained this clearly. I'd be happy to clarify.
Cheers,
Mike
First of all, I tried what you suggested. It helps with some of the noise in the deep interior but that's not really what I'm talking about. I think there must be a fundamental problem with the way dW_dx is being calculated, and also (possibly) with the non adiabatic results the code gives (i.e., the imaginary part of the frequency).
I attached some plots in a separate email to you but I don't know how to do that here. Unfortunately, they are quite central to what I'm talking about. I'm copying that email below so we'll have a record of this exchange.
I've gone and calculated the thermal timescale as a function of position for one of my models. This is plotted as the green curve in the lower panel of the plots I'm attaching, and the B-V frequency is plotted in blue. The red, dashed, vertical line denotes the point in the model where the thermal timescale equals the period of the mode (both in seconds).
The upper panels of all the plots are different, showing delta S/S, delta L/L, and dW_dx, respectively.
The first two plots, delta S/S and delta L/L, respectively, look just as I would expect: delta S/S is
large *only* above the adiabatic/nonadiabatic transition point, and delta L/L is essentially spatially
constant above this point (at least it is until it runs into the convection zone).
On the other hand dW_dx is large only below this point, in the adiabatic region. Deeper down it
certainly becomes noisy, but at the transition point and below it doesn't appear to be noisy at
all. So to me this doesn't look like a numerical problem. Rather, it looks like an error of some
sort in the computation of dW_dx. I also had the impression that by changing the theta_ad
parameter that I could change the value of the imaginary part of the frequency by factors of
several, so this problem may extend to the non adiabatic eigenfrequency calculation.
Of course, the error may somehow come from incorrect input from MESA via the GYRE output format.
I wouldn't be surprised at that at all. I've already found what I think is a small problem in "As". Using
the add_atmosphere option=True it seems that As=0 in the atmosphere when it should be non-zero.
This is why my B-V frequency, which is computed from As, dives down to zero at log dq = -15.
Anyway, please let me know if I haven't explained this clearly. I'd be happy to clarify.
Cheers,
Mike
-
rhtownsend
- Site Admin
- Posts: 397
- Joined: Sun Mar 31, 2013 4:22 pm
Re: Non-adiabatic driving
Hi Mike --
I've spent quite a bit of time looking into this issue -- hence the radio silence for the past week.
The calculation of the differential work appears to be fine; the problem really traces back to the entropy perturbation, which is not only noisy in the deep interior but also a couple of orders of magnitude too large. This in turn means that dW_dx is too large in the interior.
The root of the problem lies in dealing with almost-adiabatic equations; although the differential system is formally 6th order, it essentially behaves like a 4th-order system supplemented by a couple of algebraic equations. This makes things numerically difficult for the Magnus solver, hence the noisy entropy eigenfunction in the near-adiabatic interior.
I'm working on various approaches to fixing this; I'll be sure to let you know once I have something that works.
cheers,
Rich
I've spent quite a bit of time looking into this issue -- hence the radio silence for the past week.
The calculation of the differential work appears to be fine; the problem really traces back to the entropy perturbation, which is not only noisy in the deep interior but also a couple of orders of magnitude too large. This in turn means that dW_dx is too large in the interior.
The root of the problem lies in dealing with almost-adiabatic equations; although the differential system is formally 6th order, it essentially behaves like a 4th-order system supplemented by a couple of algebraic equations. This makes things numerically difficult for the Magnus solver, hence the noisy entropy eigenfunction in the near-adiabatic interior.
I'm working on various approaches to fixing this; I'll be sure to let you know once I have something that works.
cheers,
Rich
Re: Non-adiabatic driving
Hi Rich,
Thanks for all the great work on GYRE!
Using the most recent version of GYRE (that came with the MESA distribution) I still find problems with computing the local work in a strongly adiabatic environment (e.g., white dwarf interiors). I know this is a generic problem that all nonadiabatic pulsation codes face.
Is there any possibility of adding a feature for splitting a model into a fully adiabatic region (outer part of model) and a quasi-adiabatic region (inner part of model)? This is the solution that the Warsaw pulsation code uses and it seems to work well. I realize that this won't solve all problems in all situations (for instance, high-order g modes in the degenerate cores of evolved stars), but it would cover a considerable range of applications.
Thanks!
Mike
Thanks for all the great work on GYRE!
Using the most recent version of GYRE (that came with the MESA distribution) I still find problems with computing the local work in a strongly adiabatic environment (e.g., white dwarf interiors). I know this is a generic problem that all nonadiabatic pulsation codes face.
Is there any possibility of adding a feature for splitting a model into a fully adiabatic region (outer part of model) and a quasi-adiabatic region (inner part of model)? This is the solution that the Warsaw pulsation code uses and it seems to work well. I realize that this won't solve all problems in all situations (for instance, high-order g modes in the degenerate cores of evolved stars), but it would cover a considerable range of applications.
Thanks!
Mike
-
rhtownsend
- Site Admin
- Posts: 397
- Joined: Sun Mar 31, 2013 4:22 pm
Re: Non-adiabatic driving
Hi Mike --mikemon wrote:Hi Rich,
Thanks for all the great work on GYRE!
Using the most recent version of GYRE (that came with the MESA distribution) I still find problems with computing the local work in a strongly adiabatic environment (e.g., white dwarf interiors). I know this is a generic problem that all nonadiabatic pulsation codes face.
Is there any possibility of adding a feature for splitting a model into a fully adiabatic region (outer part of model) and a quasi-adiabatic region (inner part of model)? This is the solution that the Warsaw pulsation code uses and it seems to work well. I realize that this won't solve all problems in all situations (for instance, high-order g modes in the degenerate cores of evolved stars), but it would cover a considerable range of applications.
Thanks!
Mike
In fact, the functionality you suggest is already halfway implemented in GYRE. I'll see what I can do to get it finished, but right now I don't have too much GYRE time on my hands (teaching, tenure, etc).
cheers,
Rich