Update 9/28
Link to pdf (might ask to be logged into UR email in the same browser): https://bit.ly/3jfJmda
Coupled EBM Project Update 09/26
Outline
- Introduction
- Numerical Model
- Modeling Earth
- Analytic Model
- Dimensionless Model
- Numerical Experiment Setup
- 5 Trajectories
- Full Suite of Results
- Parameter Sweeps
- KDE plots of decline time for different dT's
- Discussion
Questions
- Formatting: for emulateapj LaTeX document class, does anyone know how to have a figure on the top/bottom of the page that spans both columns?
- When talking about CO2, should I say P(CO_2) or pCO_2?
- For our timescales, call them…
- t_c, t_g, t_t?
- Or t_C, t_G, t_T?
- Or t_clim, t_grow, t_tech?
ToDo
- Finish the dT sweeps for the new version of the model, currently running.
- Finish filling in the details of section 1-4.
- Finish writing up section 5.
- Start/Finish writing up section 6.
Current Version of the Paper is Here
Dust in AstroBEAR - Update 2020/09/28
Objectives
- Debugging + visualising
Progress:
- Gas Drag: No progress on gas drag because I was busy with other stuff. Here's a quick recap of the current state of the routines: Both plasma and collisional drag are implemented however I think that the plasma drag is too strong so I need to check that. Could be related to the grain charge routine (which I suspect is buggy and might be the cause of several issues I see in other routines as well). Here's the lastest video showing a comparison between the different drag routines.
- Sputtering: No progress on sputtering either because I was busy with other stuff. Here's a quick recap of the current state of the routines plus a new plot made with the analysis script (see below): Sputtering routines are fully implemented but need debugging. I think that currently, the growth vs destruction rates are not right, probably because of the grain charge routine. Here is a recent video showing the destruction of grains during a simulation. With the analysis script, I can also get the total dust mass over time in all bins (UPDATED PLOT) and the dust mass per bin over time. Both I think progress too quickly and are too smooth which again points to destruction/growth rates being off.
- Dust Collisions: Grain-grain collisions are now fully implemented (see here and here for theoretical background). The code currently crashes after a couple of timesteps because the collision probability starts to exceed realistic thresholds. Might be due to the grain charge routine which I'm planning to start debugging next.
- Gas Feedback: Gas feedback is done except for gas drag (which is also "done", I'm just not too happy with it and would like to redo that at some point). There's now a gas feedback flag that can be set in the physics.data input file which turns the gas feedback either on or of. If the gas feedback is off, the gas will affect the dust (depending on other flags set) but the gas will not know that there is any dust in the domain. If gas feedback is on, the gas and the dust will affect each other (which I think no previous study has been able to do within the MHD sim.). The following is currently included in gas feedback:
- for the gas drag, the energy transferred to the dust grains as they are accelerated due to interactions with the gas atoms is removed from the gas energy field
- for sputtering, the dust that is sputtered to very small sizes is assumed to become part of the gas. If gas atoms are accreted, they are removed from the gas budget in the active cell. This means the overall mass is conserved. (If gas feedback is off, the "destroyed" dust is stored in the lower collector bin. Mass is therefore not conserved in the simulation but we can check the lower collector to make sure all lost mass is accounted for.)
- for ion trapping, the gas atoms that are trapped in the dust are removed from the gas. If this is not done, the dust would just keep trapping more and more gas atoms and continue to grow.
- for grain-grain collisions the grains broken into really tiny pieces are assumed to become part of the gas as well just like for sputtering.
- Analysis Scripts: I've finished the first version of my analysis scripts that will be used to visualise the dust mass evolution in the simulations. It's a bit rough and plots don't look too neat but it will now be easier to understand the overall dust content. The script requires little user input and gets all information from the initial data files and the output. It'll be part of the dusty AstroBEAR version so people using it will have access to it.
Next steps:
- Debugging
- Improve analysis script
Update 09/21
Theoretical Predictions are Off
Got some initial results from Mach and beta parameter sweep. The theoretical equation doesn't fit the simulations very well. We will probably include this in subsequent papers.
http://www.pas.rochester.edu/~aanand/movies/fall20/
Fixed the plot from last time.
EMAC
Need to fix some links on AstroBEAR website or make the code publicly available before upload.
COMMON ENVELOPE SIMULATIONS
EoS stuff
Estimate of resource use
- I'm getting roughly 1 frame per 6 hours on 32 nodes for the test run closest to what I'm planning.
- Let's assume that the average frame rate during the course of the whole simulation (as the separation reduces) is double that, so 1 frame every 3 hours or 96 node-hours per frame.
- Since the softening radius is small, we can go up to at least 300 frames. So we are talking 28,800 node-hours per run.
- We want to do at least 4 runs (not counting any single star run). That means 115,200 node-hours.
- That is 51% of our total stampede allocation.
- The code is almost ready and there is ample time to complete these runs before the allocation expires in April.
- We have so far used 27% (15% planets project, 6% jets, 6% eos test runs).
- This would leave only 22% for other things (radiative transfer, convection, Roche lobe overflow, planets)
- Since the proposal was for the CE project, I suggest we put a hard cap on non-CE related projects at 20% of the allocation.
- If we submit the supplement proposal and it is approved for a significant amount of resources, then obviously we re-evaluate.
Work done
- More investigation of core stability
- More investigation of surface stability
- Experimentation with refinement
- Various other small improvements to the code
- Tracers
- New EoS Notes (in progress)
- New EoS tables (Yisheng)
Results
- The core gains 2% during the first frame (0.23 d) but is remarkably stable over the next three frames (Run 221), fluctuating at a level of 0.2%. For comparison, the sound-crossing time of the core is ~0.01 days.
- This 2% gain is not caused by the EoS, because I performed a run (Run 229) which is the same as 221 but for the ideal gas EoS with gamma=5/3 and the same behaviour is seen.
- It is probably caused by a combination of the following:
- discretization effects
- limited resolution in MESA
- error when interpolating MESA profile into astrobear
- inability to fully resolve the very centre of the modified Lane-Emden solution, which causes a slight reduction in the peak pressure while hardly affecting the potential energy since the latter is dominated by the particle contribution (this effect goes in the right direction but can only account for a small part of the 2%)
- interpolation error at the softening radius where I matched the modified Lane-Emden solution with the MESA profile
- slight mismatch in dP/dr at the softening radius because we chose to match rho and drho/dr as the two conditions to fit the two parameters (as in Ohlmann 2017) (and the pressure is also forced to match as an integration constant)
- possibly something to do with producing frame 0 on bluehive and then restarting from this frame on stampede (I will avoid doing this in the next stampede run, but I doubt this explains it).
- discretization effects
- For both the core and the surface, I find that the minimum number of resolution elements per pressure scale height is EIGHT. To be precise it is somewhere between 5.1 and 8.2 but let's use 8 to be conservative.
- If ~4 is used at the core, the peak pressure rapidly loses its shape. I found that to resolve the core at 8.2 cells per scale height, the maxlevel 4 refinement region should start outside 9e10 cm (1.3 Rsun). Currently (e.g. Run 221) I am using 1e11 cm (1.4 Rsun) with 16 buffer zones, so more than adequate. This core region only accounts for about 3e6 cells, or about 1% of the total cells, which suggests there is no real need to skimp here. However, doubling it to 2e11 cm does slow down the code by about 25%. I will leave it at 1e11 cm and 16 buffer zones.
- If 5 is used at the surface instead of 8, this creates an artificial wave pulse traveling outward that leads to large pressure variations and catastrophically slows the code. This catastrophic slow-down is delayed by refining at a high level out to a larger radius in the ambient, but it is not clear whether this could avert the slow-down completely. Moreover, such a wave pulse is unphysical and affects the morphology at the surface. For Pamb=1.0e5 dyn/cm2, (as used in the original RGB fiducial run, Run 143) we have 8 level 4 cells per pressure scale height at the surface. For recent trial runs with Pamb=1.2e4 dyn/cm2, we have 5 level 4 cells per pressure scale height at the surface, which is inadequate. This result is independent of the EoS used, as the same problem happens when the gamma=5/3 ideal EoS is used.
- I managed to implement a new refinement scheme such that a shell covering the surface can be refined at a higher level than the rest of the envelope. In order to use Pamb=1.2e4 dyn/cm2 and fully resolve the surface (at least 8 cells per pressure scale height), I found that the shell needs to have a thickness of at least ~4 Rsun. If I use 5 Rsun, I find that the total number of maxlevel 5 cells in this shell is 4e8, which happens to be almost the same as the number of maxlevel 4 cells needed in the rest of the initial condition. I also experimented with gradient based refinement but could not get it to work very predictably or smoothly.
- I implemented tracers for
- core, star and ambient (3 tracers)
- Initial dominance of HI or HII and HeI, HeII or HeIII (5 tracers; testing still to be done). For example HII traces the gas which has an initial ratio HII/HI>1.
- I am working on some EoS notes that Yisheng will be using to create the remaining EoS tables (particularly the one that includes helium but not hydrogen recombination energy).
- Yisheng has completed the tables for the full MESA EoS with and without radiation energy.
Discussion
The numbers below correspond to those above.
- This problem is not all that serious and although it would be nice to improve it in the future, I think the stability is good enough to move forward.
- Better to resolve the surface even if it means settling for a higher ambient pressure. Therefore, I will go back to Pamb=1e5 dyn/cm2, as in the fiducial run 143. However, I will reduce the ambient density from 6.6e-9 g/cc to 1e-9 g/cc. This will cause a slow-down, probably by a factor of ~1.5, but that should be manageable. Also, now that aT4 has been removed from the EoS, the higher ambient temperature that this causes will not lead to a high ambient energy density, which could have posed a problem.
- Another option is to refine at level 5 in a thin shell containing the surface which would allow us to run with Pamb~1e4 dyn/cm2. The slow-down caused by the extra refinement would be compensated, to some extent, by the lower ambient temperature.
- However, there are limitations:
- The star becomes aspherical due to tidal forces, but the shell is spherical, so parts of the surface will become unresolved. This would happen during the plunge-in. Refining on pressure gradient would fix this in theory, but I couldn't get it to work very well in the tests I did. It tends to refine almost everything or nothing, and for some reason the tolerance needs to be set to ~0.8 instead of 8 to get the 8 cells per scale height desired, which I don't understand.
- The file sizes would be larger by a factor of ~2 since the number of cells would go up by that factor.
- The shell is never perfectly spherical or of precisely constant thickness, which might mean that it will need to be thicker than the 5 Rsun estimate, which would require even more resources.
- So I think it is best to keep things simple and refine the whole star at level 4 (out to 5e12 cm) as in past runs (with level 7 at the core) and settle for high pressure and temperature in the ambient medium.
- However, there are limitations:
- The tracers will be useful for distinguishing ambient gas from star gas, tracking the core gas to understand how stable the core is, and for understanding how the ionization ratios of specific gas elements change with time.
- The notes will be ready soon.
- The new tables (table 3 with no recombination or radiation energy and table 4 with no hydrogen recombination energy or radiation energy) will be ready soon.
Next steps
- Test tracers
- Implement new tables, having slurm script specify which EoS tables directory to use, and properly recording this info during the run.
- Start runs 0 (gamma=5/3 ideal gas) and 1 (MESA EoS)
- Finish EoS notes
- Once I am happy with how those runs are going, I will work on doing a single star run to explore convection
- Rename module to "CE" and upload to github?
- Start writing the paper
Long Run in 3D
The following run uses a larger radius than my previous 3D runs, and runs for a significantly longer amount of time
For those wishing to create large gifs like this (particularly if you want to do so on your own computer), I have a few tips to share:
- The 'convert' command is extremely resource intensive if you want to combine thousands of images. To save space and time, consider making smaller movies first and then combining (ignore the underscores, they are their to get the wiki to accept the following code)
_convert img/img00*.gif movie_00.gif _convert img/img01*.gif movie_01.gif ... _convert movie_*.gif movie.gif _rm movie_*.gif
- I have found other gif manipulators (eg gifsicle) perform better than convert
testing
testing
A Series of runs in one dimension
I have done a sequence of runs varying beta (the power of temperature dependence) for the one dimensional problem. These runs are set with the following parameters in order to achieve a cooling length of approximately 1 cm:
1E17 particles per cubic cm
Flow temperature is 720K
Mach 15
This gives us T_0 (shock temperature) of 5.06250E4 K. To achieve the desired cooling length at this temperature we use alpha = 8.27393E-23 ergs/s per cubic cm.
Beta varies between the runs.
Beta = -1.0
Beta = -0.5
Beta = 0.0
Beta = 0.5
Beta = 1.0
Beta = 1.5
Beta = 2.0
Beta = 2.5
Beta = 3.0
Update 09/14
Forming a qual committee
- I've got all five people confirmed to be on the committee. They are Adam, Eric, Segev, Miki, and Rip. I will check in with them for scheduling (also need to work with Laura and Jean). Probably a week in February/March.
CEJet
- all runs:
- Super jet runs:
- Normal jet runs:
- Small P2 runs:
- Plot and analyze the data this week.
- What happened with the small companion runs?
Update 9/5
Since my last post, I have done the following
I have done a sequence of runs varying beta (the power of temperature dependence) for the one dimensional problem. The values of beta (in order) are -1.0, -0.5, 0.0, 0.5, 1.0, 2.0. I also ran -2.0 but decided not to analyse it yet since -1.0 and -0.5 were similar.
I have done another reading of "Numerical Analysis of the Dynamic Stability of Radiative Shocks" (Strickland and Blondin 1995) I also did my first reading of "The structure of knots in variable stellar jets - I. Symmetric knots" (Falle and Raga 1993)
Lastly, I did a few more runs in the 3D Jet collision problem. First, I did a sequence of runs in which the Jets are (anti)parallel but not aligned, with the offset decreasing with each successive run. As I expected, a portion of each jet continues to propagate forward relatively unaltered by the collision. If the separation is not to great, another portion moves diagonally having been shaped by the collision; The angle at which this portion moves however did not depend as significantly as I had expected on the separation between the beams
A final run returns the jets to axial alignment, but makes one jet significantly hotter than the other