Planet Status
Blog post summarizing everything working so far for the parameter space study.
x10 flux
2x10^{17} flux (high)
Instead of using theoretical calculation, used theoretical relationship for 1D run - so upped planet density by only 2 orders of magnitude. Looks pretty good.
Radiation Pressure
This was before implementing ramping, without setting up the steady state wind first. Currently working on running from frame 150 of Run2 with ramping of the Lyman-alpha flux.
Update 2/19
Blog post summarizing everything working so far for the parameter space study.
Fixed bug in WASP-12 run - needs to restart from frame 10.
Highest-flux case working, but not launching wind. I expect the problem is resolving the ionization front:
Run1 flux, for comparison
Run3 flux
Run3 Lyman-alpha cooling
Run3 ionization heating
Run1 Lyman-alpha cooling
Run1 ionization heating
Finally, streak images are inverted:
Just invert Y-axis in matlab?
COMMON ENVELOPE SIMULATIONS
New Work
- Plotted v_phi/v_Kepler and v_phi/c_s in frame that is corotating with instantaneous particle orbit, centered on particle 2 (the companion).
- Plotted potential in the corotating frame with coordinate origin at particle 2.
- Added a section to "en.pdf" energy notes to include the potential in a rotating frame.
Plots of velocity relative to particle 2 in corotating frame
- Plots below show a slice through the companion, parallel to the orbital plane, with companion at the center and RG core at the left.
- Each particle is surrounded by a green circle with radius equal to the instantaneous softening length.
- Plots are in the frame corotating with the orbit (instantaneous value of ).
- Pseudocolor represents where both numerator and denominator are relative to the companion (first plots), or , with relative to the companion.
Run 143 (no accretion):
- v_phi relative to particle 2 in corotating frame normalized by Keplerian circular speed (Keplerian speed is uncorrected for spline potential):
- same plot but different color scheme:
- normalized by sound speed:
Run 132 (Krumholz accretion):
- normalized by Keplerian circular speed (Keplerian speed is uncorrected for spline potential):
- same plot but different color scheme:
- normalized by sound speed:
Plots of potential in corotating frame with origin at particle 2
- Potential does not include that of the gas (Roche potential).
- Potential is normalized to , where and is the interpartical distance.
- Note that is almost equal to from Kepler's 3rd law, but not quite since that would ignore the effect of gas on the particle orbit.
- Contours are at -3 to -1.25, in spacings of 0.25.
- Run 143 (no accretion):
- Run 132 (Krumholz accretion):
Updated notes on energy, with the last section just added
en.pdf
Next steps
- Edge-on plots of density and velocity.
- ADAF papers; read and discuss.
- Begin draft of paper.
Meeting update -- 02/12/2018
- Laurence Sabin's Visit
- Time
- Projects
The main goal would be using the hydrodynamical model that you already published with Bruce and add the "magnetic component" to fit our observations. It would also be interesting to combine this with Martin's (synthetic) polarization maps to reflect what we are observing with the SMA, CARMA and ALMA.
A second point, that was raised during the PNe conference Adam and I attended last year, is the determination of the minimum field's intensity needed to trigger the shaping. This is a very important and rather unknown aspect: I am working on the measurement of photospheric magnetic fields, on Post-AGBs and PPNe, and so far the (longitudinal) values found are quite low and might not be enough to actually launch any material !!
- Current JetClump Module: 2D MHD (Toroidal magnetic field in Clump and Jet ) runs OK. 3D MHD untested.
Update 2/12
Simulations
Blog post summarizing everything working so far for the parameter space study.
Started test of radiation pressure run - first few steps are struggling (pressure protections and high CFL). Once VisIt is done with mass-loss rates, will be able to see what's up.
WASP-12 run is at frame 116. Need to look at output - console only has level 0, but when I checked previously the planet was resolved at level 2. Also, execution time shot up at frame 116.
Last but certainly not least, adding a tracer to the ambient and stopping recombination for ambient material seems to be working fairly well at solving that problem (first frame in short time, no recombination in ambient yet), but the execution time still shot up - only got 3.5 frames in 5 days. Again, once VisIt is free I'll be able to take a look at the problem.
COMMON ENVELOPE SIMULATIONS
New Work
- Plotted separation vs time for runs 132 (Krumholz accretion) and 143 (no accretion) on the same plot, with orbits as insets.
- Plotted accreted mass and accretion rates for runs 132 and 143 up to t=40 days.
- Accretion radius in Krumholz model: discussion.
Orbital separation with time and orbit
- Below is the orbital separation as a function of time for run 143 (no accretion) in solid blue, and run 132 (Krumholz) in dashed light blue.
- The solid red (dashed gold) line shows the radius of the maximally refined region for run 143 (132). Note that for run 132, the maxlevel refinement was within a sphere centered on the primary, while for run 143, it was initially within a sphere centered on the primary, but from frame 72 (t=16.7 days), the center shifts to the secondary. As well, for run 132, maxlevel refinement is also done inside a cylinder of height 20 Rsun and radius 20 Rsun centered on the secondary, until the time just before frame 215 (32.4 days). (Note that run 132 starts from frame 75 of a damping run, while run 143 starts from frame 0 without any damping having been applied.)
- The solid green (dashed light green) line shows the softening radius for the spline potential for run 143 (132).
- Inset: The trajectories of the RG core (red/gold) and companion (blue/light blue) are shown for runs 143 (left) and 132 (right). Green/light green circles show the softening radius at t=0, and for run 143, also at time t=16.7 days, when the softening radius was halved.
Accretion
- Below is the total mass contained inside a sphere centered on the companion as a function of time for run 143 (no accretion subgrid model), for spheres of four different radii (blue) and the corresponding accretion rates, obtained by differentiating the interior mass (red).
- Below is the accreted mass as a function of time for run 132 (Krumholz accretion model) in blue, and the corresponding rate in red.
Discussion about accretion radius in Krumholz accretion subgrid model (with Bo)
- In the discussion following equation (17) in Krumholz+04, we read:
- The most important sentence is: "In general, the softening length should be smaller than the size of the accretion region, to ensure that softening does not alter the rate at which gas crosses its boundary." Is this true for our run 132?
- In our case, I (unwittingly) left the accretion radius set to the default value, which is 4 grid cells.
- The softening radius is about 16-17 cells, kept constant for each simulation. This was intentionally set to a large number >10 to ensure that the flow near the particle is adequately resolved (Ohlmann+16a).
- However, we used a spline potential, while Krumholz+04 is referring to a Plummer potential.
- As mentioned in a footnote on page 64 of Ohlmann16_thesis, a Plummer softening length of unity is effectively equivalent to a spline softening length of 14/5=2.8 if we choose the softening length such that the spline and Plummer potentials have the same value at r=0. (I have checked this to be true.) This implies that the effective Plummer softening length used is actually ~6 grid cells.
- Therefore, our effective Plummer softening length is still greater (by 50%) compared to the accretion radius.
- This means that the accretion rate would be artifically reduced relative to the Krumholz accretion model.
- The Krumholz accretion is (i) based on the Bondi-Hoyle-Lyttleton theory, which is not really applicable in our case (Macleod+17b) and (ii) ignores pressure feedback onto the flow. For both of these reasons it probably overestimates the accretion rate.
- Thus, our model which (accidentally) reduces artificially the Krumholz rate by having an accretion radius < Plummer equivalent softening radius, is conservative/cautious in the sense that whatever features we find (e.g. accretion) disk would be expected to be more pronounced (stronger) if the full Krumholz model was used. Because in that case the accretion radius would be larger than the softening radius and the accretion rate would be higher.
Next steps
- Plot v_phi/v_Kepler and v_phi/c_s in the frame corotating with the orbit, and adjust vectors and binding energy contours accordingly.
- Edge-on plots of density and velocity.
- Read ADAF papers.
- Start writing paper?
Update 2/11: Wind Tunnel Runs 004AB, 005AB, and comparison to 003
Finished Runs
Just to summarize the four recent runs -
WT Run 004A
- Sim parameters:
- Spacial resolution: base 64^{3}, finest 1024^{3} (4 levels)
- Time Resolution: 0.005CU = 2.2e3s
- Sim time: 0.5 CU, 100 frames
- Accretion = none
- Poisson solver convergence tolerance: 1e-6
- Run time: less than a day.
- Softening Length: 8CU (16 cells)
WT Run 004B
- Sim parameters:
- Spacial resolution: base 64^{3}, finest 1024^{3} (4 levels)
- Time Resolution: 0.005CU = 2.2e3s
- Sim time: 0.5 CU, 100 frames
- Accretion = Krumholtz
- Poisson solver convergence tolerance: 1e-6
- Run time: less than a day.
- Softening Length: 8CU (16 cells)
- Accretion Radius (particle radius): 0.5CU (1 cell)
WT Run 005A
- Sim parameters:
- Spacial resolution: base 64^{3}, finest 2048^{3} (5 levels)
- Time Resolution: 0.005CU = 2.2e3s
- Sim time: 0.25 CU, 50 frames
- Accretion = none
- Poisson solver convergence tolerance: 1e-3
- Run time: paused after 1.5 days. Estimated next 0.25 CU takes 4.5 days.
- Softening Length: 4CU (16 cells)
WT Run 005B
- Sim parameters:
- Spacial resolution: base 64^{3}, finest 1024^{3} (5 levels)
- Time Resolution: 0.005CU = 2.2e3s
- Sim time: 0.5 CU, 100 frames
- Accretion = Krumholtz
- Poisson solver convergence tolerance: 1e-3
- Run time: less than a day.
- Softening Length: 4CU (8 cells)
- Accretion Radius (particle radius): 4CU (8 cell)
Run 005A log Rho movie
Run 005A vx movie
Run 005B log Rho movie
Run 005B vx movie
Different definitions of softening lengths
The softening function used by both Luke and me is spline interpolation. As noted by Ohlmann, the softened potential is steeper compared to Plummer's potential. By comparing the potential at r=0, Ohlmann argues for a factor of conversion 2.8. Then in Plummer's term my softening length is within 3 pixels. Yet Ohlmann claims that a 10 pixel resolution is necessary.
Drag force: a first comparison
First guesses:
- Due to self-gravity, drag force does not asymtote to stable value?
- Dynamical friction dominates over momentum accretion?
Simulation Movie Summary
Physical Parameters
1.32443x10^{-16} g/cc | |
1.350 M_{sun} | |
(0.04747, 0, 0) AU | |
1.95 rad/day | |
lRotate | T |
lStellarGravity | T |
Simulation Parameters
BaseRes | (80,80,80) |
MaxLevel | 3 |
Zones/R_{P} | ~27 |
Extent | (37, -10, -10, 57, 10, 10) |
Physical Extent | (5.55x10^{11}, -1.5x10^{11}, -1.5x10^{11}, 8.55x10^{11}, 1.5x10^{11}, 1.5x10^{11}) cm |
TimeScale | ~8.3 hours |
lScale | 1.5x10^{10} cm |
lScale = R_{P}
rhoScale = rho_{P}
Run # | M_{P} (M_{J}) | Flux (phot/cm^{2}/s) | Status |
1 | 0.07 | 2x10^{13} | Complete |
2 | 0.263435 | 2x10^{13} | Complete |
2b | 0.263435 | 2x10^{13} | Running, frame 70 |
3 | 0.263435 | 2x10^{17} | |
4 | 0.07 | 2x10^{17} | |
5 | 0.07 | 2x10^{14} | Running, frame 22 |
6 | 0.263435 | 2x10^{14} | Running, frame 66 |
7 | 0.07 | 1x10^{15} | |
8 | 0.263435 | 1x10^{15} |
2b is higher-density planet surface with low flux.
Movies
noRot, Run 1 (side view)
noRot, Run 2 (side view)
Rot, Run 1 (top view)
Rot, Run 1 (side view)
Rot, Run2 (top view)
Rot, Run2 (side view)
Mass Loss Rates
noRot, Run 2 mass loss
noRot, Run 1 mass loss
Rot, Run 1 mass loss
Streak Streamline Images
Run1
Run2
Planet Closeups (top views)
Run1
Run2
Update 2/7: Wind Tunnel Run004 & 005
Finished Runs
Finished Run
WT Run 004A
- Unless otherwise mentioned, same setup with Run003.
- From now on self gravity is enabled. Compared to Run003, this slows down the sim within order of 2, as expected.
- mu set to 0.62, average mu of the Sun. This changed the temperature, but the pressure and the sound speed is not changed. And thus Bondi radius is also unaffected.
- Slight change of timescale - one wind passing time of the box.
- I did not understand the softening length properly: it's set in CU, thought it was set in # of cells. So previously I have been using 10CU=20 cells. For this run I'm using 8CU = 16 cells. Will discuss this a bit more later.
- Simulation time = 0.5CU = 2.57 days
- Wind temp = 3.45e5K
- Sim parameters:
- Spacial resolution: base 64^{3}, finest 1024^{3} (4 levels)
- Time Resolution: 0.005CU = 2.2e3s
- Accretion = none
- Hydro solver BC: extrapolated
- Self graivty Poisson solver BC: Periodic
- Poisson solver convergence tolerance: 1e-6
- Run time: less than a day.
Run 004A logRho plot movie: Run 004A Mach number movie:
WT Run 004B
Setup is identical with Run 004A except that Krumholtz accretion is turned on. The radius of the particle (accretion radius) is set to 1. That is the main flaw of this run: the accretion radius is deep within the softening length, meaning the density sampled for the Krumholtz procedure to estimate the Bondi accretion rate is inaccurate. This motivated Run005s': to have the accretion radius be at least as large as the softening length. Ideally of course accretion radius should be larger, so that the correctly calculated densities of surrounding gas can be sampled. Within the softening length, without graivity, pressure should become equalized, meaning the density profile would be flatter than what we expect realistically. This would surely impact the performance of Krumholtz procedure.
Due to space constraints I did not make movies of 004B, though it can be easily reproduced.
WT Run 005A
- This is anagolous to Run004A, but with a factor of 2 higher resolution. This enables the softening length to be reduced to 4CU, while still resolved by 16 pixels. Bondi radius is around 18 CU.
- To save time Poisson solver's tolerance is set to 1e-3.
- After running for 1.5 days, first 50 frames (0.25CU) was produced. The remaining wall time is at 4.2 days and increasing. I thus decided to end the simulation. As a proof-of-concept this resolution would require an order of magnitude increase in computational resources, which is not unfeasible on Bluehive, but if possible to run it on a better computer, it would be great. Although a larger box is needed for a steady state to be reached.
COMMON ENVELOPE SIMULATIONS
New Work
- As last post but now for run 132 (Krumholz accretion).
New Plots of run 132
The binding energy contours are for the binding energy relative to the companion.
The contours are -0.5 (grey), -0.25 (light grey), 0 (white), 0.25 (thin line, light grey).
That is, (i) the kinetic energy density is calculated using the velocity in the frame of the companion instead of the box frame,
and (ii) the RG core is not included in the calculation of the gravitational potential energy density.
- The normalized binding energy, with extrema at -1 and 1, is given by:
- Below are the same plots but with an unsaturated color scheme:
- Below are the same plots but zoomed out, with slightly different scaling for vectors:
- Below are color plots of the binding energy relative to the companion, since color is easier to understand than contours, at different zoom levels:
- Below: Very zoomed out version, with left plot showing binding energy relative to the companion and right plot showing binding energy relative to lab frame (copied from the last post):
- Below: Same view as above but showing density (left) and temperature (right), for comparison
- Below: Density at t = 0, 10, 20, 30 and 40 days
Next steps
- Separation vs time plot for both runs on the same graph.
- Plot(s?) showing orbits of both runs.
- Polish accretion rate plots for both runs and redo plot for run 132 with higher time resolution.
COMMON ENVELOPE SIMULATIONS
New Work
- Explored the case of binding energy relative to the companion, and plotted several figures.
- Snapshots of density at various times.
New Plots of run 143
The plots below are similar to those from the last post but with one main difference:
the binding energy contours are for the binding energy relative to the companion.
The contours are -0.5 (grey), -0.25 (light grey), 0 (white), 0.25 (thin line, light grey).
That is, (i) the kinetic energy density is calculated using the velocity in the frame of the companion,
and (ii) the RG core is not included in the calculation of the gravitational potential energy density.
- The normalized binding energy, with extrema at -1 and 1, is given by:
- Below are the same plots but with an unsaturated color scheme:
- Below are the same plots but zoomed out, with slightly different scaling for vectors:
- Below are color plots of the binding energy relative to the companion, since color is easier to understand than contours, at different zoom levels:
- Below: Very zoomed out version, with left plot showing binding energy relative to the companion and right plot showing binding energy relative to lab frame (copied from the last post):
- Below: Same view as above but showing density (left) and temperature (right), for comparison
- Below: Density at t = 0, 10, 20, 30 and 40 days
Next steps
- Produce the same figures but for run 132 which had Krumholz accretion.
- Separation vs time plot for both runs on the same graph.
- Plot(s?) showing orbits of both runs.
- Polish accretion rate plots for both runs and redo plot for run 132 with higher time resolution.
Update 2/5
Simulation Status
WASP-12b, w/ stellar rotation | |
Running | 75/150 frames complete, ~3.5 days remaining |
Rotating Frame | |||
Run # | M_{P} (M_{J}) | Flux (phot/cm^{2}/s) | Status |
1 | 0.07 | 2x10^{13} | Complete |
2 | 0.263435 | 2x10^{13} | Complete |
3 | 0.263435 | 2x10^{17} | Unqueued |
4 | 0.07 | 2x10^{17} | Unqueued |
Non-rotating Frame | |||
Run # | M_{P} (M_{J}) | Flux (phot/cm^{2}/s) | Status |
1 | 0.07 | 2x10^{13} | Complete |
2 | 0.263435 | 2x10^{13} | Complete |
3 | 0.263435 | 2x10^{17} | Unqueued |
4 | 0.07 | 2x10^{17} | Unqueued |
Run 3 testing
Looks like we need to take a different tack. The ambient is still ionizing too quickly and creating a shield in front of the planet. Lowering the flux (and therefore the density requirements) is one option - to get an idea of how much, look at the recombination timescale:
For
(one unit of time), need ~5x10^{-15} g/cm^{3} of (fully ionized) hydrogen. With a flux of 2x10^{17}, the ambient medium is at about 10^{-14} g/cm^{3}. So reducing just one order of magnitude, to a flux of 2x10^{16}, may be sufficient to solve the ambient recombination problem, since it would increase the recombination timescale to about 10x the simulation timescale.Another option is to consider other adjustments to the ambient profile. If we were able to start the exponential profile farther in and still get a sufficiently low density at around R_{p}, that seems like a good option, but I can't seem to fiddle the profile into such a shape.
Run2 side view
Run2 top view
Other analysis products?
- Movies (as above)
- JC’s streak images that show flow pattern?
- Total Mass loss rates
- Back flow rates?
- Observations of the kind in Carroll-Nellenback ea 2016?
COMMON ENVELOPE SIMULATIONS
New Work
- Plotted circles for softening length around particles and new script to get softening length from chombo into VisIt.
- Calculated binding energy and plotted contours.
- Wrote general script that plots any pseudocolor, contours, and vector arrows along with softening circles.
- Polished plots and included the option of plotting in units of solar radii instead of cm.
Run 143:
- Similar to run 132 but with subgrid accretion turned off, twice larger box, resolution and softening length that evolve with time, less aggressive refinement, i.e. larger max refinement zones, and no relaxation (damping) run initially
Relaxation run: no relaxation run
First frame: 0
Last frame: 173
Total simulation time: 40 sim-days
Machine and partition: Stampde 2 normal
Number of nodes: 128 or 64, each with 68 (standard nodes) or 96 (skx nodes) cores
Total wall time: TBD
Hydro BCs: extrapolated
Poisson BCs: multipole expansion
Box size: L=8e13 cm (1150 Rsun)
Max refinement level: 4 (frames 0 to 72), then 5 (frames 72 to 117)
Base resolution: 2.25 Rsun (256^{3} cells)
Highest resolution: 0.14 Rsun (4096^{3} cells, 4 levels AMR, frames 0 to 72), and 0.07 Rsun (8192^{3} cells, 5 levels AMR, frames 72 to 173)
AMR implementation: set by hand to have max level around RG core (frames 0 to 72) or companion (frames 72 to 173)
Max resolution zone: sphere around RG core, radius 5d12 cm (frames 0 to 46), radius 4d12 cm (frames 46 to 72), radius 3d12 (frames 72 to 103), radius 2.5d12 (frames 103 to 161), radius 1.75d12 (frames 161-173)
Buffer zones: 16 cells
Softening length: 2.4 Rsun (frames 0 to 72), 1.2 Rsun (frames 72 to 173)
Ambient density: 6.7e-9 g/cc
Ambient pressure: 10^{5} dyne/cm^{2}
DefaultAccretionRoutine=0 (no accretion)
New Plots of run 143
- All plots below show a slice through the companion, parallel to the orbital plane, with companion at the center and RG core at the left.
- Each particle is surrounded by a green circle with radius equal to the instantaneous softening length.
- Pseudocolor represents where both numerator and denominator are in the frame of the companion (first plot), or ,
with
in the frame of the companion (second plot).- Velocity vectors shown in the frame of the companion, with length scaled to vector magnitude.
- Contours show the values -0.5 (dark grey) and -0.25 (light grey) of the quantity
is the potential energy density, including self-gravity of the gas (which comes with a factor of ½ multiplying the gas potential to avoid double-counting), potential energy of gas-RG core system, potential energy of gas-companion system, but NOT potential energy of RG core-companion system. In other words, includes all potential energy EXCEPT for particle-particle term. Note that the potential being used for the particles is the actual spline potential, which had to be entered by hand. Also is the bulk kinetic energy density in the lab frame, while is the internal energy density in the lab frame. Note that the script also plots the contours -0.75, 0, 0.25, 0.5, 0.75, but they do not appear on the graph so the smallest value is >-0.75 while the largest is <0 for this part of the simulation box.
- Below is a similar plot to the first plot, showing , now zoomed out (with vectors and colors scaled differently).
- The contours for binding energy are -0.75 (dark grey), -0.5 (grey), -0.25 (light grey), 0 (white), 0.25 (dashed light grey).
- Below are color plots of the binding energy, since color is easier to understand than contours.
- Below is a zoomed out version.
- Below is a more zoomed out version, different color scheme.
- Below is a very zoomed out version, with two plots showing two different color schemes. Note that box was rotated so that the RG core was situated to the left of center.
Next steps
- Further polish the above plots (more contour levels? improve resolution of vector arrows?).
- Same plots for run 132 which had Krumholz accretion submodel.
- Separation vs time plot for both runs on the same graph.
- Plot(s?) showing orbits of both runs (possibly as inset on separation vs time graph)
- Plots showing face-on slice and edge-on (slice/projection?) of snapshots of density at various times.
- Polish accretion rate plots for both runs and redo plot for run 132 with higher time resolution.
- Line Integral Convolution using Jonathan's MatLab code.