Comparison with AMR
Both with recombination/cooling, 2e13 flux, no Lyman alpha cooling
323 + 2 levels AMR
1283 fixed grid
Lineouts
Circles are 323 + 2 levels AMR, triangles are 1283 fixed grid. The 1283 fixed grid is 200 frames shorter - the last frame is repeated. Day side looks nearly identical - night side differs, probably due to lower resolution.
Density
Ionization fraction
Optical depth
Temperature
Heating/cooling time scales
Heating/cooling rates
Planet To-do
- Add diagnostics for cooling and recombination. ✓
- Subcycling time for planet. ✓
- Add flag for Lyman cooling. ✓
- Check cooling terms for correctness. ✓
- Ramp up flux over time. ✓
- Check recombination on night side. ✓
- With stored optical depth, calculate ionization and ionization heating rate (can't be done straightforwardly with a diagnostic variable). ✓
- Jonathan will be adding optical depth variable. ✓
Runs
- Lyman cooling, recombination/cooling, 2e13 flux ✓
- Recombination/cooling, 2e13 flux ✓
- Recombination/cooling, no flux ✓
- Lyman cooling, recombination/cooling, 2e10 flux ✓
- Calculate new profile for planet unaffected by Lyman cooling ✓
- Lyman cooling, recombination/cooling, ramping flux 1.2e13 → 2.02e13 ✓
- Recombination/cooling, ramping flux 3.6e12 → 2e13 - running
Comparison with O+16a, attempt 2
Recap of Last Post
- In the last post, I performed a run with the same parameters as the Ohlmann+16a run.
- However the simulation went from 2-body to 3-body when the secondary was mistakenly reintroduced after about 7 days.
- This was not a bug in the code, just me forgetting to flip a 'switch'.
New Work
- I flipped the switch and reran the code up to 20 days.
Summary of New Results
- The results match qualitatively those of O+16a. Differences are likely caused by remaining differences in the numerical setups, including softening radius, ambient density, resolution, and degree of corotation.
Detailed Results
I) Circular orbit as in Ohlmann16a
Damp088) Extrapolated hydro BCs, Multipole expansion Poisson BCs, ambient dyne/cm
(comet compute 1728 cores, 2 cores per task to increase memory per task)
( cm, , 5 levels AMR)
Restarted from run Damp062, at s, just after Damping stopped, to s
2d density
2d density (zoomed)
2d density and velocity (zoomed)
2d density edge-on slice (zoomed)
2d particle velocity (zoomed)
2d particle position (zoomed)
2d Particle mass (zoomed)
Here are PROJECTIONS, instead of slices, zoomed in with a better color scheme, and extended to go up to 20 days:
2d density projection
2d density projection edge-on
Here is a version of the SLICE in the orbital plane, going up to 20 days, rotated 180 degrees to match the orbit of O+16a, color range shifted to match O+16a, and with a better color scheme:
2d density
…or with blue labels for better visibility:
2d density
2d density (zoomed)
Here is temperature:
2d temperature
…and Mach number:
Mach number
Below I've plotted Fig 1 of Ohlmann+16a and the equivalent figure with this simulation for comparison. In my plot of separation, I've shown the softening radius as a dotted horizontal line. In the inset showing particle trajectories, I've inverted the axes and shifted the origin to allow direct comparison (our simulation is rotated by 180 degrees with respect to theirs). Circles represent the softening radius, while the green dot shows the initial center of mass of the two stars. The very small green square in the upper right shows the smallest resolution element.
Below are snapshots of density in the orbital plane from t = 10 days and t = 20 days. Note that to compare directly, I have rotated each figure by 180 deg. In the O+16a figure, the '+' denotes the primary core particle, while the 'x' denotes the secondary point particle. The core and secondary are denoted as '0' and '1' respectively in my plot.
Discussion
- The results match qualitatively those of Ohlmann+16a, but differ in detail.
- The orbit of the secondary is larger, while that of the primary is slightly smaller, compared to O+16a. The separation does not become as small at the first minimum or as large at the first maximum, compared to O+16a. This may be due to the ~5 times larger gravitational softening length used.
- At t = 10 days, the density in the orbital plane appears similar to that of O+16a, but the spiral arm is less pronounced and features are less detailed. This could be due to the resolution being lower by a factor of somewhere between 2 and 12. It may be partly caused by the weaker gravitational interaction due to larger softening radius mentioned above.
- Also, the common envelope has not expanded as much into the ambient medium. This could be due to the >107 times larger ambient density used compared to that of O+16a.
- Recall the differences between our setup and O+16a, outlined in the Table 2 of the following pdf file. The main differences between the two setups are:
- 5x larger softening radius
- 107 times larger ambient density
- up to 12 times smaller resolution
- 9x smaller box size
- no initial rotation vs 95% solid-body corotation with orbit
- some differences in the relaxation run to prepare the red giant star (detailed in the table and in previous posts)
- extrapolated instead of periodic BCs
Next Steps
- The easiest and most obvious thing to do is to reduce the softening radius (= cutoff radius defening RG core. Recall that for r smaller than this radius, we replace the core with a modified Lane-Emden solution). However, reducing this softening radius to 1 Rsun and keeping the same max resolution would mean that we resolve the softening radius by only 3.5 resolution elements, whereas O+16a claim that 10 resolution elements is needed.
- That suggests that we should increase the max resolution by a factor of 2 or 4, but this would slow down the code considerably. It would be best to first lower the softening length without changing the resolution.
- It would also be worth testing lower ambient densities. In past tests the code has crashed, but it should be possible to lower the ambient density by at least a factor of 10 or 30.
- So this suggests the following 3 runs (starting from the beginning of the relaxation run):
- As 088 but reduce the ambient density by 10 or 30
- As 088 but reduce the softening radius by 4.8 to equal that of O+16a
- As 088 but with both of the above changes
Reproducing Ohlmann+16a results: 1st attempt
Recap of Last Post
- In the last post, I performed several runs with varying separation and secondary mass.
- However it was determined that point particles were not feeling the gravity of the gas, so the results were wrong.
New Work
- I wrote a script in IDL which calculates the Kepler orbit given the initial velocities, and used it to show that the orbit obtained was just what would be expected if the point particles were not feeling the gas gravity, further confirming the above conclusion.
- This problem with the code was corrected by Baowei and Jonathan.
- The base case (run 084 from last blog) was rerun now with the code corrected, with the goal of reproducing the basic results from Ohlmann+16a.
Summary of New Results
- The simulation now looks "reasonable" at the beginning (first ~dynamical time) but then deviates from the results of Ohlmann+16a.
- I realized while writing this blog that the problem is caused by the mistaken re-introduction of the secondary (hence the introduction of a third particle with mass equal to that of the secondary) on restarts. In fact, there was a switch in the problem.data file called "secondary_already_present" to deal with restarts, but I had forgotten to change it to .true. before restarting.
Results
I) Circular orbit as in Ohlmann16a
Damp088) Extrapolated hydro BCs, Multipole expansion Poisson BCs, ambient dyne/cm
(comet compute 1728 cores, 2 cores per task to increase memory per task)
( cm, , 5 levels AMR)
Restarted from run Damp062, at s, just after Damping stopped, to s
2d density
2d density and velocity
2d density and velocity (zoomed)
2d temperature (zoomed)
2d particle velocity (zoomed)
2d particle position (zoomed)
2d Particle mass (zoomed)
Particle trajectories, (with circles representing softening radius)
Particle separation with time
Particle trajectories up to mistaken introduction of third particle
Particle separation with time up to mistaken introduction of third particle
Discussion
- What has happened is that another secondary got introduced on each restart. This happens at frame 107. It also happens at frame 92 but this won't matter much because I had to redo that frame only on bluehive because there had been an I/O problem at that frame on comet. So the simulation is correct up to frame 106, or elapsed time t = 2.12e6 - 1.5e6 = 6.2e5 s = 7.2 days, about 1 dynamical time, but is wrong after that.
- Results up to the kink seen in the plot of separation vs time are very close to those of Ohlmann+16a (see below). The introduction of the third point particle occurs before the kink. So the results so far are consistent with those of Ohlmann+16a. (If one looks closely there is another point in the plot of separation vs time that is a little off, and this is due to the one frame that was done on bluehive where the 3rd particle also appeared.)
Next Steps
Once I've edited the problem.data file to make "secondary_already_present=.true.", I will redo the simulation from frame 106. I will then again compare the result with Ohlmann+16a
Figure 1.
Here are the other figures from that paper for reference:
Figure 2
Figure 3
Figure 4
Update on common envelope simulations
Recap of Last Post
- In the last post I presented a simulation of a RG translating across the grid.
- I also presented a common envelope simulation but later realized that the initial velocities of the primary and secondary were wrong (about a factor of 2 too small). I corrected the script to calculate the initial velocities and tested it against the Sun-Earth system.
New Work
I corrected the initial velocities and simulated the following cases:
- RG star of mass 1.956 Msun including 0.369 Msun core (point mass), with secondary point mass of 0.978 Msun (half the RG mass), at initial separation 49 Rsun, just outside of the RG outer radius of about 48 Rsun. Very close to parameters of Ohlmann+16a.
- As above but larger separation of 60 Rsun instead of 49 Rsun.
- As above but even larger separation of 98 Rsun instead of 49 Rsun.
- As the first case but now make secondary mass equal to only 0.001 Msun so that problem reduces to "1-body".
- Replace RG with a point particle of the same mass and run in low resolution with uniform background to test orbital dynamics.
- Same but with original point masses of the first case, i.e. the RG core and the secondary, without RG envelope.
Summary of New Results
- We do not achieve a circular orbit in the CE sims. In each case the orbit is looser than it should be. The binary separation grows with time and the speed of the secondary reduces with time.
- When the RG is replaced by a point mass then the correct circular orbit is obtained (in low res with a uniform ambient medium).
- If the RG envelope is excluded and the point masses representing the core and secondary are made to orbit in a uniform ambient medium, we obtain almost the identical (non-circular) orbit as for the full CE simulation.
- The natural explanation for these results is that the point masses can feel each other's gravity, the gas can feel the gravity of the point masses, but the point masses do not feel gravity from the gas! This feature must obviously be included in the next set of runs.
Results
I) Circular orbit as in Ohlmann16a
Damp084) Extrapolated hydro BCs, Multipole expansion Poisson BCs, ambient dyne/cm
(comet compute 1728 cores, 2 cores per task to increase memory per task)
( cm, , 5 levels AMR)
Restarted from run Damp062, at s, just after Damping stopped, to s
2d density
2d density and velocity
2d density and particle velocity
2d density edge-on
2d Temperature
2d Particle position
2d Particle mass
II) As above but wider orbit (binary separation of 60Rsun instead of 49Rsun)
Damp086) Extrapolated hydro BCs, Multipole expansion Poisson BCs, ambient dyne/cm
(bluehive standard 120 cores)
( cm, , 5 levels AMR)
Restarted from run Damp062, at s, just after Damping stopped
2d density
2d density and velocity
III) As above but even wider orbit (binary separation of 98Rsun instead of 49Rsun)
Damp087) Extrapolated hydro BCs, Multipole expansion Poisson BCs, ambient dyne/cm
(comet compute 1728 cores, 2 cores per task to increase memory per task)
( cm, , 5 levels AMR)
Restarted from run Damp062, at s, just after Damping stopped
2d density and velocity
IV) As above but very low mass secondary (0.001Msun instead of 0.978Msun)
Damp085) Extrapolated hydro BCs, Multipole expansion Poisson BCs, ambient dyne/cm
(stampede normal 1024 cores)
( cm, , 5 levels AMR)
Restarted from run Damp062, at s, just after Damping stopped
2d density andparticle velocity
V) As above but RG represented by a point mass, so two point masses orbiting on uniform background
Kepler083) Extrapolated hydro BCs, Multipole expansion Poisson BCs, ambient density and pressure set to 1d-10
( cm, , 0 levels AMR)
(Started from t=0)
2D projected particle positions
particle velocities
2D projected particle positions, edge-on
VI) As (V) above but only mass of RG core is included in primary point mass, so same point masses as in
Kepler083core) Extrapolated hydro BCs, Multipole expansion Poisson BCs, ambient density and pressure set to 1d-10
(stampede normal 1024 cores)
( cm, , 0 levels AMR)
(Started from t=0)
2D projected particle positions
particle velocities
Comparison between and (VI)
Particle velocities
COMMENTS: Clearly, the secondary does not feel the gravity of the RG envelope. Thus, point masses do not feel gravity of the gas. This obviously needs to be corrected for the next set of runs.
Next Steps
Once this problems is resolved, we need to redo run in particular. If everything is in order, we can then plot the orbital separation as a function of time, and compare the result with Ohlmann+16a
Figure 1.
Here are the other figures from that paper for reference:
Figure 2
Figure 3
Figure 4
Aside from resolution and other details, the biggest difference between our setup and that of Ohlmann+ is that their RG is initiated with 95% corotation with the orbit, whereas ours is not rotating initially.
Other points worth mentioning
- My quota on bluehive has been reduced from 24 TB to 16 TB. I will look to reduce it more in the next two weeks. My quota on blue streak has been reduced from 5 TB to 2 TB.
- There is a nice paper by MacLeod et al. about accretion onto the secondary during the common envelope that I'm reading at the moment. Aside from being interesting, it is quite pedgogically written. I'm doing a literature review to prepare for the upcoming conference.
Update 6/14
- Implemented multispecies gas for Roe solver in AstroBEAR, but need a solution for density or sound speed for van der Waals gas to implement it. Relevant code snippet:
! Entropy fix !!! Implement a different entropy fix to use van der Waals gas - density equations have complicated/no solution lambdaLL = WL(ivx) - PrimSoundSpeed(WL) lambdaRR = WR(ivx) + PrimSoundSpeed(WR) if (lambdaLL < 0 .or. lambdaRR > 0) then call StarPU(WL((/ irho, iE, ivx /)), WR((/ irho, iE, ivx /)), PrimSoundSpeed(WL), PrimSoundSpeed(WR), pstar, ustar) ! Passively advected (tracers; y,z momentum) ExactL = WL ExactR = WR ! Solution in interaction region ExactL(iE) = pstar ExactL(ivx) = ustar ExactR(iE) = pstar ExactR(ivx) = ustar ! Densities differ if (mpi_id == 0 .and. iEOS == EOS_vanDerWaals) write (*,*) 'This solution for density will not be accurate. The attraction parameter is not yet included in these calculations.' ExactL(irho) = PrimExactDensity(WL,pstar) ExactR(irho) = PrimExactDensity(WR,pstar) lambdaLR = ustar - PrimSoundSpeed(ExactL) lambdaRL = ustar + PrimSoundSpeed(ExactR) end if if (lambdaLL < 0 .and. lambdaLR > 0) then ! Check for transonic rarefaction on left side print *, 'Using entropy fix.' lambda(1) = lambdaLL*((lambdaLR - lambda(1))/(lambdaLR - lambdaLL)) end if if (lambdaRL < 0 .and. lambdaRR > 0) then ! check for transonic rarefaction on right side print *, 'Using entropy fix.' lambda(nDim+2) = lambdaRR*((lambda(nDim+2) - lambdaRL)/(lambdaRR - lambdaRL)) end if
- Wrote WAF scheme for own solver. Working, but results are not as good as I expect based on Toro's. Not yet sure where the problem lies.
- Related: Transmissive boundary conditions for second set of ghost cells don't quite make sense yet. Why move back inward to get wave to flow out?
- Finished edits of WASP-12 paper; reprocessed absorption.
- Still getting explosion (but better than implosion) with adiabatic atmosphere.
Should be the same as Jonathan's result - not sure why it appears to not be working.
CE
In this simulation, I placed an outflow sphere at the center. The boundary condition is 30000K, 1e-13 g/cc and 300km/s from 0 to 100 day (mass loss rate is 7.2e-6 solar mass / yr) at the radii of 50 R_sun then the outflow speed changes to 10km/s. This outflow is not doing much at this timescale.
I launched another outburst from 1000 to 1100 day with 30000K, 1e-10 g/cc and 1000km/s. Thus the mass loss rate is 0.024 solar mass / yr.
Stable Atmosphere
I've added some data and a Matlab script for creating images such as those below. Density is on a linear plot- looks pretty good overall, although there's a strange wiggle near the inner boundary. Pressure is on a semilog plot, and is orders of magnitude off (I believe these should be in the same (computational) units).
EDIT: After running for longer time, not as stable as it first appeared (although I'm not seeing the same values in the lineouts as before, either). Perhaps the outer boundary is not as balanced as it should be?
Density plot (same scale as below):
Lineouts:
Pressure lineout scale changes dramatically after first frame, but can't get VisIt to save window correctly.
Meeting Update --6/2/17
- Disk Space
- archived MachStems data. Currently on bluehive, /scratch/afrank_lab has 4.4TB available.
- WireTurbulence
- redid Mach plot in blog:bliu05152017
- High res runs with tracers on Bluestreak: hydro 104/200 frames done, mhd 76/200 frames.
- Helped Bo working on his WindTunnel & StarAmbient modules. Some of the test results: wind tunnel; star ambience test1; star ambience test2
Update 6/2
- Roe solver implemented in AstroBEAR:
Still need to implement van der Waals gas and figure out how to implement multispecies (and find tests).
- Testing charge exchange:
Seems to be working pretty much as expected, although the proportions of cold HII to hot H seem a little wonky - likely because of the other physics present, though.
- WASP-12 paper nearly there. Got some comments back from Luca, so I'll look at those in the next few days.