Mapping a modified RGB profile to the grid: Increasing time, resolution, or box size, and smoothing P at the surface

Recall from the last post that the RGB star was most stable for case (D), which was the case where the ambient pressure is set to P_\mathrm{ambient}=10^7 dyne/cm^2, and the ambient density parameter is set to a very low value such that the actual ambient density is set automatically by AstroBear. I now extend the study in four ways:
1) Increase the run time from t=4\times10^5 seconds (about 1 dynamical time) to t=2\times10^6 seconds and then again to t=8.52\times10^6 seconds.
2) Make the pressure profile smooth at the transition from the MESA profile to the ambient value, as the discontinuity in the pressure gradient could be causing unwanted effects.
3) Increase the resolution from 256^3 to 512^3 without changing the box size.
4) Increase the box dimension from 10^{13} cm (with 256^3) to 2\times10^{13} cm (with 512^3).

Below I present the results of each of these experiments in turn. 2D plots are slices through the center of the Y-axis, while 1D plots show an extra slice through the center of the Z-axis.

1) Increasing the run time

Model (D) of last blog post ran for 4\times10^5 seconds, and is labeled (A) below. Model (B) below had the same parameter values but ran for 5 times as long, with the same number of snapshots (101). Model © below was meant to run for 25 times as long as the original run and 5 times as long as run (B), but stopped just short of completion due to bluehive crashing (presumably not because of this run!). It contains 214 snapshots.

A) t=4\times10^5 seconds: ​3d density ​1d density ​3d pressure ​1d pressure

B) t=2\times10^6 seconds: ​3d density ​1d density ​3d pressure ​1d pressure

C) t=8.52\times10^6 seconds: ​3d density ​1d density ​3d pressure ​1d pressure

Discussion:

• The core density and pressure decrease somewhat but then remain fairly stable.
  
• However, the density and pressure near the surface oscillate with time.
  
• This appears to be caused by a pressure wave which starts near the surface, moves inward, reflects off the core, moves outward, reflects off the box boundaries, and moves inward again toward the core. This is most evident from the 1D plots of pressure for the longer run times.
  
• The four-fold symmetry of the grid becomes apparent in the 3d density profile by about 1\times10^6 seconds, or one full oscillation. This is about 2.5 dynamical times.
  
• Both the density and pressure eventually become completely unstable at about 4\times10^6 seconds, or about 10 dynamical times.

Conclusions:
The star becomes unstable after about 10 dynamical times. Even before this, it shows oscillations during which the density profile experiences kinking (cuspiness). The oscillations appear to be caused by pressure waves reflecting off the core and grid boundaries. Below we explore three possible ways to help reduce these oscillations and cuspiness: smoothing the pressure profile near the surface, increasing the resolution, and using a larger box.

2) Smoothing the pressure profile near the stellar surface (with run time t=2\times10^6 seconds)

I realized that the simplest way to smooth the pressure profile to avoid a discontinuous pressure gradient while going from the surface to the ambient 10^7 dyne/cm^2 pressure was to make P\rightarrow P+P_\mathrm{ambient}. Because the pressure spans several orders of magnitude, this prescription hardly affects the pressure near the core, and only begins to be important about 5-10 solar radii from the `surface,' where surface here means the radius at which P=10^7 dyne/cm^2. However, this prescription changes the pressure and pressure gradient near the surface, so although it avoids a sudden change in pressure gradient, hydrostatic equilibrium is not expected to be satisfied quite as accurately, at least not initially. This run took about 6 hours to complete on bluehive.

​3d density ​1d density ​3d pressure ​1d pressure

​3d density P non-smooth(left) vs smooth(right)
​1d density P non-smooth(green) vs smooth(purple)
​3d pressure P non-smooth(left) vs smooth(right)
​1d pressure P non-smooth(green) vs smooth(purple)

Discussion:

• Four-fold symmetry in the density still appears but at somewhat later times, so the star is rounder than in the unsmoothed case.
  
• The star is somewhat larger than in the unsmoothed case.
  
• The core behaves similarly to the unsmoothed case.
  
• The frequency of oscillations is reduced compared to the unsmoothed case.
  
• The cuspiness produced during the oscillations is reduced compared to the unsmoothed case.
  
• The amplitude of oscillations remains about the same as in the unsmoothed case.
  
• The star oscillates between a state similar to the initial state and a state where it is larger; i.e. it pulsates.

Conclusions:
Smoothing the profile helps to avoid cuspiness and reduces the frequency of oscillations. However, with increased pressure in the outer regions, the star expands to be greater than its original size, before contracting again to its original size and undergoing fairly regular oscillations between these states. Overall, the tradeoff between less rapid less cuspy oscillations and slightly larger star is probably worth making, so (for now at least) we adopt this smooth pressure profile prescription in the runs described below.

3) Increasing the resolution (with smooth P and run time t=2\times10^6 seconds)

A natural step is to increase the resolution, as some of the effects discussed above may be caused by lack of resolution of the pressure scale height, either near the surface or near the core. Therefore, we next doubled the resolution to 512^3, keeping the box size constant, and retaining the smoothed pressure profile. This run took about 3.5-4 days to complete on bluehive.

​3d density ​1d density ​3d pressure ​1d pressure

​3d density-pressure comparison ​1d density-pressure comparison

​3d density 256$^3$ (left) vs 512$^3$ (right)
​1d density 256$^3$ (left) vs 512$^3$ (right)
​3d pressure 256$^3$ (left) vs 512$^3$ (right)
​1d pressure 256$^3$ (left) vs 512$^3$ (right)

Discussion:

• The initial central density and pressure are maintained much more faithfully in the high resolution run (though they still decrease and then quasi-stabilize at slightly lower values than the initial values). This is as expected because the relatively small pressure scale height near the center is better resolved.
• The amplitude, frequency and morphology of the oscillations is very similar to the 256^3 run, so the oscillations are not caused by lack of resolution.
• At the end of the simulation, the 512^3 run is fairly circular, and symmetric in the 1D plots, while the 256^3 run is boxier and shows asymmetry in the 1D plots. The higher resolution thus seems to make the solution more stable. However, the 3D density profile shows that the star seems to become slightly `diamond' shaped closer to the beginning of the run, and then recovering to a more circular morphology.

Conclusions:
Increasing the resolution makes the star remain stable for a longer time (though our run was not long enough to know when the star becomes unstable, if it becomes unstable). The core density and pressure are better preserved (to within about 15% instead of to within about 40%). The star is clearly more circular by the end of the run compared with the 256^3 run, though after the first 1-2 dynamical times, the morphology seems to be slightly more diamond-shaped and slightly less circular than the 256^3 run (but eyeballing it isn't easy!).

4) Increasing the box size (with smooth P and run time t=2\times10^6 seconds)

Since the pressure wave discussed above seems to reflect off the boundaries of the box, perhaps such waves could be reduced, in frequency and/or amplitude, by pushing out the boundaries. Therefore, we now make the grid dimension twice as large. This run took about 3 days to complete on bluehive.

​3d density ​1d density ​3d pressure ​1d pressure

​3d density 256$^3$ with $L=1\times10^{13}$ cm (left) vs 512$^3$ with $L=2\times10^{13}$ cm (right)
​1d density 256$^3$ with $L=1\times10^{13}$ cm (left) vs 512$^3$ with $L=2\times10^{13}$ cm (right)
​3d pressure 256$^3$ with $L=1\times10^{13}$ cm (left) vs 512$^3$ with $L=2\times10^{13}$ cm (right)
​1d pressure 256$^3$ with $L=1\times10^{13}$ cm (left) vs 512$^3$ with $L=2\times10^{13}$ cm (right)

Discussion:

• The frequency of the oscillations is clearly reduced, and just less than one full oscillation is completed before the end of the run (I use the word `oscillation' assuming that it will turn out to be periodic if run for longer). For the standard box size we were getting almost two full oscillations by the end of the run. So the frequency is reduced by about half. This is evidence that these waves are due to reflections off the walls.
• Moreover, the amplitude of the oscillations is drastically reduced.
• However, the star morphology is slightly `boxier' in the large grid run. This boxiness becomes apparent after about two dynamical times or about 8\times10^5 seconds. This could perhaps be due to an extra time lag between reflections from the wall and expected reflections from the corner of the grid.
• By the end of the run there is less asymmetry in the 1D density and pressure profiles in the larger box run, which seems to indicate that the star remains stable for longer.
• The decrease and quasi-stabilization of the core density or pressure are similar to the smaller box case.

Conclusions:
Increasing the grid size helps to reduce the oscillations, both in terms of frequency and especially in terms of amplitude. However, the density profile still becomes boxy, perhaps even boxier than with the smaller grid. The oscillations and boxiness are probably caused by pressure waves reflecting off the walls.

Overall discussion and conclusions:

• I adopted three measures designed to increase the stability of the star: smoothing of the pressure profile, increase of the resolution, and increase of the grid size. Smoothing the profile leads to oscillations that are less cuspy, and a density profile that is less boxy, but the star inflates more significantly during the oscillations. The core is hardly affected.
• Increasing the resolution does not change the frequency or amplitude of oscillations, but the star remains stable for longer. It also retains its round shape better though it seems to develop a diamond-shaped morphology. The core density and pressure are better maintained, likely because the scale height at the center is better resolved.
• Increasing the grid size dramatically reduces the amplitude of oscillations and also reduces their frequency approximately in proportion to the increase in grid size (doubling the grid seems to halve the frequency). However, the star still becomes boxy, and slightly boxier than for the smaller grid.
• The smoothing has already been incorporated into the runs with increased resolution or box size. Using 1024^3 with the larger box of 2\times10^{13} cm is too demanding computationally. Therefore it is best to switch to AMR at this point.

Next steps:
I) AMR
I am now running a case with 5 levels of AMR corresponding to 32^3 to 1024^3 resolution of a fixed grid. The simulation will take 9-10 days to complete in total, or about 2-3 more days from today. I will be presenting the results here after it finishes.

II) More realistic atmosphere
The strong pressure waves possibly rely on the ambient pressure being as large at the grid boundaries as it is at the stellar surface. This high ambient pressure is rather unrealistic anyway. What would be more realistic is that rather than assuming a uniform pressure ambient medium, to instead insert an extended isothermal atmosphere for which pressure and density decay with distance. The pressure and density profiles can be set by solving the equation of hydrostatic equilibrium in spherical symmetry neglecting self-gravity, which should anyway be very small for this atmosphere. The solution could be subjected to the following constraints:
1) Assume P is continuous at the stellar surface (e.g. defined to be the radius at which P=10^7 dyne/cm^2), which will give the amplitude P_0.
2) Assume an isothermal equation of state P=c^2\rho, with c^2 a constant. This constant is obtained by requiring that \rho is also continuous at the stellar surface.

Alternatively, but less simply, one could perhaps retain the old profile for \rho, and solve for P (assuming hydrostatic equilibrium) assuming now that c^2 is a function of r. This would then allow one to match both P and dP/dr at the surface.

I feel it is worth trying such atmospheres before resorting to implementing an artificial forcing term to damp motions, as done in other studies such as Ohlman et al.