Posts for the month of January 2017

Plan for next step of mapping the modified RGB profile to the grid

Now that we have been able to generate a modified MESA profile, we must map it to the AstroBear grid. The steps are as follows:

1) Determine the desired cutoff radius and resolution.
2) Solve the modified Lane-Emden equation as explained in the last blog entry, including an iteration over the mass of the gravitation-only particle.
3) Produce an input file for AstroBear with tabulated , , and values.
4) Perform multiple runs for comparison.
5) Assess the level of hydrostatic equilibrium for each run.
Let us consider item (1). For now let us consider a uniform grid. Let the resolution be . The outer radius of the star is , and the cutoff radius is . As in Ohlmann, the softening length for the spline potential of the gravitation-only particle is also . (For the MESA solution will be replaced by the modified Lane-Emden solution for , and for will be determined by re-integrating the hydrostatic equilibrium equation inward using , with equal to the MESA value.) Further, we define , that is, the ratio fo cutoff radius to outer stellar radius. The number of resolution elements over a softening length is given by , where is the grid element size. Finally, let the box size be .

Putting this together, we find

This formula gives the resolution required for a given , , and . Ohlmann+16a states "…we find that a resolution of about 10 cells per softening length is required to ensure energy conservation during the in-spiral." Also, "…in a sphere of five softening lengths of the gravitation-only particle, the maximum cell radius was bound to a tenth of the softening length." This suggests using , so we set . For the RGB star, . To be consistent with our previous simulation, we may choose . As in Ohlmann+16c, we can try , , or . This results in required resolutions , , or , respectively, or, in physical units, , , or . Ohlmann+16a states "The smallest cells near the RG core have a radius of about at the beginning and about at the end of the simulation."

We now consider item (4). It would be useful, initially, to compare the following runs:
a) , .
b) , (i.e. lower resolution).
c) , (i.e. even lower resolution).
d) , , without mass iteration step in item (2), so using a slightly larger particle mass.
e) , , with mass iteration step but without modifying the MESA profile (i.e. skipping item (2) but using the particle mass from model (a)).
f) , (i.e. smaller cutoff radius).
g) , with ambient pressure set to higher value to effectively cut off very outer layer of star where scale height cannot be adequately resolved.

Resolving the scale height

One important additional point is that one should adequately resolve not only the softening length , but also the local pressure scale height . See page 3 of the past presentation rgb.pdf for a plot of scale height vs. radius. Modifying the MESA profile for will lead to larger scale heights there, but we must still check whether is large enough. This will likely be most problematic at the stellar surface , where . Here we will not be resolving the scale height so we would expect velocity perturbations to arise. Ohlmann+16c calculates a necessary but not sufficient expression for to ensure that Mach number fluctuations are below a specified level. Assuming the Mach number at , after the first part of a time step (see their Sect. 2.3 and Appendix A for details) one obtains

,

where is the Courant-Friedrichs-Levy constant and is the adiabatic index, so that

.

E.g. for , , and , we obtain . However, Ohlmann+16c somehow obtain the slightly larger value . In any case from this condition we require , and possibly , to avoid Mach number fluctuations greater than . However Ohlmann+16c comments "Even if the resolution requirement is met, other sources of numerical error (e.g. interpolation error, errors from the gravity solver) still introduce spurious velocity fluctuations. Thus, an appropriate relaxation procedure is necessary when mapping stellar models to hydrodynamical grids: the velocity fluctuations have to be damped."

It is important to replot the scale height vs radius plot of rgb.pdf for the modified RGB profile, as well as for the case of a higher ambient pressure.

An updated README file with instructions on how to run the code, along with the relevant files, is available here.

Professor De Marco's project and orbital decay

10 AU model will see separating binary stars and others will have orbital decay. 3 AU model will shrink 10% of the separation in 4000 yrs.

Thomas, Professor De Marco's student, communicated with me before. Jonathan and I took a look at their data this morning and find that the effective resolution is 32000. The size of the simulation box is and the finest cell is . It need at least 6 levels of AMR.

This job is very computationally intense.

Meeting Update --01/26/17

  • Contact User of Toronto?
  • Book Vista for the new semester?
  • Wire Turbulence
    1. Al Cooling: tried new parameters based on the aluminum table (density and temperature range). 2D runs show that the cooling intensity is too small comparing with the post shock energy. So the cooling is too slow comparing with the downstream velocity. Details can be found on this page.
    2. Tracers and Gaussian2 fit: Added tracers for grid and wind material. Gaussian 1 won't do a good fit for the PDF of either grid or wind material only. Details can be found on this page
http://www.pas.rochester.edu/~bliu/wireTurbulence/Tracers/logWG0200.png http://www.pas.rochester.edu/~bliu/wireTurbulence/Tracers/GridWind.png
  1. Paper: reading and writing..

Update 1/26

  • Working on (done with?) application for Blue Waters fellowship
  • WASP-12b:

Ran for 10 more orbits - planetary mass loss rate not asymptoting yet, but it does appear to be becoming more variable (just numerical?).

Qualitatively, disk is identical:

Also running with

  • MHD: Running with betas for comparison with Owen and Adams - time growing exponentially again for stronger B, so need to look at that.

Modified giant profiles (continued from last blog)

I continued to work on generating the modified RGB profile to use in AstroBear, as discussed in the previous blog.

As I had been confused by the mismatch between the modified Lane-Emden (MLE) solution and the MESA solution at the transition radius , I wrote S. Ohlmann to ask about this. He explained that the pressure profile of the MLE solution is not actually used. Instead the equations for and are integrated again, this time backwards starting from some point on the MESA profile, to obtain and for . Thus, only the MLE profile for is used. (So the final profile for will not be a true polytrope, but that is okay.) This made sense but then at least the MLE profile for should almost match the MESA profile at , I realized, because both solutions satisfy hydrostatic equilibrium, after all, and by construction is almost the same for both solutions. The reason it is almost the same and not exactly the same will be discussed below.

This led me to look at my code and sure enough I discovered that the solution was not properly converging. It was oscillating between two apparent solutions, which led to the error in the profile. Therefore I set on fixing the PDE solver. The method I had been using was a shooting method (of my own design). I consulted Numerical Recipes (Press et al.) and found the relevant programs, but eventually realized I would have to start from scratch if I did it that way. I modified my own program and it now converges properly, iterating over 2 free parameters to satisfy the two boundary conditions at r=h. (It could instead be solved as a boundary value problem but that is more complicated and unnecessary.) A useful "trick" is to first run the code in low resolution, which converges in a few seconds, and then you have a good estimate of the parameters which you can use to run in high resolution.

As can be seen from this new version of profile.pdf, the profiles are reasonable and match those of Ohlmann. The MLE profile for does not precisely match MESA at (as explained last blog there is no reason it should) but is now only off by about 20% for the polytrope solution with . This seems very reasonable. The only concern I had (which I kept revisiting) is that obtained from the hydrostatic equilibrium equation by first integrating to obtain , does not perfectly match MESA at the boundary. After studying various solutions (with or or ) I believe this is just because the IDL differentiation and smoothing functions used to obtain from the MESA profile is imprecise and sensitive to the number of points chosen for the smoothing, etc. I don't think this will cause problems in AstroBear. Note that the profile matches perfectly at the boundary. But this is by design: IDL was used to calculate which was then used to set the boundary condition that should equal the MESA value. Note that MESA apparently does not give the option of directly outputting .

There is, however, one slight inconsistency in the method is the following (as alluded to above). The value of of the MESA profile is equal to the point particle mass . However, after solving the MLE, we have , where is the gas profile obtained by solving the MLE. By integrating backward and resetting , a discontinuity is avoided. But this means that the core at is being replaced in AstroBear by a mass that is greater than the MESA value of . In other words, we have added the extra mass to the core, so the MESA profile for will no longer be in perfect hydrostatic equilibrium. I discussed this issue with Ohlmann who was aware of it. He says he tried setting the particle mass to some value less than and then iterating over (yet a third parameter over which to iterate). He said it did not make enough of a difference to worry about. Tomorrow morning I will try a manual iteration and report back here.

As promised I have now tried the iteration procedure mentioned. This was done by modifying the code to include an additional parameter: which is not equal to . I obtain in the usual way from MESA but then instead of setting to equal , I set , with a fraction less than (but close to) unity. I first set to a reasonable value 0.96 given that for , we had . After a few manual iterations, I was able to get a converged solution such that . The value of turns out to be 0.963 for our fiducial profile ( polytropic index and ). I have now added the plots to the end of profile.pdf (Figs. 14-16). I have compared these plots against those where the iteration over is not performed (Figs. 2-5). The profiles are almost indistinguishable. However, the new method including this manual iteration is more self-consistent. It is not difficult to implement nor does it require significantly extra time. Therefore I see no reason why it should not be used when constructing the profiles to be used in AstroBear. Including this last step, not included by Ohlmann, will probably not typically make an important difference. Nevertheless, it is an improvement that I think is worth noting in the paper.

Ionization Test Problem

WASP-12 column density

Need 2x1017 /cm2 column density of magnesium, so ~2x1022 of hydrogen. Currently too small - by calculating average density in a box around the disk and multiplying by length of box, I find ~2x1018 /cm2 of hydrogen. Would need to increase density by 4 orders of magnitude - plotting average density in this box over time, estimate 105 orbits, or about 300 years, to build up sufficient density.

They estimate that a mass loss rate of 3x1010 g/s for WASP-12b will give an upper limit of 4x1021 /cm2 for Mg column density.

At 3.3 Rp, mass loss rate of planet is around 8x109 g/s (appears to be decreasing over time, though), which is about 3.75 times lower than the mass loss rate used to estimate the column density - linear in mass loss rate, so this still gives an estimate of 1021 /cm2. A rough estimate of the mass loss rate of the disk gives 5.1x107 g/s, so two orders of magnitude smaller than the mass loss rate of the planet - which correlates well with an increasing overall disk mass (using slope of disk mass over time, get ~2.4x109 g/s as a very rough estimate - within an order of magnitude).

De Marco's project

Fluid Code

Meeting Update --01/06/17

  • Visitors and Users
    1. UCSB visitor next week: volunteers needed
    2. RIT user: compiling issues with *.cpp.f90 files with gnu fortran 4.8?. Possible problems for our visitor and other users
    3. Bruce: OH231 module and Hen3–401 paper
    4. Toronto user requesting for Binary code: NK cooling table.
  • XSEDE allocations
    1. 2M CPU-hours on XSEDE machines (1M on Stampede and 1M on Comet) available for production runs. Current allocations can be found here
    2. Comparison of the Stampede and Comet machines can be found on this page
  • Wire Turbulence (
    1. Gaussian 2 fit for density PDF: using wire and wind materials using temperature. Details can be found here.
http://www.pas.rochester.edu/~bliu/wireTurbulence/seperatePDF/Hist_2X4.png
  1. Tracers (to do)
  2. Al Cooling (testing)
  3. Paper (working on)

Modifying stellar profiles from MESA for input into AstroBear

I continued working toward getting a giant (RGB or AGB) onto the grid.

Recall that when trying to translate the MESA stellar profile to the grid, the compact core gets washed out due to the grid not being fine enough to adequately resolve the local scale height. Last time, I showed that adding a sink particle with the mass of the unresolved RGB core can help to slow down the (unphysical) expansion of the star. But the star is still unstable. This is at least partly due to the fact that the inputted stellar profile has not been modified to take into account the modification to the potential due to the replacement of the central part of the star by a point particle. As explained in Ohlmann 2016 (phd thesis), the asymptotic behaviour of the pressure gradient as tends to 0 must be correctly reproduced ( as ) to avoid unphysical sound waves eminating from the center. This is achieved in Ohlmann 2016 by truncating the MESA profile at some radius (the higher the resolution, the smaller may be chosen to be) and replacing the inner part of the profile with the solution to a MODIFIED Lane-Emden equation that takes into account not only the potential of the gas (as in the standard LE equation), but also that of the central point particle. Doing this involves a few steps:

1) Choose a truncation radius and obtain the interior mass from the MESA profile. This is the mass of the central point particle .

2) Modify the potential to account for the central point particle using a softening to prevent it from blowing up as . Here we follow Ohlmann 2016 (see also Ohlmann et al. 2016, arXiv:1612.00008).

3) For simplicity, approximate the stellar profile for as a polytrope , with the polytropic index a parameter that Ohlmann sets to . The pressure , with constant, so for , we obtain , as for a relativistic white dwarf. Assuming mass continuity and hydrostatic equilibrium leads to the modified Lane-Emden equation derived in Ohlmann.

4) Now solve the modified Lane-Emden (MLE) equation by first breaking the equation up into a system of two coupled first order PDEs and then solving the system numerically. I have done this using a 3rd order Runge-Kutta "time"-stepping routine in fortran (Brandenburg 2003). For the simplest case of , the MLE equation reduces to the standard LE equation. I've tested my code for this case against a simpler, 1st order code downloaded from the web.

Solving the MLE equation is actually much more challenging than solving the LE equation. This is because aside from , there are now two additional parameters: the central gas density , and the characteristic radius ( can be written in terms of ). Basically, these parameters must be iterated over until a solution is found which satisfies the BOUNDARY CONDITIONS at the transition radius . (Ohlmann mentions using a non-linear root-finder to accomplish this step, but does not provide details.) I have used a simple "home-built" recipe for the iteration that leads to a converged (rapidly enough) solution.

Following Ohlmann, we use the following BCs:
i) matched to MESA profile
ii) matched to MESA profile

6) This approximate (polytrope) solution for the modified profile for can then be combined with the MESA profile, so that for , one retains the more realistic MESA profile. Now the modified profile can be compared with that of Ohlmann.

Refer to the file profile_blog1.pdf.

Figure 1 is the relevant figure from Ohlmann et al. 2016 with which results may be compared.

Figure 2 shows the mass profile of the RGB star as obtained from MESA, and the chosen transition radius , where is the outermost radius of the star. The central mass is about , in agreement with Ohlmann.

Figure 3a shows the density obtained from MESA (thick blue) and from the MLE equation (thin black). We have adopted and the polytropic index .
Figure 3b shows the density gradient . Note that and of the MLE solution are equal to the MESA values at . In other words, the boundary conditions have been correctly implemented. Note also that the MLE solution for agrees with Ohlmann's solution from Figure 1.

Figure 4a shows the pressure . Note that the MLE solution for (black line) DOES NOT AGREE with the MESA solution at . However, in Figure 1, Ohlmann obtained agreement. Have I made a mistake? Actually, there is no reason to expect perfect agreement, since the MESA solution is not an polytrope; MESA employs a more complicated equation of state known as OPAL. However, we can RE-NORMALIZE the MLE solution for to make it match the MESA solution at (red line). This re-normalized solution now agrees with Ohlmann's solution of Figure 1. Ohlmann mentions "integrating" the stellar structure equations to obtain , presumably while using the MLE profile for . But the relation between and has already been specified when solving the MLE equation, with the constant (or ) set to satisfy the BC on smoothness of across the transition radius. Thus, after re-normalizing, hydrostatic balance will no longer be perfect. On the other hand, the pressure is now continuous across the transition radius.
Figure 4b shows the pressure gradient . Again, we must re-normalize at . Doing so brings the solution into close agreement with that of Ohlmann, seen in Figure 1. The re-normalization factor is close to but not precisely the same as for .

Figures 5 and 6 show the profiles for , , , and for the case , while Figures 7 and 8 show the same for . The solutions obtained are quite similar to that obtained for the case.

Figures 9, 10, and 11 show the case . That is, with double the transition radius (and ). The MLE solutions can be compared with the green lines in Figure 1 and show good agreement.

Figures 12 and 13 show the fiducial and , but now with the boundary conditions modified so that the solution to the MLE equation has equal to the MESA value (as before) and equal to the MESA value (instead of ). We see that this results in a solution for which and do not match the MESA solution at the transition radius , as would be expected.

Most of the relevant fortran, IDL and tex files for the above work including a README file can be found here.

Remarks:
Matching all quantities , , , and to the MESA solution at AND satisfying hydrostatic equilibrium is not possible to do when approximating the profile for as a polytrope. This is because there are only 2 free parameters ( and ), whereas there are 4 boundary conditions, so the problem is overdetermined. Thus, either smoothness or hydrostatic equilibrium must be sacrificed. To the best of my understanding at least, Ohlmann et al. 2016 chose to sacrifice perfect hydrostatic equilibrium to make and smooth accross the transition radius.

Some natural next steps would be:
1) Plot the quantity for the modified profile to assess the level of agreement with hydrostatic equilibrium. Do the same for the unmodified MESA profile and for the pure polytrope solution of the MLE equation.
2) Map the modified solution to the AstroBear grid and perform a run with the appropriate sink particle mass () to assess the level of stability.