Posts for the month of March 2017

How cooling is called/used in the code

The code's main program is scrambler.f90, of course. In it are the initialization sequence calls, like 'physics init'. One of the things physics init does is call cooling init.

Cooling Init is in Cooling.f90. It is commented to say "initializes cooling table", but actually doesn't do anything… (It is a stub). Perhaps it is the relic from an older version of the code.

Now, the main meat of the cooling functionality goes like this.

In source_control.f90, there is a function called "SrcDerivs". This function's main goal is to populate the dqdt array with the code's source functions. Dqdt array contains geometrical source terms (like cylindrical corrections), energy source terms (like cooling), etc. At the end of this function, we have, "SrcDerivs=dqdt" (arrays).

If cooling is turned on, this function SrcDerivs calls the function "Cooling" (located in the cooling module: cooling.f90). "Cooling" returns a dqdt for SrcDerivs, specific for the type of cooling that is switched on.

For instance, take IICooling. "Cooling" selects the type of cooling that is present with a case statement:

If Case(IICool), Call II_Cooling

Now, IICooling is located in !IICooling.f90 and it calculates a dqdt(iE) based on a heating/cooling rate which themselves are functions of the II Cooling Parameters (stored in the array !IICoolPars). Once dqdt(iE) is calculated for this cooling choice, it is stored in the iE's element of the array SrcDerivs.

SrcDerivs array is used to update the q array in the subroutine, "RKorder2" back in source_control.f90. This is something like: q(i)=SrcDerivs(i)*dt + q(i).

The procedure would be similar (i.e. SrcDeriv and q would be updated in a similar manner) if another cooling process was selected instead of IIcooling.

Thus, it would seem, the only things necessary for turning cooling on in the code are the following two steps:

  1. Select the type of cooling in physics.data, e.g. iCooling=3 (for II cool) (as long as iCooling/=0, cooling should be turned on in the code)
  2. Define any cooling params to be used by cooling module by adding these to both the problemdata namelist in problem.f90 (e.g. IICoolPars() array for the case of II Cooling), as well as define their corresponding values in problem.data.

You should now be ready to cool!

Meeting Update --3/23/17

  • Wire Turbulence
    1. Schematic Diagram: http://www.pas.rochester.edu/~bliu/wireTurbulence/schematic3.png
    2. redid Spectra analysis with wave number range [2 20] comparing [2 40]. linear fit slope along x direction doesn't change much while y & z direction are goes much lower (-1.2). Details.

Meeting Update

Looking into tighter coupling between hyperbolic fluxes and source terms

https://astrobear.pas.rochester.edu/trac/wiki/u/johannjc/scratchpad24

And ways to improve corotating frame and possibly implement cylindrical geometry

https://astrobear.pas.rochester.edu/trac/blog/johannjc10222013

Update: Trying AMR, extending the profile, exploring fluctuations, including velocity-damping

I) Invoking AMR (with smooth and run time seconds)

Motivation:
Ultimately runs will be done with AMR so it is worth trying AMR at this stage and comparing with uniform grid runs.

Setup:
I introduced AMR with 5 levels of refinement for a run with dyne/cm, with refinement 'level 0' resolution equal to grid cells, in a box of physical dimension cm. This implies the equivalent of for level 5, which translates to a physical resolution equal (in principle) to that of Run (3) from the last blog post, which used grid cells for a box of dimension cm. So in a sense, this run combines the higher resolution of Run (3) with the larger box of Run (4). This run should in principle take about 10 days to complete on bluehive, but it did not actually complete (see below).

Results:
2d density slice full box 2d pressure slice full box 2d density slice central 2d pressure slice central 1d density central 1d pressure central

2d density fluctuations slice $6.6\times10^{-9}$ to $6.8\times10^{-9}$ g/cm$^3$ 2d pressure fluctuations slice $9.98\times10^{6}$ to $1.04\times10^{7}$ dyne/cm$^2$
1d density fluctuations $10^{-8.3}$ to $10^{-8.0}$ g/cm$^3$ 1d pressure fluctuations $10^{6.95}$ to $10^{7.1}$ dyne/cm$^2$

  • The run got killed, initially because the time allotment was exceeded. But after restarting the run, it again got killed, possibly because of a problem in the calculation (unclear). But clearly the run has become unphysical by this time anyway. 80 of 101 frames completed up to a time of s.
  • The pressure remains remarkably stable. The core pressure retains its initial value marginally better than Run (3), which had high resolution in a smaller box. Global oscillations are similar in amplitude and frequency to Run (4) (the larger box run), as would be expected, but seem to be less regular. This is best seen in the 1D pressure movie. This is not surprising given that the AMR grid is non-uniform, and refined patches are not spatially distributed in a regular manner.
  • The density, on the other hand, becomes unstable on a time scale of about 1 dynamical time. The star first develops small indentations at the points on the surface closest to the boundary and then develops a boxy morphology with the orientation aligned with the grid. Finally, the instabilities develop at the corners of the boxy stellar profile in the 2d slice.
  • The oscillations in the pressure stabilize somewhat after about half the simulation time (a few dynamical times). In the 1d plots that zoom in on the boundary, there seems to be at least two changes to the profiles of density and pressure: first at about 1 dynamical time ( s) and then again at a few dynamical times, when the density is seen to become unstable in the 2d density plot.

Discussion:

  • The fact that the density becomes unstable but the pressure distribution remains stable is consistent with what we found happens last blog post when the resolution is reduced (comparison plots of Run C). So it seems likely that somewhere in the box, there is a lack of resolution that leads to the instabilities in density. The instability seems to be strongest at the "corners" of the boxy profile that develops, which suggests that the resolution is too low to resolve these sharp corners.
  • It is worthwhile looking at the low-level fluctuations for the previous simulations as a comparison:

For the small box run we have (in 1d plots purple=smoothed pressure, green=unsmoothed pressure):
2d density slice boundary $6.6\times10^{-9}$ to $6.8\times10^{-9}$ g/cm$^3$ 2d pressure slice boundary $9.98\times10^{6}$ to $1.04\times10^{7}$ dyne/cm$^2$
1d density boundary $10^{-8.3}$ to $10^{-8.0}$ g/cm$^3$ 1d pressure boundary $10^{6.95}$ to $10^{7.1}$ dyne/cm$^2$

For the large box run we have (in 1d plots purple=smoothed pressure):
2d density slice boundary $6.6\times10^{-9}$ to $6.8\times10^{-9}$ g/cm$^3$ 2d pressure slice boundary $9.98\times10^{6}$ to $1.04\times10^{7}$ dyne/cm$^2$
1d density boundary $10^{-8.3}$ to $10^{-8.0}$ g/cm$^3$ 1d pressure boundary $10^{6.95}$ to $10^{7.1}$ dyne/cm$^2$

Conclusions:
The AMR fails to resolve the sharp density gradients that appear at certain points, causing the simulation to crash. Instead of developing 'clean' global pressure waves, like in the uniform grid simulations, the waves have a complex structure, with many small-scale fluctuations embedded within them. This probably helps to damp regular global osciallations that develop in the uniform grid models.

II) Adding an isothermal atmosphere (with smooth and run time seconds)

Motivation:
The profiles used in previous runs are characterized by a constant ambient density and pressure. The pressure profile was smoothed at the star-ambient boundary, but not the density profile. Adding a thick hydrostatic atmosphere for which the density and pressure decline with radius may be more realistic. In its simplest form, such an atmosphere would be isothermal, continuous in pressure, density and their respective gradients across the star-ambient boundary, and in hydrostatic equilibrium if self-gravity is neglected. The derivation of the profile, with plots, is available in the file profile_atmosphere.pdf.

Setup:
The atmosphere is added by simply modifying the density and pressure profiles inputted. It is ensured that these new profiles are defined over the entire box, so there is no need for any ambient values (either defined explicitly or imposed implicitly by astrobear).

Results:
Below are movies from two runs with the same physical resolution but different box size:

a) small box (with smooth and run time seconds)

2d density slice extended color bar 2d pressure slice extended color bar 2d density slice old color bar 2d pressure slice old color bar 1d density 1d pressure

b) large box (with smooth and run crashes after seconds)

2d density slice extended color bar 2d pressure slice extended color bar 2d density slice old color bar 2d pressure slice old color bar 1d density 1d pressure

  • Clearly, boundary effects lead to sharp gradients in density, which probably is what causes the code to crash in (b).
  • Morphologies appear unphysical, if one looks at small densities and pressures, but at the density/pressure range we were looking at before, the density slice remains quite circular and stable until a few dynamical times, when the density becomes suddenly unstable.

Conclusions:
These runs are somewhat less stable than the constant ambient density/pressure runs above, as sharp gradients develop more easily, probably as a result of reflection from the boundary. If this problem could be overcome, then this idea of an atmosphere may be more useful, because it seems to succeed (at least transiently) in transfering the disturbances to very low density and pressure where they are negligible.

New results (post-March23 meeting):
c) small box (with smooth , Multipole expansion BCs for Poisson solver, and run time seconds)

2d density slice extended color bar 2d pressure slice extended color bar 2d density slice old color bar 2d pressure slice old color bar 2d density slice extended color bar with velocity (scaled) 1d density 1d pressure

d) small box (with smooth , Multipole expansion BCs for Poisson solver, initial values in ghost zone held fixed, and run time seconds)

2d density slice extended color bar with velocity (scaled) 2d density slice extended color bar with velocity (unscaled) 2d density slice old color bar 1d density

III) Implementing damping (with smooth and run time seconds)

Motivation:
We wish for the star to be stable. As we have seen this is rather difficult to achieve. We must then resort to damping the motions by, in effect, adding a damping source term to the momentum equation. We follow the idea discussed by Ohlmann+17. That is to add a term of the form

,

where is a parameter.

Setup:
The most direct way to implement this form of damping would be to add the extra term directly to the momentum equation. However this would involve changing the code at the lowest level. For the time being we try an alternative prescription, whereby we set

,

with the timestep. Differentiating with respect to time, we obtain

,

since and hence are small. Alternatively, we can use

,

which leads to

.

However, it was unclear to me exactly how to implement this in the code, so I tried two different methods below.

Preliminary Results:

a) small box, with resolution , smooth and run time seconds, s, implementation method 1

2d density 2d pressure 1d density 1d pressure 2d hydrostatic equilibrium 2d sound speed 2d Mach number

b) small box, with resolution , smooth and run time seconds, s, implementation method 2
2d density 2d pressure 1d density 1d pressure 2d hydrostatic equilibrium 2d sound speed 2d Mach number

c) small box, with resolution , smooth and run time seconds, s, implementation method 2
2d density 2d pressure 1d density 1d pressure 2d hydrostatic equilibrium 2d sound speed 2d Mach number

…and the fiducial run for comparison (small box, , smooth , run time seconds, no damping):
2d density 2d pressure 1d density 1d pressure 2d hydrostatic equilibrium 2d sound speed 2d Mach number

Note that the fractional hydrostatic equilibrium

does not include the gravity of the point source (I must correct this).

d) small box, with resolution , smooth and run time seconds, s, implementation method 3
2d density 2d density and velocity 2d pressure 1d density 1d pressure

e) small box, with resolution , smooth and run time seconds, s, Multipole expansion BCs, implementation method 3
2d density 2d density and velocity 2d pressure 1d density 1d pressure

f) small box, with resolution , smooth and run time seconds, s, Multipole expansion BCs, implementation method 3
2d density 2d density and velocity 2d pressure 1d density 1d pressure

g) small box, with resolution , smooth and run time seconds, Multipole expansion BCs, NO DAMPING
2d density 2d density and velocity 2d pressure 1d density 1d pressure

Update 3/23

  • Finished correct WASP run - time extended by ~100000 years, so not terribly significant. Updating paper figures (still need to find good output format).

Graph Volume rendering

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

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 dyne/cm, 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 seconds (about 1 dynamical time) to seconds and then again to 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 to without changing the box size.
4) Increase the box dimension from cm (with ) to cm (with ).

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 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) seconds: 3d density 1d density 3d pressure 1d pressure

B) seconds: 3d density 1d density 3d pressure 1d pressure

C) 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 seconds, or one full oscillation. This is about 2.5 dynamical times.
  • Both the density and pressure eventually become completely unstable at about 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 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 dyne/cm pressure was to make . 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 dyne/cm. 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 and run time 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 , 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 run, so the oscillations are not caused by lack of resolution.
  • At the end of the simulation, the 512 run is fairly circular, and symmetric in the 1D plots, while the 256 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 run, though after the first 1-2 dynamical times, the morphology seems to be slightly more diamond-shaped and slightly less circular than the run (but eyeballing it isn't easy!).

4) Increasing the box size (with smooth and run time 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 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 with the larger box of 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 to 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 is continuous at the stellar surface (e.g. defined to be the radius at which dyne/cm), which will give the amplitude .
2) Assume an isothermal equation of state , with a constant. This constant is obtained by requiring that is also continuous at the stellar surface.

Alternatively, but less simply, one could perhaps retain the old profile for , and solve for (assuming hydrostatic equilibrium) assuming now that is a function of . This would then allow one to match both and 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.