Questions about outer boundary conditions / asymptotic behaviour of lag_P / eul_P

[vincent]

Dear Rich,
dear all,

while computing solar p-mode eigenfunctions and trying to understand the different outer boundary conditions (BCs), I came across a few difficulties.

I computed eigenfunctions for l=2 with gyre 5.1, JCD outer boundary conditions, and Model S, see the enclosed input file.

1.) The first question is whether it is possible to obtain eigenfunctions with the outer boundary condition of a vanishing Eulerian pressure perturbation? I am asking because I observe that eul_P = (P_1 /P) is steeply increasing at the outer boundary (where P_1 is the Eulerian pressure perturbation and eul_P is the output variable from gyre).

2.) After going through the description of the outer BCs implemented in gyre (in the gyre papers and equations.pdf), I am still not sure exactly which conditions from the ones discussed in Christensen-Dalsgaard (2008, "JCD") and Unno (1989, "UNNO") are actually implemented. My guess is that in the JCD condition, the Lagrangian pressure perturbation vanishes (first outer BC mentioned in JCDs paper). In the UNNO condition, I guess that in addition to zero Lagrangian pressure perturbation, also P=\rho=0 is required at the outer boundary (Unno 1989, p. 105, see text above eq. 14.10)? However, in addition, JCD (2008, page 3) discusses matching to an isothermal atmosphere, he also cites Unno (1989), and he mentions also a vacuum outer BC. In Unno (1989), there is also a discussion that does not require P = 0 and a discussion of some other options in section 21.2.

3.) I also observe that lag_P, which is \delta P / P caused by the p mode (\delta P being the Lagrangian pressure perturbation), takes its largest values at the outer boundary, where it steeply increases, see the enclosed figure. For this figure, I normalized the mode to a surface amplitude of 1cm. As I understood the JCD outer boundary condition to be \delta P = 0, I expected \delta P / P to go to zero, too, because P does not vanish in model S at the outer boundary.

What is the reason for this behaviour? Maybe there is a good reason for this which I am missing? Is it possible or does it make sense to request \delta P / P = 0 at the outer boundary? Of course, I am aware that 1/P strongly increases at the surface.

Looking at lag_P*P, which is Lagrangian \delta P caused by the p mode, it looks in principle fine, see the bottom panel in the enclosed figure.

In the plots, you find always two nearly identical curves. One is from the output lag_P from gyre, and the other one I computed from the displacement eigenfunctions to check. I used this way to check the behaviour for the modes that I have from adipls, too, and they look the same.

4) Finally, I tried the other outer BCs available, while everything else was as in the enclosed input file. UNNO gives the same answer as JCD within 1E-4 / 1E-5 in relative numbers. Is that because Model S does not vanish at the outer boundary? Running the code for VACUUM and DZIEM did not work for some reason, I posted a bug report in an additional thread.

Looking forward to any answers. Thanks in advance for any help or advice anyone can give!

Best regards,

Vincent Böning

[rhtownsend]

Dear Rich,
dear all,

while computing solar p-mode eigenfunctions and trying to understand the different outer boundary conditions (BCs), I came across a few difficulties.

I computed eigenfunctions for l=2 with gyre 5.1, JCD outer boundary conditions, and Model S, see the enclosed input file.

1.) The first question is whether it is possible to obtain eigenfunctions with the outer boundary condition of a vanishing Eulerian pressure perturbation? I am asking because I observe that eul_P = (P_1 /P) is steeply increasing at the outer boundary (where P_1 is the Eulerian pressure perturbation and eul_P is the output variable from gyre).

2.) After going through the description of the outer BCs implemented in gyre (in the gyre papers and equations.pdf), I am still not sure exactly which conditions from the ones discussed in Christensen-Dalsgaard (2008, "JCD") and Unno (1989, "UNNO") are actually implemented. My guess is that in the JCD condition, the Lagrangian pressure perturbation vanishes (first outer BC mentioned in JCDs paper). In the UNNO condition, I guess that in addition to zero Lagrangian pressure perturbation, also P=\rho=0 is required at the outer boundary (Unno 1989, p. 105, see text above eq. 14.10)? However, in addition, JCD (2008, page 3) discusses matching to an isothermal atmosphere, he also cites Unno (1989), and he mentions also a vacuum outer BC. In Unno (1989), there is also a discussion that does not require P = 0 and a discussion of some other options in section 21.2.

3.) I also observe that lag_P, which is \delta P / P caused by the p mode (\delta P being the Lagrangian pressure perturbation), takes its largest values at the outer boundary, where it steeply increases, see the enclosed figure. For this figure, I normalized the mode to a surface amplitude of 1cm. As I understood the JCD outer boundary condition to be \delta P = 0, I expected \delta P / P to go to zero, too, because P does not vanish in model S at the outer boundary.

What is the reason for this behaviour? Maybe there is a good reason for this which I am missing? Is it possible or does it make sense to request \delta P / P = 0 at the outer boundary? Of course, I am aware that 1/P strongly increases at the surface.

Looking at lag_P*P, which is Lagrangian \delta P caused by the p mode, it looks in principle fine, see the bottom panel in the enclosed figure.

In the plots, you find always two nearly identical curves. One is from the output lag_P from gyre, and the other one I computed from the displacement eigenfunctions to check. I used this way to check the behaviour for the modes that I have from adipls, too, and they look the same.

4) Finally, I tried the other outer BCs available, while everything else was as in the enclosed input file. UNNO gives the same answer as JCD within 1E-4 / 1E-5 in relative numbers. Is that because Model S does not vanish at the outer boundary? Running the code for VACUUM and DZIEM did not work for some reason, I posted a bug report in an additional thread.

Looking forward to any answers. Thanks in advance for any help or advice anyone can give!

Best regards,

Vincent Böning

Hi Vincent --

Thanks for your post -- I'll see if I can help answer:

1) That would require a modification of the outer boundary condition code in GYRE. I can certainly show you how to do this, but it will require a bit of work on my part; so I'll respond with a later post about this topic.

2) The 'UNNO' boundary condition in GYRE, which is described in Unno et al. (1989), does not make any assumptions about the behavior of the surface pressure, or delta P; rather, it performs a local dispersion analysis at the stellar surface to find the local wave solutions: either upward- and downward-propagating waves (if we're above the acoustic-wave cutoff frequency, or below the gravity-wave cutoff frequency), or upward- and downward-decaying evanescent solutions. It then applies the constraint that there is no downward-propagating wave (if we're above/below the cutoffs), or that the energy density decays outward.

I don't have a copy of Unno et al. (1989) with me right now (I'm at the TASC4/KASC11 conference in Aarhus), but you might want to look at Townsend (2000, MNRAS, 319, 289) to see how I go about choosing the correct solution.

The 'JCD' boundary condition in GYRE is basically a cut-down version of the 'UNNO' condition, made under the assumption that the atmosphere is isothermal. A full description of the BC is given in section 2.4.3 of Jørgen Christensen-Dalsgaard's user guide for his ADIPLS code; you can grab a copy of this guide from http://www.astro.wisc.edu/~townsend/res ... -guide.pdf

3) The behavior you see in delta P/P is because the JCD condition does *not* enforce delta P/P -> 0 at the surface

4) Because the atmosphere in model S *is* close to isothermal, there's almost no difference between the JCD and UNNO boundary conditions. That's why they produce similar results.

For the DZIEM boundary condition, the crash occurs because in your GYRE inlist file, you've specified the upper search frequency (freq_max) in units of the acoustic cutoff frequency. This cutoff frequency is formally infinite in the DZIEM BC (there *is* no upper cutoff), and this causes GYRE to crash. In a future release, I'll make the code behave more gracefully; but in the meantime, you can get the DZIEM BC working by chosing something other than 'ACOUSTIC_CUTOFF' for the freq_max_units parameter.

For the VACUUM boundary condition, there is a bug in the code which stops it from properly recognizing that 'VACUUM' is a valid choice. I'll fix this in the next release. However, even if 'VACUUM' worked, it would run into the same problem as 'DZIEM'. Therefore, the workaround is the same as in the 'DZIEM' case.

Hope this answers your questions!

cheers,

Rich

[vincent]

Dear Rich,

thanks a lot for your detailed answers that make sense and thanks for pointing me to the right locations in the references! Also the issue with the VACUUM and DZIEM BCs makes very much sense now and they give nice results now.

Thanks also for your offer to help modify the code:

1) That would require a modification of the outer boundary condition code in GYRE. I can certainly show you how to do this, but it will require a bit of work on my part; so I'll respond with a later post about this topic.

I think that I can live without this option for now, so I would prefer not to bother you with that work. I would get back to you if it got relevant for my work!

Greetings,

Vincent