For my masters level thesis, I plan complete and analyze MHD simulations of an accretion disk around the center of a RG and AGB star to investigate the stability and magnetic field generation . This is in an attempt to computationally verify the suggestion that High Field Magnetic White Dwarfs are the result of a Common Envelope, Nordhaus et al. 2011:

Link:​https://arxiv.org/abs/1010.1529

To do this I have looked to Ohlmann’s thesis and previous blog posts made by Luke:

1.) Run The Mesa Stellar Evolution to Generate density, pressure, and temperature profiles at different points in time.

2.) The idea in Nordhaus et al. 2011 is that magnetic field is built up by the engulfment and the accretion of a companion. This is most likely when the radius is at a maximum, So I identified these maximum. (Radius (log[R/Rsun]) vs Time (years) for a 2 solar mass star)

3.) Modify those profiles with the modified-lane-emdan equations in accordance with Ohlmanns thesis (Link:​https://archiv.ub.uni-heidelberg.de/volltextserver/21513/1/thesis-sohlmann.pdf) (Lukes Code)

Link:​http://www.pas.rochester.edu/~lchamandy/Presentations/profile_blog1.pdf

4.) Set the central mass by hand such that the mass contained within the smoothing radius is the same as before modification. Place onto grid with at least 10 cells across smoothing radius

5.) If the central density is stable then I can move on to adding an accretion disk to the central region of the simulation.

Disk parameters: R ~ 10e10 cm

H ~ 10e9 cm M ~10 Jupiter Masses

Radiation Pressure

• Static radiation pressure test looks good. Setting up some dynamic tests. (Normalize min vel to sound speed?)

• Radiation/recombination-limited simulation (same as previous, but with 2x10¹⁷ max flux) needs to have ramping adjusted - larger flux just obliterates the atmosphere. I'm thinking start from 10^-4 or so lower, and adjust the ramping speed and half time accordingly.

• Need to work out photon-limited parameters — mass of planet as given is too low for the radius and temperature and results in an infinitely extended atmosphere.

New Work

1. Addressing a difficult issue with making VisIt plots.
2. More analysis if high res run 132
   • Movies centered on P2 that include velocity vectors (in inertial frame).
   • Graph of accretion rate vs time.

Summary of New Results

1. Extracting information from one of the particles to use in the variable expressions in VisIt is proving to be extremely non-trivial. Baowei is working on this, but it probably merits a workaround from within astrobear (defining each particle to be a separate subset of the mesh 'particles' might do the trick). Otherwise, as it stands now, one cannot plot within VisIt, for example, the gravitational potential due to a particle, or the gas velocity vectors in the reference frame of the particle.
2. Results:
   • Edge-on movies show a transient global vertical flow through the center of the torus that changes direction with time.
   • The accretion rate from the Krumholz prescription is of order 1 Msun/yr, which seems too large.

Run 132 that uses twice as high resolution and 10x lower ambient pressure as previous runs:

Binary run 132 with double max resolution, lower resolution in ambient medium, 10x smaller ambient pressure than run 116
Relaxation run: 129
First frame: 75 (5 RG freefall times, when velocity damping ended)
Last frame:
Total simulation time:
Machine and partition: Stampde 1 normal/Bluehive 2.5 standard/Bluestreak standard (completed up to frame 375, or 70 sim-days)
Number of cores: 1024
Total wall time: around 14 days (starting from frame 75 of relaxation run)
Hydro BCs: extrapolated
Poisson BCs: multipole expansion
Box size: L=4e13 cm (575 Rsun)
Base resolution: 2.25 Rsun (256³ cells)
Highest resolution: 0.14 Rsun (4096³ cells, 4 levels AMR)
AMR implementation: set by hand to have max level around point particles
Max resolution zone: within 5e12 cm (71.87 Rsun) of primary center and within a cylinder of radius 20 Rsun and height 20 Rsun around secondary center. After t~11d (~frame 123) refinement radius around primary was halved to 2.5e12cm (36Rsun). After t~31d (frame 210) halved again to 1.25e12cm (18Rsun).
Buffer zones: 2 cells
Softening length: 2.4 Rsun
Ambient density: 6.7e-9 g/cc
Ambient pressure: 10⁵ dyne/cm²
DefaultAccretionRoutine=2 (Krumholz)

New Movies of run 132
1) Edge-on , centered at P2, with velocity vectors (relative to inertial frame). Left: slice through both particles. Right: slice through P2 as viewed from direction of P1.
​edge-on with velocity vectors (relative to inertial frame)
2) Face-on , centered at P2, with velocity vectors (relative to inertial frame).
​face-on with velocity vectors (relative to inertial frame)

Comments:

• Movie (1) shows that there is a significant vertical flow of gas of order up to a few times 10⁷ km/s.
  • The velocity vectors are in the inertial frame whereas the movie is in the frame of P2. However, the orbital velocity of P2 is at all times less than 10⁷ km/s so the largest vectors shown would be almost unchanged if they were plotted with respect to the velocity of P2. Moreover P2 has close to zero vertical velocity.
  • The vertical flow is present only when the torus is present.
  • Even when the torus is present, the vertical flow is intermittent.
  • When the vertical flow is present it is located in the 'hole' of the torus where the density is low.
  • When the vertical flow is present the flow is either upward below and above the orbital plane, or downward below and above the orbital plane, unlike a jet, which would be upward above the midplane and downward below.
  • Roughly speaking, it switches from downward (frames 36-40) to non-existent to upward (frames 75-99) to non-existant to downward (frames 119-152) to non-existent to downward (frames 164-168).
• Movie (2) shows rotation around the secondary.

Accretion rate

Next steps

• Compute accretion rate using alternative prescription using mass within a sphere centered at the secondary.
• Compute drag force.
• New simulation with a jet.
  • Use outflow object in astrobear.
  • Can take mdot from RT08/12.
  • Velocity of order 100 km/s (or escape veloc at ~1 solar radius say).
  • Half opening angle ~10 degrees.
• Plot CM position of system as a function of time.

First Set of Parameters

See here for an enumeration of the possible parameter space with orbital distance, flux, and planet mass variable. What follows is the top left of the first table.

See here for other simulation parameters.

Top view

Side view

Radiation Pressure

Testing with slab of density 1 (in CU), all neutral, embedded in ambient of density 0.01, all ionized. Setting J = \frac{n_{H} g}{\sigma_H p_\gamma}, the whole thing should be in a stationary state (the domain is thin enough that \tau over the grid is close to 0). Some problems with initial test conditions (I don't think I've sufficiently isolated the effects of radiation pressure).

I've restarted test with no hydrodynamics and constant temperature rather than pressure (\sigma_H depends on T), and the entire domain neutral. Seems like having hydrodynamics turned off removes gravity, though. From comparing diagnostic variables, seems like scaling is off on the radiation pressure term.

Scaling:

…and fixed.

New Work

1. Notes about extracting diagnostics on energy from the simulation.
2. Analysis of runs:
   • 133 (as low res fiducial run 116 and more aggressive refinement run 125 but with even more aggressive refinement algorithm that de-resolves more of the `ambient' envelope).
   • 135 (as 133 but with accretion turned off).
   • 136 (as 133 but without any relaxation run, so no initial damping of velocities).
   • 132 (similar to 125/133 but with maxlevel increased by 1 so twice higher resolution).
   • 120 (similar to 116 but ambient density reduced by a factor of about 67).

Summary of New Results

1. Procedure to extract energy diagnositics seems quite clear.
2. Key results from runs are as follows:
   • 133: Results are generally consistent with run 116 but the radius of the maxlevel refinement region around the primary point particle should be kept at a minimum of 2 times the particle separation.
   • 135: Results are very similar to 133 with `accretion torus' morphology preserved in spite of accretion being turned off.
   • 136: Results are very similar to 133. This means damping is probably unnecessary and can be omitted in future runs!
   • 132: Results of this first higher resolution run are closer to the result of O+16a. Qualitatively similar to run 133.
   • 120: Results are similar to those of run 116, but expansion and outer shock morphology of expanding envelope are slightly affected.

Notes on energy

I compiled the following notes. Thanks to Bo for a discussion on this:
Please see ​en.pdf

Notes on accretion and drag force (from last blog post)

Please see ​df.pdf

Analysis of runs 133, 135, 136, 132, 120

Old run 116 (for comparison)

Relaxation run: 096
First frame: 75 (5 RG freefall times, when velocity damping ended)
Last frame: 168
Total simulation time: 3.36e6 s or ~4.2 RG sound-crossing times or 21.5 days (93 frames)
Machine and partition: Stampede 1 normal
Number of cores: 1024
Total wall time: 96 hours
Hydro BCs: extrapolated
Poisson BCs: multipole expansion
Box size: L=4e13 cm
Base resolution: 9.0 Rsun (64³ cells)
Highest resolution: 0.29 Rsun (2048³ cells, 5 levels AMR)
AMR implementation: set automatically by AstroBear
Softening length: 2.4 Rsun
Ambient density: 6.7e-9 g/cc
Ambient pressure: 10⁶ dyne/cm²
DefaultAccretionRoutine=2 (Krumholz)

Old run 125 (for comparison)

Binary run 125 with longer simulation time and lower resolution in ambient medium
Relaxation run: 096
First frame: 75 (5 RG freefall times, when velocity damping ended)
Last frame: 335
Total simulation time: 5.2e6 s or 60 days ~6.5 RG sound-crossing times or 60.2 days (260 frames)
Machine and partition: Bluehive standard
Number of cores: 120
Total wall time: about 8.5 days
Hydro BCs: extrapolated
Poisson BCs: multipole expansion
Box size: L=4e13 cm (575 Rsun)
Base resolution: 9.0 Rsun (64³ cells)
Highest resolution: 0.29 Rsun (2048³ cells, 5 levels AMR)
AMR implementation: set by hand to have max level around point particles
Max resolution zone: within 50 Rsun of primary and within a cylinder of radius 50 Rsun and height 50 Rsun around secondary
Buffer zones: 0 cells (no buffer zones)
Softening length: 2.4 Rsun
Ambient density: 6.7e-9 g/cc
Ambient pressure: 10⁶ dyne/cm²
DefaultAccretionRoutine=2 (Krumholz)

Run 133 to test making maxlevel refinement zone decrease with time:

Binary run 133 similar to run 125 but now the refinement zone changes with time
Relaxation run: 096
First frame: 75 (5 RG freefall times, when velocity damping ended)
Last frame:
Total simulation time: 84 days (up to frame 439)
Machine and partition: Bluestreak standard
Number of cores: 8192 (2 cores per task to increase memory)
Total wall time: 6 days
Hydro BCs: extrapolated
Poisson BCs: multipole expansion
Box size: L=4e13 cm (575 Rsun)
Base resolution: 9.0 Rsun (64³ cells)
Highest resolution: 0.29 Rsun (2048³ cells, 5 levels AMR)
AMR implementation: set by hand to have max level around point particles
Max resolution zone: within min(5e12cm, 1.5*particle_separation)
Buffer zones: 2 cells
Softening length: 2.4 Rsun
Ambient density: 6.7e-9 g/cc
Ambient pressure: 10⁶ dyne/cm²
DefaultAccretionRoutine=2 (Krumholz)

Comparison between runs 125(top left), 133(top right), O+16a(bottom) (FROM LAST BLOG POST)

Comparison between run 125(left) and 133(right)
​face-on zoom-in
​face-on slice through secondary, extra zoom-in
​viewed from P1 with P2 at center
​slice through P1 (left side) and P2 (center)

Comments:
FROM LAST BLOG POST (BASED ON ORBITS ONLY):

• The runs are very similar but there is a slightly smaller final separation for run 125(less aggressive refinement)
• This suggests that our refinement criteria for run 133 was probably slightly too aggressive.
• However, the otherwise close agreement tells us that decreasing the refinement zone with time to be a sphere centred on p1 with radius ~2 times particle separation is reasonable.

NEW COMMENTS (BASED ON MOVIES):

• The runs are very similar but having max refinement in a sphere centered around the primary with radius 1.5 times the separation is a bit too aggressive. Could probably get away with about 2 times the separation.

Run 135 to test case where no accretion onto secondary is permitted:

Binary run 135 similar to run 133 but now DefaultAccretionRoutine=0 instead of 2. Also suppress generation of new sink particles.
Relaxation run: 096
First frame: 75 (5 RG freefall times, when velocity damping ended)
Last frame: 900
Total simulation time: 191 days
Machine and partition: bluehive2.5 standard
Number of cores: 120
Total wall time: 2.8 days
Hydro BCs: extrapolated
Poisson BCs: multipole expansion
Box size: L=4e13 cm (575 Rsun)
Base resolution: 9.0 Rsun (64³ cells)
Highest resolution: 0.29 Rsun (2048³ cells, 5 levels AMR)
AMR implementation: set by hand to have max level around point particles
Max resolution zone: within min(5e12cm, 1.5*particle_separation)
Buffer zones: 2 cells
Softening length: 2.4 Rsun
Ambient density: 6.7e-9 g/cc
Ambient pressure: 10⁶ dyne/cm²
Default Accretion Routine=0 (No accretion)

Comparison between runs 133(top left), 135(top right), O+16a(bottom)

Comparison between run 133(left) and 135(right)
​face-on zoom-in
​face-on slice through secondary, extra zoom-in
​viewed from P1 with P2 at center
​slice through P1 (left side) and P2 (center)

Comments:
FROM LAST BLOG POST (BASED ON ORBITS ONLY):

• The results of 135 and 133 are similar, except:
  • Run 135 (no accretion) has slightly smaller separation at a given time.
  • Run 135 (no accretion) has slightly larger frequency of oscillations.
• Thus accretion cannot fully explain the discrepancy with the O+16a results, but removing accretion does make results slighly closer to those of O+16a, who did not have accretion.
• We must think more about accretion, but it is not making a huge difference at present.

NEW COMMENTS:

• Results are very similar. Most noticeable difference is that density is higher around the point particles for run 135, which does not have accretion. This is most clearly seen in the first, most zoomed out movie. This would be expected since gas is not removed by the particle.
• Particle creation was turned off for run 135 and this may explain some differences.
• Bottom line is that the gas flow around the secondary is very similar regardless of whether accretion is turned on. Even without accretion turned on, there is a torus structure with bipolar outflows.
• Note that for both runs the flow around the secondary undergoes a transition after about 75 days (frame 375) and the torus morphology is no longer present at this time. This corresponds to a particle separation of about 9 Rsun.

Run 136 to test how much relaxation makes a difference:

Binary run 136 similar to run 135 but now start from frame 0 of relaxation run 096 instead of frame 75 (so no damping).
Relaxation run: 096 (but no relaxation! Just start with initial profile
First frame: 0
Last frame: 18, running
Total simulation time: 4 days
Machine and partition: bluehive2.5 standard
Number of cores: 120
Total wall time: 21 hours, running
Hydro BCs: extrapolated
Poisson BCs: multipole expansion
Box size: L=4e13 cm (575 Rsun)
Base resolution: 9.0 Rsun (64³ cells)
Highest resolution: 0.29 Rsun (2048³ cells, 5 levels AMR)
AMR implementation: set by hand to have max level around point particles
Max resolution zone: within min(5e12cm, 1.5*particle_separation)
Buffer zones: 2 cells
Softening length: 2.4 Rsun
Ambient density: 6.7e-9 g/cc
Ambient pressure: 10⁶ dyne/cm²
Default Accretion Routine=0 (No accretion)

Comparison between runs 133(left) and 136(right)

Comparison between run 133(left) and 136(right)
​face-on zoom-in
​face-on slice through secondary, extra zoom-in
​viewed from P1 with P2 at center
​slice through P1 (left side) and P2 (center)

Comments:

• The results are very similar. This suggests that relaxation (velocity damping before the run) is NOT necessary. (It also confirms that most of the differences between runs 133 and 135 or 125 and 133 above are not due to the particle creation on run 133.)

Run 132 that uses twice as high resolution and 10x lower ambient pressure as previous runs:

Binary run 132 with double max resolution, lower resolution in ambient medium, 10x smaller ambient pressure than run 116
Relaxation run: 129
First frame: 75 (5 RG freefall times, when velocity damping ended)
Last frame:
Total simulation time:
Machine and partition: Stampde 1 normal/Bluehive 2.5 standard/Bluestreak standard (completed up to frame 375, or 70 sim-days)
Number of cores: 1024
Total wall time: around 14 days (starting from frame 75 of relaxation run)
Hydro BCs: extrapolated
Poisson BCs: multipole expansion
Box size: L=4e13 cm (575 Rsun)
Base resolution: 2.25 Rsun (256³ cells)
Highest resolution: 0.14 Rsun (4096³ cells, 4 levels AMR)
AMR implementation: set by hand to have max level around point particles
Max resolution zone: within 5e12 cm (71.87 Rsun) of primary center and within a cylinder of radius 20 Rsun and height 20 Rsun around secondary center. After t~11d (~frame 123) refinement radius around primary was halved to 2.5e12cm (36Rsun). After t~31d (frame 210) halved again to 1.25e12cm (18Rsun).
Buffer zones: 2 cells
Softening length: 2.4 Rsun
Ambient density: 6.7e-9 g/cc
Ambient pressure: 10⁵ dyne/cm²
DefaultAccretionRoutine=2 (Krumholz)

Movies of run 132
​face-on full box
​face-on zoom-in
​edge-on full box
​edge-on zoom-in
​face-on slice through secondary, extra zoom-in
​viewed from P1 with P2 at center
​slice through P1 (left side) and P2 (center)

Comparison of run 132 with Ohlmann+16a run:
​Table

Comparison between run 133(left) and 132(right)
​face-on zoom-in
​face-on slice through secondary, extra zoom-in
​viewed from P1 with P2 at center
​slice through P1 (left side) and P2 (center)

Comparison between runs 133(top left), 132(top right), O+16a(bottom)

Comments:
FROM LAST BLOG POST (BASED ON ORBITS ONLY):

• The 1st minimum, 2nd maximum, and 2nd minimum are located roughly at separations, with softening length denoted as h:
  • Run 088 (res=0.29Rsun, h=4.8Rsun, cells/h=17): 14 , 21 , 12 Rsun
  • Run 125 (res=0.29Rsun, h=2.4Rsun, cells/h= 8): 13.5, 18 , 14.5 Rsun
  • Run 132 (res=0.14Rsun, h=2.4Rsun, cells/h=17): 14 , 18.5, 11 Rsun
  • O+16a (see above table) : 10 , 23.5, 11 Rsun
• These are found roughly at times:
  • Run 088 (res=0.29Rsun, h=4.8Rsun, cells/h=17): 12.5, 15.5, 18 days
  • Run 125 (res=0.29Rsun, h=2.4Rsun, cells/h= 8): 12.5, 14.5, 17 days
  • Run 132 (res=0.14Rsun, h=2.4Rsun, cells/h=17): 12.5, 14.5, 17 days
  • O+16a (see above table) : 13 , 16 , 19 days
• From the difference between 125 and 132, we conclude that h is not adequately resolved in run 125 (8.5 cells/h vs 17 cells/h)
• O+16a advocates >10 cells/h, so this is consistent with their result.
• In our case, we still cannot tell if 17 cells/h is sufficient.
• We could run 088 with double the resolution to see if it converges.
• The softening length of the point particles decreases with time in the O+16a simulation.
  • It is not clear whether this procedure is justified.
  • The main reason to improve the resolution with time may be to resolve the decreasing softening length.
  • Therefore, if we keep the softening length constant, improving the resolution with time may not be necessary or productive.
• O+16a initially resolves the softening length by 20 cells, whereas we resolve it by 17 cells. They advocate >10 cells.
• At 120 days, O+16a resolves the softening radius by 35 cells.

NEW COMMENTS:

• Runs 133 and 132 are very similar for the first 15 days, or when the separation does not dip below 14 Rsun. This suggests that run 133 has inadequate resolution (or resolution per softening length) when the separation is <14 Rsun. Recall that both runs have softening length 2.4 Rsun but that run 133 has best resolution of 0.29 Rsun while run 132 has best resolution of 0.14 Rsun.
• The amplitude of the oscillations in the separation is about the same for run 132 and 133, both lower than that of O+16a.
• The final asymptotic separation is about 6 Rsun for run 132, compared to about 9 Rsun for run 133 and 4 Rsun for O+16a.
• The transition of the flow around the secondary happens earlier in run 132, at about 55 days (frame 320) compared to about 75 days (frame 375) in run 133. For run 132, this corresponds to a mean particle separation of about 7 Rsun, and for run 133, to a mean particle separation of about 9 Rsun. This supports the idea that the transition to a flow without the torus morphology is probably related to the separation. Further, it also supports the idea that such a transition is an artifact of having a too-low resolution per softening length, because the separation at onset is lower when the resolution is higher.
• There is a tendency for the particles to migrate in the direction of larger x and y. This is visible to a lesser extent in O+16a. This should be explained.

Run 120 to test importance of density of ambient medium

Relaxation run: 096
First frame: 75 (5 RG freefall times, when velocity damping ended)
Last frame: 119
Total simulation time: 3.36e6 s or ~4.2 RG sound-crossing times or 10.2 days (48 frames)
Machine and partition: Stampede 1 normal
Number of cores: 1024
Total wall time: 6 days
Hydro BCs: extrapolated
Poisson BCs: multipole expansion
Box size: L=4e13 cm
Base resolution: 9.0 Rsun (64³ cells)
Highest resolution: 0.29 Rsun (2048³ cells, 5 levels AMR)
AMR implementation: set automatically by AstroBear
Softening length: 2.4 Rsun
Ambient density: 1e-10 g/cc
Ambient pressure: 10⁶ dyne/cm²
DefaultAccretionRoutine=2 (Krumholz)

Comparison between runs 116(left) and 120(right)

Comparison between run 116(left) and 120(right)
​face-on
​face-on different color bar limits
​face-on slice through secondary, extra zoom-in

Comments:

• Results of the two runs are very similar.
• Differences are apparent at the interface of the ambient medium and envelope. Specifically, whisps of low density (Kelvin-Helmholtz unstable?) is visible outside the shock in low ambient density run 120. Also, by the end of the run, the envelope has expanded to a slightly larger size in the low ambient density case (the larger upper lobe is about 10% wider). The outer shock structure also shows minor differences.
• In the region around the point particles, differences are very minor. If our focus is on this region, we should not have a problem with the higher ambient density.

Discussion & conclusions

• The energy contributions can be readily analyzed from the chombo files, now that we have written down the important equations.
• Damping is not necessary (and introduces artifacts through grid effects). Can avoid it in future.
• Turning accretion off does not have a big effect on the flow. Since we do not understand how to model the accretion very well, perhaps we should keep it turned off for now. This would avoid the unphysical removal of pressure during the accretion process, as pointed out by Eric. In addition, the analysis of the drag force and energy budget would be simplified, and we could use the chombo files without having to output additional information. But we need to get an accretion rate to use for the jets, so ultimately we need to decide on an accretion model.
• We are able to be reasonably aggressive with refinement, but max level refinement in a sphere around the RG core with radius of 2 times the particle separation seems to be the limit, and trying to make this region even smaller introduces too much error.
• Ambient density does not seem to matter very much, but not sure yet with regard to binding and ultimate escape of envelope.
• High resolution run is consistent, qualitatively, with lower resolution run, but shows quantitative results closer to those of O+16a. Namely smaller final separation and orbital period than our lower resolution runs, but still not as small as O+16a, and amplitude of oscillations still lower than O+16a.
• There is a transition from a torus/thick disk structure to flow which does not show this structure. This transition seems to happen when the inter-particle separation reduces past a certain threshold. Futhermore, the threshold seems to be smaller for a better resolution, suggesting that this transition may be a numerical artifact due to a lack of resolution (or resolution per softening length).
• We should keep track of the center of mass of the system (gas+particles) which should not change. (Some gas leaves the box during the simulation and it might be worth increasing the box size or at least testing the effect of increasing it.)

Next steps

• We should now start to perform the analysis of the drag force and energy budget using results of run 132. We don't have accreted momentum outputted from that run, but we can still make do with the chombo files for most of the analysis.
• Simultaneously, we should be doing runs on stampede 2 and comet with our remaining resources, as well as on bluestreak and bluehive:
  • Rerun from middle of 132 but with softening length now able to decrease with time. Need to increase AMR max level with time accordingly.
  • New run like 132 but with initial RG spin equal to 95% of corotation to see what effect this has.
  • New run like run 088 (low res, twice the softening length of more recent runs) but with twice as high resolution (i.e. same resolution as run 132) to test the convergence of results as the number of resolution cells per softening length increases (currently we know that 8 is insufficient and 17 may or may not be sufficient…here we would bump it up to 34).

Appendix

Movies of run 133
​face-on zoom-in
​face-on slice through secondary, extra zoom-in
​viewed from P1 with P2 at center
​slice through P1 (left side) and P2 (center)

Movies of run 135
​face-on zoom-in
​face-on slice through secondary, extra zoom-in
​viewed from P1 with P2 at center
​slice through P1 (left side) and P2 (center)

Movies of run 136
​face-on zoom-in
​face-on slice through secondary, extra zoom-in
​viewed from P1 with P2 at center
​slice through P1 (left side) and P2 (center)

Movies of run 120
​face-on
​face-on different color bar limits

I put a 1 solar mass object at the origin of the simulation box. There is no outflow from 0 to 500 day. From 500 to 1000 day, there is outflow for the boundary of the object. I list the physical quantities below:

Radius of the object that has outflow: r=46.875R_{\odot}

Mass of the object: 1M_{\odot}

Density of the outflow: \rho=10^{-11}g\cdot cm^{-3}

Temperature of the object: T_{1}=30000K

The outflow is asymmetric, in the polar direction (theta=0 and theta=pi), the outflow velocity is the maximum; in the equatorial plane, the outflow velocity is 0. The outflow velocity is a function of polar angle, there are two kinds of functions, linear and Gaussian - to vary the open angle. The Gaussian is narrower.

Linear: v_{r}(\theta)=v_{max}\times\frac{\|\pi-2\theta\|}{\pi}

Gaussian: v_{r}(\theta)=v_{max}\times\exp(-10(\frac{\|\pi-\|2\theta-\pi\|\|}{\pi})^{2})

Spherical: v_{r}(\theta)=v_{max}

Below shows the linear and Gaussian velocity v.s. angle.

The outflow temperature: 30000K

Outflow duration: 500 day - 1000 day

Mass loss: \dot{M}=\int_0^{\pi}2\pi\rho r^{2}\sin\theta v_{r}(\theta)d\theta. However, much of the outflow will be impeded by the surrounding gas, the actual mass loss will be lower than the calculated ones.

ID	v_{r}(0)(km\cdot s^{-1})	\rho (g\cdot cm^{-3})	\dot{M}(M_{\odot}\cdot yr^{-1})	simulation	approximate expansion velocity
1	300 (linear)	10^{-11}	2.31\times10^{-4}		120(km\cdot s^{-1})
				movies:​density and temperature
2	400 (linear)	10^{-11}	3.07\times10^{-4}
3	400 (gaussian)	10^{-12}	1.00\times10^{-5}
4	400 (gaussian)	10^{-11}	1.00\times10^{-4}		310(km\cdot s^{-1})
				movies:​density and temperature
5	300 (spherical)	10^{-11}	6.35\times10^{-4}		240(km\cdot s^{-1})
				movies:​density and temperature

Simulations that has indefinite outflow duration. (it is always on)

Some questions: where exactly does the jet emerge? What are the open angle and outflow velocity. What are the outflow density and temperature?

New Work

1. Notes about extracting diagnostics on accretion and drag force from the sim.
2. Work on Xsede proposal.
3. Improvement of high res run 132 to suppress unwanted extra refinement.

Summary of New Results

1. Procedure to extract accretion and drag force diagnostics is now quite clear.
2. Xsede proposal is coming along okay.
3. Mysterious refinement to the max level was happening outside the volume for which the error flags are set to 1, slowing down the code. I couldn't prevent it by making changes to global.data so I decided to manually suppress it in problem.f90 (fortunately buffer zones turn out to be preserved).

Notes on accretion and drag force

Please see ​df.pdf

Discussion

• There are a few issues that need to be discussed going forward:
  1. Type of accretion (​Federrath+10 seems more appropriate to me than ​Krumholz+04)
  2. I've turned off particle creation—is this fine?
  3. Should we reduce the softening length and enhance the resolution with time as the particles get closer, as done by ​Ohlmann+16a? This seems reasonable except
     • do the gains justify the increase in resources?
     • is it really meaningful/beneficial/physical to change a "fuzzy" point particle into a less fuzzy point particle while the sim is in progress?
     • is it at least worth doing some tests to see what effect this would have?
  4. How important is it to use a larger box?
     • Maybe it would be enough to test this by doing one large-box run for comparison?

Next steps

• Do for energy and angular momentum what has been done for accretion rate and drag force (see above notes).
• Analysis of runs 136 (no relaxation), 120 (low density ambient medium) and 132 (high-res run, in progress)
• New run with initial spin of RG at 95% Keplerian to test what difference this makes.
• New run like run 088 (low res, twice the softening length of more recent runs) but with twice as high resolution (i.e. same resolution as current run 132) to test the convergence of results as the number of resolution cells per softening length increases (currently we know that 8 is insufficient and 17 may or may not be sufficient…here we would bump it up to 34).

Timing

See previous post for most simulation parameters.

Timing estimates that we discussed at Monday's meeting are inaccurate - I think there may have been a soft crash.

Wind

Line transfer evolving steadily. Still looks good at about two crossing times further than last post.

Top

Side

Fixed boundary temperature profile is backwards, but I don't think I'd expect that to affect the evolution other than reflecting it over the y=x line.

Also still working on box size/flux ramp test, but now I'm getting dt=0 on BlueStreak. Will be looking into this.

1D

After quick discussion with Jonathan, added spherical geometric terms (dilution and cooling) to code. Queued on BlueHive.

Line transfer

Stellar gravity and Coriolis test went well.

Parameters are the same for linked line transfer run, except:

Side view

​mesh

Top view

​mesh

Contour key:

Next

Need to test bigger box on BlueStreak (and perhaps get OpenMP working?). Since line transfer step happens at maximum resolution, the resolution of the box plays only a small role until a large portion becomes resolved. Pace is ~50 frames in 20 hours.

Fixed boundary

The equilibrium state is too close to the corrected Parker solution for any evolution to be seen. At 10^-11 CT, the wind is the same as at 1 CT. Nothing more to be done here.

1D comparison

Physical Parameters

Simulation Parameters

lScale = R_P

Got a run that seems to have reached steady state. Stays well below supersonic, though. Still need to get lineouts to compare with Ruth's '09 paper, but the problem seems to be that as the wind launches, the \tau=1 surface moves out so that the wind is launching from what seems to be a very different planet.

Next

Comparison of the ionization diagnostic variables should show why the wind isn't staying sufficiently ionized for an insignificant optical depth. Results should lead toward either futzing with the code or altering the simulation parameters.

New Work

1. Analysis of runs 132 (high res), 133 (like 125 of last post but max refinement volume that decreases with time), 135 (like 133 but no accretion).
2. New run 136 to test importance of RG damping (in progress).

Summary of New Results

1. Each of the runs analyzed gives useful information:
   • Run 132 shows that doubling the resolution matters, e.g. increases frequency of oscillations, decreases final separation.
   • Run 133 shows that it is quite reasonable to use max refinement within a radius min(5e12cm, 1.5*particle_separation) about the primary.
   • Run 135 shows that accretion has a small effect on particle separation (decreases) and orbital frequency (increases).
2. New run 136: first 4 days so far consistent with runs that start from "relaxed" RG.

Analysis of runs 133, 135, 132

Old run 088 (for comparison)

Binary run 088 with longer simulation time and lower resolution in ambient medium
Relaxation run: 062
First frame: 75 (5 RG freefall times, when velocity damping ended)
Last frame: 161
Total simulation time: 20 days
Machine and partition: Comet normal
Number of cores: 1728 (2 cores per task to increase memory)
Total wall time: 4 days
Hydro BCs: extrapolated
Poisson BCs: multipole expansion
Box size: L=4e13 cm (575 Rsun)
Base resolution: 9.0 Rsun (64³ cells)
Highest resolution: 0.29 Rsun (2048³ cells, 5 levels AMR)
AMR implementation: set internally by astrobear
Max resolution zone: n/a
Buffer zones: n/a
Softening length: 4.8 Rsun
Ambient density: 6.7e-9 g/cc
Ambient pressure: 10⁶ dyne/cm²

Old run 125 (for comparison)

Binary run 125 with longer simulation time and lower resolution in ambient medium
Relaxation run: 096
First frame: 75 (5 RG freefall times, when velocity damping ended)
Last frame: 335
Total simulation time: 5.2e6 s or 60 days ~6.5 RG sound-crossing times (260 frames)
Machine and partition: Bluehive standard
Number of cores: 120
Total wall time: about 8.5 days
Hydro BCs: extrapolated
Poisson BCs: multipole expansion
Box size: L=4e13 cm (575 Rsun)
Base resolution: 9.0 Rsun (64³ cells)
Highest resolution: 0.29 Rsun (2048³ cells, 5 levels AMR)
AMR implementation: set by hand to have max level around point particles
Max resolution zone: within 50 Rsun of primary and within a cylinder of radius 50 Rsun and height 50 Rsun around secondary
Buffer zones: 0 cells (no buffer zones)
Softening length: 2.4 Rsun
Ambient density: 6.7e-9 g/cc
Ambient pressure: 10⁶ dyne/cm²

Run 133 to test making refinement zone decrease with time:

Binary run 133 similar to run 125 but now the refinement zone changes with time
Relaxation run: 096
First frame: 75 (5 RG freefall times, when velocity damping ended)
Last frame:
Total simulation time: 84 days (up to frame 439)
Machine and partition: Bluestreak standard
Number of cores: 8192 (2 cores per task to increase memory)
Total wall time: 6 days
Hydro BCs: extrapolated
Poisson BCs: multipole expansion
Box size: L=4e13 cm (575 Rsun)
Base resolution: 9.0 Rsun (64³ cells)
Highest resolution: 0.29 Rsun (2048³ cells, 5 levels AMR)
AMR implementation: set by hand to have max level around point particles
Max resolution zone: within min(5e12cm, 1.5*particle_separation)
Buffer zones: 2 cells
Softening length: 2.4 Rsun
Ambient density: 6.7e-9 g/cc
Ambient pressure: 10⁶ dyne/cm²
DefaultAccretionRoutine=2 (Krumholz)

Comparison between runs 125(top left), 133(top right), O+16a(bottom)

Comments:

• The runs are very similar but there is a slightly smaller final separation for run 125(less aggressive refinement)
• This suggests that our refinement criteria for run 133 was probably slightly too aggressive.
• However, the otherwise close agreement tells us that decreasing the refinement zone with time to be a sphere centred on p1 with radius ~2 times particle separation is reasonable.

Run 135 to test case where no accretion onto secondary is permitted:

Binary run 135 similar to run 133 but now DefaultAccretionRoutine=0 instead of 2. Also suppress generation of new sink particles.
Relaxation run: 096
First frame: 75 (5 RG freefall times, when velocity damping ended)
Last frame: 900
Total simulation time: 191 days
Machine and partition: bluehive2.5 standard
Number of cores: 120
Total wall time: 2.8 days
Hydro BCs: extrapolated
Poisson BCs: multipole expansion
Box size: L=4e13 cm (575 Rsun)
Base resolution: 9.0 Rsun (64³ cells)
Highest resolution: 0.29 Rsun (2048³ cells, 5 levels AMR)
AMR implementation: set by hand to have max level around point particles
Max resolution zone: within min(5e12cm, 1.5*particle_separation)
Buffer zones: 2 cells
Softening length: 2.4 Rsun
Ambient density: 6.7e-9 g/cc
Ambient pressure: 10⁶ dyne/cm²
Default Accretion Routine=0 (No accretion)

Comparison between runs 133(top left), 135(top right), O+16a(bottom)

Comments:

• The results of 135 and 133 are similar, except:
  • Run 135 (no accretion) has slightly smaller separation at a given time.
  • Run 135 (no accretion) has slightly larger frequency of oscillations.
• Thus accretion cannot fully explain the discrepancy with the O+16a results, but removing accretion does make results slighly closer to those of O+16a, who did not have accretion.
• We must think more about accretion, but it is not making a huge difference at present.

Run 132 that uses twice as high resolution and 10x lower ambient pressure as previous runs:

Binary run 132 with double max resolution, lower resolution in ambient medium, 10x smaller ambient pressure than run 116
Relaxation run: 129
First frame: 75 (5 RG freefall times, when velocity damping ended)
Last frame:
Total simulation time:
Machine and partition: Stampde 1 normal (running, completed up to frame 215, or 32 sim-days)
Number of cores: 1024
Total wall time: 8 days (starting from frame 75 of relaxation run)
Hydro BCs: extrapolated
Poisson BCs: multipole expansion
Box size: L=4e13 cm (575 Rsun)
Base resolution: 2.25 Rsun (256³ cells)
Highest resolution: 0.14 Rsun (4096³ cells, 4 levels AMR)
AMR implementation: set by hand to have max level around point particles
Max resolution zone: within 5e12 cm (71.87 Rsun) of primary center and within a cylinder of radius 20 Rsun and height 20 Rsun around secondary center. After t~11d (~frame 123) refinement radius around primary was halved to 2.5e12cm (36Rsun). After t~31d (frame 210) halved again to 1.25e12cm (18Rsun).
Buffer zones: 2 cells
Softening length: 2.4 Rsun
Ambient density: 6.7e-9 g/cc
Ambient pressure: 10⁵ dyne/cm²
DefaultAccretionRoutine=2 (Krumholz)

NOTE: For the last 4 days of wall time the code ran very slow, probably because it was creating many many low mass particles due to the Jeans criterion. I have now commented out this part of the code, so it does not create particles. It is pending on stampede from chombo file 208, before extra particles were created.

Comparison of run 132 with Ohlmann+16a run:
​Table

Comments:

• The softening length of the point particles decreases with time in the O+16a simulation.
  • It is not clear whether this procedure is justified.
  • The main reason to improve the resolution with time may be to resolve the decreasing softening length.
  • Therefore, if we keep the softening length constant, improving the resolution with time may not be necessary or productive.
• O+16a initially resolves the softening length by 20 cells, whereas we resolve it by 17 cells. They advocate >10 cells.
• At 120 days, O+16a resolves the softening radius by 35 cells.
• It may be worth doing a run with slightly larger or smaller resolution to ensure that we are adequately resolving the softening length.

Comparison between runs 088(top left), 125(top right), 132(bottom left), O+16a(bottom right)

Comments:

• The 1st minimum, 2nd maximum, and 2nd minimum are located roughly at separations, with softening length denoted as h:
  • Run 088 (res=0.29Rsun, h=4.8Rsun, cells/h=17): 14 , 21 , 12 Rsun
  • Run 125 (res=0.29Rsun, h=2.4Rsun, cells/h= 8): 13.5, 18 , 14.5 Rsun
  • Run 132 (res=0.14Rsun, h=2.4Rsun, cells/h=17): 14 , 18.5, 11 Rsun
  • O+16a (see above table) : 10 , 23.5, 11 Rsun
• These are found roughly at times:
  • Run 088 (res=0.29Rsun, h=4.8Rsun, cells/h=17): 12.5, 15.5, 18 days
  • Run 125 (res=0.29Rsun, h=2.4Rsun, cells/h= 8): 12.5, 14.5, 17 days
  • Run 132 (res=0.14Rsun, h=2.4Rsun, cells/h=17): 12.5, 14.5, 17 days
  • O+16a (see above table) : 13 , 16 , 19 days
• From the difference between 125 and 132, we conclude that h is not adequately resolved in run 125 (8.5 cells/h vs 17 cells/h)
• O+16a advocates >10 cells/h, so this is consistent with their result.
• In our case, we still cannot tell if 17 cells/h is sufficient.
• We could run 088 with double the resolution to see if it converges.
• I will also continue to run 132 but it will be slow as the queue time in stampede 1 has increased as nodes are being transferred to stampede 2.

Run 136 to test how much relaxation makes a difference (IN PRORESS):

Binary run 136 similar to run 135 but now start from frame 0 of relaxation run 096 instead of frame 75 (so no damping).
Relaxation run: 096 (but no relaxation! Just start with initial profile
First frame: 0
Last frame: 18, running
Total simulation time: 4 days
Machine and partition: bluehive2.5 standard
Number of cores: 120
Total wall time: 21 hours, running
Hydro BCs: extrapolated
Poisson BCs: multipole expansion
Box size: L=4e13 cm (575 Rsun)
Base resolution: 9.0 Rsun (64³ cells)
Highest resolution: 0.29 Rsun (2048³ cells, 5 levels AMR)
AMR implementation: set by hand to have max level around point particles
Max resolution zone: within min(5e12cm, 1.5*particle_separation)
Buffer zones: 2 cells
Softening length: 2.4 Rsun
Ambient density: 6.7e-9 g/cc
Ambient pressure: 10⁶ dyne/cm²
Default Accretion Routine=0 (No accretion)

Discussion

• Possible reasons for discrepancy with O+16a in (my guessed) order of importance
  1. resolution (especially number of cells per softening length h) (a new run will test this).
  2. softening length.
  3. different RG relaxation runs with different type codes.
  4. initial spin of RG (a new run will test this).
  5. very different ambient densities (I have an old run that we can use to check this).
  6. different accretion algorithms.
  7. small differences in particle masses.
  8. some reflection off of boundaries at late times in our smaller box.
• Probably run 132 is at the highest resolution we're going to get for the paper (especially if we need to do multiple runs).
• But we need to make sure that the softening length h is adequately resolved by checking that results are converged.

Next steps

• Continue to run high res run 132 (with particle creation turned off!).
• Continue to run 136 (no relaxation run).
• Make separation plot for old run 120 to see if ambient density matters.
• New run with initial spin of RG at 95% Keplerian to test what difference this makes.
• New run like run 088 but with twice as high resolution to test the convergence with cells/h.
• Continue to improve and extend the outputting of relevant data to various files.
• Clean up code and divide into different modules for relaxation runs and binary runs.
• Think more about the accretion routine, Krumholz vs. Federrath vs. no accretion….

Drag analysis for Run 001:

Were able to test runtime for gamma = 1, 1.2, 1.4 & 1.67. All of them in the order of a couple hours & affordable.

Log density plot for Run 002:

X direction velocity plot for Run 002: