Summary of last week's work

1. Worked on alpha prescription stuff for energy paper and discussed with Eric.
2. Continued running simulation with secondary mass halved (run 149) ​separation vs time plot
3. Started running new simulations to explore effects of softening length and resolution.
   • This consisted of 3 new simulations restarted from frame 72 (t=16.7 days) of run 143 (Model A of Paper I), where the softening length and smallest resolution element had been halved.
     1. do not halve softening length nor double resolution (Run 152)
     2. halve softening length but but do not double resolution (Run 153)
     3. do not halve softening length but double resolution (Run 154)
4. Discussed jets with Jonathan.
5. Ran low resolution simulations to explore the possibility of simulating Roche lobe overflow (leading up to common envelope evolution).

Energy paper

Ivanova+2013 equation 3:
\frac{\mathrm{G} M_1 M_\mathrm{1,e}}{\lambda R_1} = \alpha_\mathrm{CE} \frac{\mathrm{G} M_\mathrm{1,c} M_2}{2}\left(\frac{1}{a_\mathrm{f}} -\frac{M_1/M_\mathrm{1,c}}{a_\mathrm{i}}\right)

New equation suggested by us:
\mathrm{G} M_\mathrm{1,e}\left(\frac{M_1}{\lambda R_1} +\frac{M_2}{2a_\mathrm{i}}\right) = \alpha_\mathrm{CE}\frac{\mathrm{G} M_\mathrm{1,c} M_2}{2} \left[\frac{1}{a_\mathrm{f}} -\frac{1}{a_\mathrm{i}}\left(2-\frac{M_1^2/M_\mathrm{1,c}+M_2}{M_1+M_2}\right)\right]

Run 143 (Model A of Paper I) has initial separation a_\mathrm{i}=49R_\odot. The Roche limit radius for q=2 as in our case is 109 R_\odot.

With a_\mathrm{final}=7 R_\odot, which is the final separation at t=40 days, we get:

So one would never expect the envelope to be completely unbound at this point in the simulation because this would require \alpha_\mathrm{CE} > 1. Therefore, that the envelope is not completely unbound at the end of the simulation is consistent with the alpha prescription.

Then at what value of a_\mathrm{f} should the envelope be completely unbound, for a given value of \alpha_\mathrm{CE}?

So if we say that there are no extra energy sources or sinks and take the most optimistic (but not very realistic) case \alpha_\mathrm{CE}=1, then the final separation is predicted to be about 3R_\odot. Compare this to the softening radius at the end of the simulation of 1.2R_\odot.

In the probably more realistic case of \alpha_\mathrm{CE}=0.1, we would have a final separation of about 0.3R_\odot. This would imply a merger if the secondary is a main sequence star, but not necessarily if the secondary is a white dwarf.

An AGB star is less tightly bound so would be more promising for ejecting the envelope.

Roche lobe overflow tests

​Face-on density (zoomed in) (Run 156 a=109Rsun), equal to theoretical Roche limit separation for this system.
​Face-on density (zoomed in) (Run 157 a=98Rsun)
​Face-on density (zoomed in) (Run 159 a=73.5Rsun)
​Face-on density (zoomed in) (Run 158 a=49Rsun)

The sound-crossing time for the RG is about 8 days or ~35 frames. This is the approximate time the star would take to deform to fill the Roche lobe.

While this deformation is happening, the secondary accretes from the ambient medium. The freefall time onto the point mass m_2 is t_\mathrm{ff}\sim\frac{\pi}{2}\frac{R^{3/2}}{\sqrt{2\mathrm{G}m_2}}. So in 5 days, or about 20 frames, this corresponds to gas being accumulated out to a sphere of radius \sim40\mathrm{R}_\odot from the secondary, corresponding to about 3\times10^{-4} \mathrm{M}_\odot of gas falling onto the secondary from the ambient medium (\rho_\mathrm{amb}=6.7\times10^{-9} grams/cc).

There are numerical problems in the low resolution runs that cause the RG to bulge out opposite to the direction of motion along the orbit.

It is probably worth starting a full resolution run to see what happens. But at which separation?

To-do list

• recombination equations (priority)
• sort out high momentum wind result
• continue running low momentum wind

Low Momentum Wind

Speed = 300 km/s, density 5e-13 g/cm³

Run to 7400 days (1400 days after quiescent phase) Density movie:
​http://www.pas.rochester.edu/~yzou5/PNe_simulation_movies/Old_movies/Low_momentum_Day7400_rho.gif

Next, adjust refinement region and level to continue.

High Momentum Wind

Run to 8000 days (2000 days after quiescent phase). Will update a link to the results as soon as they are done.

Maybe redo the runs for better consistency in frame rate and refinement?

Edge-on through both particles

Density: Edge-on through both particles (full box): ​https://www.pas.rochester.edu/~yishengtu/CEE_gp_meeting_09242018/Zoom_center_p1_edge_on_through_par_rho/D570.gif

Density: Edge-on through both particles (Zoom-in): ​https://www.pas.rochester.edu/~yishengtu/CEE_gp_meeting_09242018/Zoom_center_p1_edge_on_through_par_rho/D150.gif

Normalized energy: Edge-on through both particles (full box): ​https://www.pas.rochester.edu/~yishengtu/CEE_gp_meeting_09242018/Zoom_center_p1_edge_on_through_par_En/Etot_gas_normal/D570.gif

Normalized energy: Edge-on through both particles (Zoom-in): ​https://www.pas.rochester.edu/~yishengtu/CEE_gp_meeting_09242018/Zoom_center_p1_edge_on_through_par_En/Etot_gas_normal/D150.gif

absolute value of energy: Edge-on through both particles (full box): ​https://www.pas.rochester.edu/~yishengtu/CEE_gp_meeting_09242018/Zoom_center_p1_edge_on_through_par_En/Etot_gas_abs/D570.gif

absolute value of energy: Edge-on through both particles (Zoom-in): ​https://www.pas.rochester.edu/~yishengtu/CEE_gp_meeting_09242018/Zoom_center_p1_edge_on_through_par_En/Etot_gas_abs/D150.gif

Parallel to edge-on through particles but center at center of box

Density: center of box parallel to Edge-on through both particles (full box): ​https://www.pas.rochester.edu/~yishengtu/CEE_gp_meeting_09242018/Zoom_center_of_box_edge_on_through_par_rho/D570.gif

Density: center of box parallel to Edge-on through both particles (Zoom-in): ​https://www.pas.rochester.edu/~yishengtu/CEE_gp_meeting_09242018/Zoom_center_of_box_edge_on_through_par_rho/D150.gif

Normalized energy: center of box parallel to Edge-on through both particles (full box): ​https://www.pas.rochester.edu/~yishengtu/CEE_gp_meeting_09242018/Zoom_center_of_box_edge_on_through_par_En/Etot_gas_normal/D570.gif

Normalized energy: center of box parallel to Edge-on through both particles (Zoom-in): ​https://www.pas.rochester.edu/~yishengtu/CEE_gp_meeting_09242018/Zoom_center_of_box_edge_on_through_par_En/Etot_gas_normal/D150.gif

absolute value of energy: center of box parallel to Edge-on through both particles (full box): ​https://www.pas.rochester.edu/~yishengtu/CEE_gp_meeting_09242018/Zoom_center_of_box_edge_on_through_par_En/Etot_gas_abs/D570.gif

absolute value of energy: center of box parallel to Edge-on through both particles (Zoom-in): ​https://www.pas.rochester.edu/~yishengtu/CEE_gp_meeting_09242018/Zoom_center_of_box_edge_on_through_par_En/Etot_gas_abs/D150.gif

Radiation pressure

Communicating with TACC about missing tmi fabric. While loop appears to be working (should be blocking, but it's hard to tell when it immediately errors…), so should be working once the execution problem is straightened out.

BlueHive predicts 2.3 years, before being killed (due to long timestep), so it's not going to be effective to run even with a large reservation.

HD209458b short paper

Postprocessing of HD209458b run in progress. I'll finish some of the analysis in the next few days. Things to look at right now:

• Tau/xi terms (ratios of dimensionless Coriolis, tidal, wind accelerations)
• Regimes (should be energy-limited; probably type II from Matsakos?)
• Velocities and densities along sub-/anti-stellar ray, with markers for sonic and Hill radius (Coriolis?)
• Synthetic observations - high v?
• Compare all of these with most similar (high mass, low flux) planet from parameter space paper

Plan for XSEDE proposal 2018

Summary of current allocation

• A bit less than 130,000 node hours remain.
• It is estimated (still very rough though) that we can complete about 15 runs like the non-accretion run from paper I (run 143).
• We have thus far used about 23% of the allocation so about 77% remains.
• To get this down to <50% by Oct 15, we need to use up about 27% of the full allocation or a little over 1/3 of what remains.
• This translates to about 5-6 runs by Oct 15.

Plan for now until Oct 15

The following runs do not involve major changes to the code so can hopefully be done in this time frame:

1. Test dependence on secondary mass (parameter regime 1).
   1. Run 149: As 143 but with secondary mass smaller by a factor of 2, i.e. m2=0.5 Msun instead of 1.0 Msun (running).
2. Test dependence on secondary mass (parameter regime 2).
   1. Run 151: As 143 but with secondary mass smaller by a factor of 4, i.e. m2=0.25 Msun instead of 1.0 Msun (submitted).
3. Test convergence with respect to softening length and resolution.
   1. Run 152: As 143 restarted from frame 72, which is the frame where the softening radius and smallest resolution element were halved in 143, but now do not halve the softening length nor the resolution element (submitted).
   2. As 152 but halve softening length only.
   3. As 152 but halve smallest resolution element only.
4. Roche lobe overflow.
   1. Put secondary at Roche limit separation.
   2. Must think about refinement.
   3. Increase q to 1 both to reduce Roche limit separation and also increase |\dot{a}/a|.
5. Test dependence on initial spin of the RG.
   1. Repeat run 143 but now initialize RG in solid body rotation at the initial orbital angular velocity.
6. Test convergence with respect to size of refinement region (parameter regime 1).
   1. Repeat run 143 using a more liberal choice for the refinement radius for each interval (c.f. Fig 3 in Paper I)
7. Test convergence with respect to size of refinement region (parameter regime 1).
   1. Based on the results of the above, repeat run once more using either an even more liberal refinement radius or else a refinement radius that is more conservative than run 143.

• Another possibility is to restart run 143 from the end of the simulation at t = 40 days. This is easy to do but the drawback is that material will flow out of the box soon after.

Plan for Oct 15 to Dec 31

We can expect about 10 more runs equivalent to run 143 but I provide 15 possibilities below:

Questions to answer for jet project:

• How does the jet affect the morphology of the envelope?
• How does the jet affect the ejection of the envelope?
• How does the jet evolve, does it get quenched?
• What is the dependence of these questions on when the jet gets turned on?
• What is the dependence on accretion rate?
• Opening angle?

1. Simulate CEE with a jet from the secondary (parameter regime 1) (Restart from run 143)
2. Simulate CEE with a jet from the secondary (parameter regime 2) (Restart from run 143)
3. Simulate CEE with a jet from the secondary (parameter regime 3) (Restart from run 143)
4. Simulate CEE with a jet from the secondary (parameter regime 4) (Restart from run 143)
5. Simulate CEE with a jet from the secondary (parameter regime 5) (Restart from run 143)
6. Simulate RGB/AGB star in improved simulation including (part I of a long run)
   1. larger box (but smaller ambient resolution): the larger the box the longer we can run it before material starts leaving the box,
   2. ambient consisting of hydrostatic atmosphere surrounded by very low density and low pressure gas,
   3. start from Roche limit separation (optional)
7. Simulate RGB/AGB star in improved simulation including (part II of a long run)
8. Simulate RGB/AGB star in improved simulation including (part III of a long run)
9. Simulate RGB/AGB star in improved simulation including (part IV of a long run)
10. Simulate RGB/AGB star in improved simulation including (part V of a long run)
11. Triple system involving a planet or tertiary star orbiting the secondary… (parameter regime 1)
12. Triple system involving a planet or tertiary star orbiting the secondary… (parameter regime 2)
13. Triple system involving a planet or tertiary star orbiting the secondary… (parameter regime 3)
14. Triple system involving a planet or tertiary star orbiting the secondary… (parameter regime 4)
15. Triple system involving a planet or tertiary star orbiting the secondary… (parameter regime 5)

Storage

• The size of run 143 data is about 10 TB.
• To do 15 more runs like 143 would require 150 TB.
• This would cost an additional $15,000 per year on bluehive which is unaffordable.
• We need to look for other options like free space somewhare on blue hive.
• We probably need to keep only ½ of the data, i.e. every second chombo file.

self-gravity

To do list:

• sharp density gradient
• try to make the Poisson solver work properly
• ambient density

ratio between Laplacian of potential and density

wind

• started low momentum wind
• high momentum wind extended to 7500 days

Radiation Pressure

No real progress to report, other than that stampede is beginning to test my patience. But I believe I've solved all of the problems (except the MPI issue - more below), so we should at least start getting some frames.

MPI issues

Located all of our calls to IRECV and ISEND, and all matching WAITs are there. No MPI errors are going ignored. Is there anything else we can do to be convinced that the problem's not on our end? (Here's ​one example of a case where irecv is called by RECV, which could potentially cause the issues we're seeing but not be anything we can affect).

I've found the environment variable that controls the request descriptors (not sure why I couldn't find it before, but I distinctly remember not being able to set it), and upped it to 10¹⁴, so hopefully that will be sufficient to allow us to run for the full two days we're allowed.

HD209458b ionizing flux only

Got some analysis done here. Mass loss rates are about an order of magnitude lower, which is to be expected given the tighter binding of the atmosphere.

Mass loss rate: ~5.3x10⁹ g/s

It is (almost?) just Roche lobe overflow. Here's the Hill sphere (in red) overlaid on the wind, with contours of the dimensionless Coriolis length.

MHD planets

Here's the global steady state.

Global simulations are actually much worse than local. The run stalls out after 6 frames. Running for 2 more days, we get three more frames; here's the last one output:

Stellar wind implementation

I have a Matlab script that's solving for the (radial) velocity, pressure, and density of a pseudo-stellar wind with adiabatic gamma. I have noticed a couple of odd jumps that seem to be artifacts of the numerical solver, however. Also, I can't solve for any points inside the reference radius. This seems like a Matlab (or, specifically, vpasolve) limitation, but I haven't looked into other numerical methods in Matlab or otherwise. I've attached a sample solution file​.

Current list

I'll just keep this as a running tally, to keep us all organized a bit:

To Do

1. Continue radiation pressure runs
2. MHD planet runs (use the larger local runs as good-enough?)
3. Invert subcycling for line transfer
4. Add MHD to line transfer run
5. Test John’s stellar wind for charge exchange
6. Test charge exchange
7. Run charge exchange
8. Development for AMR line transfer/point source line transfer
9. Development for metastable He

Complete

Debug MPI issue on stampede (ish)
Set up global MHD run (w/ stellar wind)

Summary of last week's work

1. Reran first part of simulation with secondary mass halved, after deriving orbital parameters and carefully testing the initial orbit ​p_mult_143_149.png
2. Worked on energy paper (on the section dealing with the energy prescription with \alpha formulation) ​en.pdf.
3. Worked with Thomas (Orsola's student) to set up an AGB star for CE simulations.
4. Worked with Amy to debug the CE-wind with self-gravity simulation.

All of these are ongoing—-I expect to report more progress next Monday.

Working on first draft of energy paper

​http://www.pas.rochester.edu/~lchamandy/CE_papers/Energy/en.pdf

Project/To-Do List

• Radiation pressure
  • Run high/med/low flux
  • Debug MPI issue (using too many descriptors)
• MHD planets
  • Fill parameter space
  • Add fields to parameter space run
• Charge exchange
  • Get initial conditions working (stellar wind without charge exchange populations)
  • Test charge exchange (have an attempt to compare to Christie et al [with fixed boundary code, though] - any other ways?)
  • Run (with suite of stellar wind strengths as in John's paper, perhaps? Should maybe look at what's predicted to affect the efficacy of charge exchange)

MHD planets

Refer to ​Jonathan's page for simulation parameters.

Using box of effective resolution 512³ for size (1 CU)³, with sphere of radius 0.25 CU forced to be resolved to level 1 (so it's equivalent to Jonathan's previous runs). Currently running corotating sims (do we want non-rotating as well?). Beta of 0.1 corresponds to surface magnetic field of 0.23 G.

Corotating

Have hydro steady state:

Also have a late state of the beta=0.1 run. May actually be a couple of frames too long, with some of the field lines being distorted down-orbit by stellar wind/planetary wind interactions? Running time is increasing, but only 4 days remaining when I restart (worth trying to get a few more frames?).

Beta = 10 run hit a wall at frame 224 - 1.9 months remaining.

Radiation Pressure

Run with no radiation pressure appears to have reached steady state. Actually not too different from frame 180 where radiation pressure is introduced in new runs.

Have a couple of frames from the high flux and medium flux cases. Just had stampede queue access restored, so will be restarting them.

High

Med

Debugging

Fixed makefile issue for debugging on stampede. Now that I have queue access again, should be able to use DDT to see if there's a leaky thread somewhere.

Movies

Outflow

370 days after quiescent phase ​http://www.pas.rochester.edu/~yzou5/PNe_simulation_movies/Old_movies/High_momentum_Day6370_rho.gif

Quiescent phase

6000 days ​http://www.pas.rochester.edu/~yzou5/PNe_simulation_movies/Old_movies/High_momentum_D6000_Qphase.gif

Refinement region

for frames 380-392 (t = 6370 days). it takes ~6 hours to generate 1 frame. dimensions are in solar radius.

level 1-4

(x, y, z) = +-(5000, 5000, 10000)

level 5-7

(x, y, z) = +-(1000, 1000, 6000)

Instead of making the wedge tip same as the ambient as in blog:bliu08012018 , tried with Bruce's idea of "t=0 simply replicate the flow at the edge of the launch surface into the wedge". Seems working much better than the scheme in blog:bliu08012018

t=0
t=400y
Movie up to 400y	​movie