Published Papers

• The external dust-processing code PAPERBOATS which my internal processing routines are based on is introduced and described in ​Kirchschlager et al. 2019

Processing

The following publications are used as references for the various dust processing routines:

Drag

• The drag consists of two terms: Collisional drag and plasma drag. Descriptions are available in ​Baines et al. 1965 and ​Draine & Salpeter 1979

Oxygen and Dust Ionisation

• for the gas atom ionisation see figure 1 in ​Kirchschlager et al. 2019 which shows the charge number of the atom as a function of the gas temperature calculated using ​CHIANTI
• for the dust ionisation degree, we use the formalism described in ​Fry et al. 2018

Sputtering

• Collisions between gas atoms and dust grains can result in the erosion of the grains (if collisional energy is above the threshold energy) or grain growth (if collisional energy is below the threshold energy). If gas and grains are moving with the same velocity, collisions occur due to the temperature-dependent Brownian motion of the gas in which case we use a Maxwell distribution for the relative velocity. If gas and dust move with different velocities we use a skewed Maxwellian instead.
• The following papers were used as resources: ​Dwek & Arendt 1992, ​Tielens et al. 1994, ​Nozawa et al. 2006, ​Micelotta et al. 2016, ​Barlow 1978, ​Draine & Salpeter 1979, ​Bocchio et al. 2014, ​Shull 1978, ​Nozawa et al. 2006, ​Serra Diaz-Cano & Jones 2008

Collisions

• Collisions between grains can i) vaporise the grains ii) break the grains apart iii) grow the grains if they stick together
• The equations are derived in section 4.5 of ​Kirchschlager et al. 2019
• The following papers were used as resources: ​Bocchio et al. 2016, ​Borkowski & Dwek 1995, ​Hirashita & Yan 2009, ​Tielens et al. 1994, ​Jones et al. 1996, ​Dominik & Tielens 1997, [Blum et al. 2000], ​Yan et al. 2004

Ion Trapping

• Derivation of equations are in the paper that was just submitted, will add once it's on the arxiv!
• Penetration depth is calculated based on ​Bethe 1930 & ​Bloch 1933 (articles are in German)

sigma=B*B/rho*v*v

	sigma=B*B/rho*v*v=0
	sigma=6e-7

Radiation Pressure

Still waiting for responses to comments. Should probably send around again.

MHD

The clump test is also clearly not working. Additionally, the advection test is behaving very differently on the development branch than on the line_transfer branch. All of this suggests that there's a bug in one or more places.

Also, here's the rotating run test:

Charge Exchange

Here's frame 250.

Looks very similar to the last frame of the movie I showed previously. Once I get a movie up to the latest frame, I'll see if we can stop.

Scaling Test

Haven't gotten the KNL nodes working yet.

Objectives

• Fix advection routines
• Test drag & sputtering routines
• Implement grain-grain collisions
• Set up new problem

Progress:

• Drag Source routines finished: Dust is now advected using collisional drag (energy transfer when gas atoms collide with dust grains) and plasma drag (Coulomb force interactions between charged gas atoms and dust grains)
• Dust Sputtering routines finished: Dust sputtering routine is done and includes dust destruction due to collisions between grains and gas atoms (either due to relative motion between gas and dust or Brownian motion if both gas and dust are at rest) and dust growth if the collisional energy is below a certain threshold. Two additional feature still need checking but are basically finished as well: Trapping of ions inside the grain (if collisional energy is high enough so that the gas atom penetrates a certain number of atomic layers in the dust grain) and gas feedback (if gas is accreted or if the dust is destroyed material should be either removed or added from the gas density field).
• Started work on grain-grain Collisions but that will take a lot more work and I'd like to prioritize advection again
• Currently setting up a new problem to replicate simulations from a recent publication

Issues:

• Not really an issue but drag and sputtering routines need excessive testing now
• Advection: Didn't do much work on this as I was focusing on sputtering but I think I've now tracked it down to an issue with wrongly implemented dust velocities which then result in enormous wl and wr values.

New Work

• Advanced Run 183
• Estimated remaining cost of simulation assuming evolve to t=700 days (same as low res run 164) and assuming no changes to refinement (Answer: 30% of the total allocation or about half of the remaining node-hours, and about 40 TB of additional storage in addition to 30 TB currently)
• Modified slurm script to exit loop if error occurs (to avoid wasting node-hours)

Results

• Energy gain in 183 is again too fast! Now crossed 5% threshold and energy gain is accelerating.

Figures

Interpretation

• The most likely reason for the energy gain is (again) lack of resolution around the point particles
• In particular, I suspect that as particle 2 accretes more and more mass, the pressure gradients become too high to resolve
• Other, less likely (non-mutually exclusive) possibilities are:
  • insufficient resolution around particle 1 (but we know it was adequate at earlier times, so less likely)
  • insufficient resolution in the ambient gas, where Rayleigh Taylor instabilities are occurring (but this material has much lower density, and development of instabilities does not seem to correlate very well with increase in energy, so less likely)

Strategy

• I plan to add two AMR levels of resolution around particle 2 and see if this helps.
• If yes, then I need to perform tests to determine a more optimal refinement strategy.
• If we find that the refinement needs make the simulation too expensive to perform, then it will be time to write up the results for a paper.

Next steps

• More analysis to better understand what is causing the energy gain (extend 1D movies of pressure profiles around particles, and make movies to understand better the development of instabilities in the ambient and to what extent this correlates with energy gain)
• Test with higher resolution around particle 2 (restarting from an earlier time in the simulation)
• If the test shows improvement, then need to perform additional tests to optimize refinement strategy

MHD

Here's a side by side comparison of the same frame in the MHD and non-MHD runs. This is with a stellar beta of 0.01, which gives a beta of a bit less than 1 near the planet. Note that the field is clearly pretty stiff.

Do we want to see how high a magnetic field we require in order to suppress the up-orbit arm? Can calculate the R_M compared to R_t and R_w, to see if the approximate field value is even feasible.

Charge Exchange

Not quite steady state, but does look like it's almost there. Here's a movie up to frame 22:

Scaling on Stampede

Think I have it working now. Should be done in a few days.

Objectives

• Finish gas drag routines
• Fix advection routine
• Implement sputtering

Progress:

• Fixed the source term issue where dqdt was too large per time step resulting in a RK45 error (see update from the 4th of November). dqdt is now in an acceptable range but fluxes for the dust momenta go towards infinity → suspect there's an issue with the dust advection (this advection type hadn't been tested before)
• Finished implementing a second gas drag option: plasma drag (two options here: traditional plasma drag or a fancier model based on Fry et al. 2019). The dust is now accelerated due to a) collisions with the protons/electrons in the gas (collisional drag) and b) Coulomb interactions with charges gas particles if the grain itself is charged (plasma drag)
• Started work on sputtering routine

(Note: Since advection isn't working properly I can't properly test the gas drag and sputtering routine. Some debugging might be necessary once I can test them.)

Issues:

• The wl and wr values go to infinity when the fluxes are calculated for the dust momenta. I'm still tracking down what exactly goes wrong (potentially even a setup issue with faulty initial dust momenta?)

Objectives

• Fix problems with collisional drag routine (dqdt too large)
• Set up setup for SN shock CSM/ISM intreactions

Progress:

• Drag Advection Bug: The way I was using the routine implied that each collision between a dust grain and a gas particle is taking place with a proton. However, the gas is ionised so many collisions should be with electrons. Currently implementing an ionisation subroutine for oxygen and altering the drag routine to take this into account
• Have set up a smooth and clump version of the SN shock CSM interaction model (aiming to simulate 1987A)

Videos:

• Smooth CSM (¼): ​https://www.youtube.com/watch?v=qaXIDxjR5so&feature=youtu.be
• Smooth CSM: ​https://www.youtube.com/watch?v=Itf2gPRpeLg&feature=youtu.be
• Clumpy CSM (¼): ​https://www.youtube.com/watch?v=job2iSM_Brc&feature=youtu.be
• Clumpy CSM: ​https://www.youtube.com/watch?v=CY90pZh8psk&feature=youtu.be

Issues:

• Numerical artefacts at boundaries → Discussion outcome: unavoidable though higher resolution might improve it

New Work

• Advanced Run 183 and produced figures and movies
• Progress on deriving expressions for total energy conservation prescription
• Plotted convergence test 185 to confirm that changing refinement region center from particle 1 to particle CM has negligible effect on energy conservation (confirmed)
• Encountered (and solved) particle creation issue that was crashing code (code is now running on stampede2)

Results

• Rate of energy gain for Run 183 is now quite small and apparently on track for staying below 10% of the old AGB run (Run 164).
• Run 183 crashed because of particle creation….I created a switch that turns off particle creation.

Figures

Movies

​Run 183 Density
​Run 183 Density with mesh
​Run 183 Pressure profiles (sliced along x-axis in simulation frame of reference)

Next steps

• Continue Run 183 (and increase speed by using more nodes)
• Further modify new energy conserving routines to take into account particle-gas and particle-particle energy terms
• Further automate and improve simulation:
  • predefine refinement radius as a function of time (should vary smoothly rather than in sudden steps)
  • center of refinement shape should transition smoothly in time from particle 1 to particle CM (instead of suddenly)
  • experiment with using true AMR, e.g. refine on density or pressure gradients within refinement shape

Energy Balance Model (EBM)

The EBM reaches equilibrium at Teq. Below are plots of this for varying distances and partial CO2 pressures (pCO2):
​3D Plot (top view)
​3D Plot (side view)
​Contour Plot

Numerical Model

​Code Location on GitHub (Initial Conditions: t=1880, pCO2=284ppm, N=1.13 Billion People)

General Simulation Steps

1. Run program with Earth like inputs as they were before the Anthropocene (in 1880), until an equilibrium temperature is reached.
   
2. Set the populations peak growth rate to this temperature ​as shown here.
   
3. Then turn on population growth; exponential rise in population results in exponentially increasing partial CO2 pressures, ultimately triggering a global climate change. Once the global temperatures pass a certain point, the growth rate of the population becomes negative, which results in mass extinction.

Goal #1

Repeat steps 1-3 for four different values of initial partial CO2 pressure, ranging from a tenth of ours to one hundred times ours. At each new value of pCO2, determine the habitable zone for that point, defined here to be where water can be liquid (273-373K).
Then repeat steps 1-3 for 4 different values of orbital distance, spanning the entirety of that habitable zone.
​Results Shown Here (as 4x4 grid)

Goal #2

Choose 3 different orbital distances. For each one, increase pCO2 until the equilibrium temperature reaches 287.09K (the temperature of earth at 1AU).
​Results Shown Here (Earth included for reference)

Next Steps

1. Model 1: Introduce constant carrying capacity (maximum population)
2. Model 2: Introduce a temperature-dependence to the carrying capacity, to emulate reduction of habitable regions due to global climate change.
3. Model 3: Introduce a time-dependence to the carrying capacity, to emulate environmental sensitivity to climate change. So environments with low sensitivity could exceed their carrying capacities for a longer time than ones with higher sensitivity before climate change begins.

Analytical Model

​Current state of the Equations and Nondimensionalization Process
​Phase Diagrams of Analytic Model

Discs From Tidally Disrupted Companions:

Animation:

Tidal Disruption Simulations:

Goals →

1.) Investigate Bonnor-Ebert Sphere stability in a hot wind-tunnel

2.) Will a disk form from the cold tidally disrupted material?

3.) What is the fraction of material that forms a disk?

Just to mention: I wrote my own RK4 code in python that can modify the AGB profiles from MESA. Iterates over point mass, alpha, and central density, and integrates the pressure.

Animation:

R sculptoris:

I am working with Zhuo to simulate the system R Sculptoris, This will be similar to the L2 Puppis Paper where we attempt to predict the binary parameters and generate synthetic images with RADMC-3D. As of now Zhuohas prepaired a module to run and I am sorting out how the refinement should work.

Animation: