Progress:

For a while now I've been working on a new setup based on papers investigating dust destruction in so-called "windblown bubbles" around SNe. Most studies use several codes for this (one to simulate the immediate explosion aftermath to generate the shock structure, one to simulate the expansion of remnant for a few decades, one to simulate the expansion for several thousand years after & interaction with ISM). I'm trying to do the whole within one simulation. The most crucial part seems to be reproducing the correct shock speeds and structures. I think I'm pretty close now but will probably increase the initial velocity a bit to better match observations of Cas A. There are also some cosmetical issues (instabilities forming in the CSM where AMR levels change) plus issues at the boundaries. Will re-introduce dust once that's done.

Here's a gif of the current simulation going to almost 2k years: ​gif

And the individual outputs: ​1, ​2, ​3, ​4, ​5, ​6, ​7, ​8, ​9, ​10, ​11, ​12, ​13, ​14, ​15, ​16, ​17, ​18, ​19, ​20.

Two questions about this:

• are the instabilities developing along the y-axis caused by numerics?
• I'm getting string reflections from the centre which I'm not sure are physical, any way to suppress that?

Objectives

• Thesis stuff (main focus)
• HPC proposal (due this month so will need to focus on that starting next week)
• new problem setup runs

Progress:

Advice from last time solved the issue I was having with the ridiculously long run times. I had to redo the problem set up a bit to clean it up but it's now working fine. Looking at dust processing next. (​here is a screenshot of the current initial setup with a "better" random distribution of clumps)

Objectives

• Thesis stuff (main focus, still very busy with this, less progress on other stuff)
• Debugging & Testing

Progress:

Some updates on the new problem setup: I've added a random distribution of clumps to the inner ejecta, and the 1D density profile from ​last time now looks like ​this in 2D.

As mentioned last time, I have no information about the temperature or velocity profile so I've made some guesses there. The temperature is less problematic, the velocity is a bit difficult. The idea is to find an initial velocity that will give us a forward shock at 2.5 pc after 330 yrs (suggested by observations).

I Started a simulation to test an initial velocity estimate but I've run into a problem with simulations taking ages. The current estimate for the simulation to finish is 8 years… I feel like I might have made a mistake during the setup. Some settings:

• 32 initial cells with 8 AMR levels giving us the same resolution as the Slavin+2020 paper
• HLLC solver and parabolic interpolation though if I go to HLL and constant it doesn't change much.

And some outputs from the simulation:

• ​Hydro CFL issues
• Simulation restarts several times until ​no more CFL issue but I guess the timestep is just tiny now

Potential error sources:

• I removed the smoothing length h from the clumps because Slavin+2020 doesn't use that
• scaling???
• something else?

UPDATE: Advice helped, simulation is now running quickly! The expansion velocity is currently a bit fast so I will need to adjust that but it does form the ​forward shock close to the expected radius.

Objectives

• Thesis stuff
• Debugging & Testing

Progress:

• still mainly busy writing/correcting thesis chapters
• setting up new problem setup to reproduce ​Slavin+2020. Simulation is currently in the queue, no videos to show yet but I've plotted the initial density profile ​here. Velocity and temperature profile are still a work in progress, currently trying to find the best fits to reproduce observations.

Next Steps:

• writing, writing, writing
• potentially further work with the setup, have to discuss with supervisor

Objectives

• Debugging
• Production runs for some setups
• Thesis writing

Main progress over break:

Writing: I've mainly been writing over the break and will probably be quite heavily focused on that for the next few weeks as well

Time-stepping issue with the Gas Drag: I'm pretty sure the time-stepping issue for the gas drag is NOT due to the equations being too stiff (not yet certain about the grain-grain collisions, still working on other stuff there). It looks like one cell eventually has a temperature that messes stuff up and leads to a runaway increase of number density (see screenshots: ​1 ​2 ​3 ​4). Need to go back to the equation next week to see if there are any caveats for certain temperature regimes that might cause this.

Next steps:

• Continue with above

​Link to the code flow chart with suggestions re sub-cycling. (EDIT 1: pdf was locked, should work now)

EDIT 2: I have added a schematic outline of the code for the gas drag routine to showcase the dependence on the time step and the move between different FORs. It can be found ​here. Some more notes of the equations can be found in ​this blog post from May.

Very short post this week:

Good News:

I've compared my intermediate and end results per cell/bin/timestep with the external dust processing code and they match up. Meaning I get the same rates for dust advection/destruction!

Bad News:

It seems that for the Spitzer Plasma Drag and the Grain-Grain Collisions the rates with which the dust densities change from one timestep to the other are too large and AstroBEAR crashes with a "minimum timestep"-error…. not sure what to do here, would be great to discuss a bit.

Objectives

• Debugging

Progress:

Gas Drag: Did more benchmarking of intermediate calculation steps and they all look fine so still not sure what's going on. I also ran the simulations discussed last time:

• no wind, just a clump with initial velocity in x → looked fine, did not see the same effect as in the wind+drag simulations
• Fixed temperature to check the temperature dependence → that's noticeably different. Made a couple of videos: ​1. Normal temperature dependence, ​2. Fixed high temperature (around 2E8K), ​3. Fixed medium temperature (around 2E6K), ​4. Fixed low temperature (around 2E4K). To summarise: Movement looks fine for large T (though I suspect the drag force is still too strong), the movement looks very strange for small T. For small T is looks like the attraction between charged ions in the wind and charged dust grains is so large, that the grains are sucked towards the wind. Potential error sources here: Faulty gas ionisation degree (have checked this before but will again just in case), maybe T is off, drag force too strong, something else? (made a couple of screenshots for the low T case as well with more output steps: ​1 ​2 ​3 ​4 ​5 ​6 ​7

Grain-Grain Collisions: Made quite a bit of progress here, found a lot of smaller errors and the routine now works (previously it crashed because the collision probabilities exceeded 1) but I'm gaining mass here as well just as with the sputtering routine. → Will check the routine which sorts the grains into bins next, if mass conservation is broken, it's most likely because something goes wrong in there. UPDATE: Might have found the bug here, still needs testing but hoping to show some first simulations next week.

Next steps:

• Continue debugging

Very short update this week:

Objectives

• Debugging

Progress:

Gas Drag: I check the error sources I had identified and fixed some bugs (mainly typos with some constants) and both Spitzer1941 and Fry2020 grain charges now produce realistic charges. The simulations do look different but I'm very convinced now that there's something wrong still (​weird finger-like structures)

Sputtering and Grain-Grain Collisions: Got some benchmarking simulations set up for comparison to check whether the rates of destruction are reasonable or not.

Next steps:

• Continue debugging hopefully shifting focus to sputtering soon

Objectives

• Debugging

Progress:

I'm gonna try and do one dust physics process per week now so I'll be done by the end of October, this week it's the gas drag routine.

Gas Drag: I have found the following issues so far and am currently working to fix them:

• My grain charge is (as expected) wrong by quite a substantial amount. As grain charge affects all dust processes except collisional drag, this is a major bug that probably messes up a lot.
• ​One bin is a bit strange with a lot of dust seemingly unaffected by the drag until ​much later in the simulation. Not sure if this is a bug or not… need to check if this behaviour continues once the grain charge is fixed/when I change the distribution.

Sputtering: Some corrections/updates regarding my last blog post:

• It looks like the destruction rate, in general, is somewhat ok contradicting what I said last time. With my current initial grain size distribution, no dust is expected to survive beyond 20 years. Here are the updated plots with a new x axis: Total mass over time ​m(t) and change in mass over time ​dm/dt (note: the y-axis is in code units which doesn't fit the x-axis. Will fix this soon but for now you can clearly see different dust-gas interaction phases which are not easy to identify in the other plot).
• However, I have found that I am leaking dust mass somewhere so that might be an issue with my grain sorting routine.

Next steps:

• Gas Drag this week, then sputtering, etc.

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

This is a mid-week summary of the past two weeks as I won't be around much for the next couple of weeks:

Objectives

• Debugging the dust processing routines (Gas Drag, Sputtering, Grain-Grain Collisions)

Progress:

• Gas Drag: I made a comparison between the different gas drag methods (collisional drag, plasma drag, collisional + plasma drag) which is available ​here. Bin 1 is doing weird stuff, not sure why so much material is left behind, need to check what's going on there. I'm also quite certain that the plasma drag is a bit too large. This might be due to the grain charge which I think is a factor 10 larger than it should be. Still tracking down why that is (probably a messed up constant somewhere).
• Sputtering: I'm done with the first big debugging round of the sputtering routines. The bug that causes the crashes is fixed but I'm pretty sure I'm destroying too much dust at the moment (meaning my destruction rate is too high - probably again some faulty constant), and I'm still trying to fix the upper and lower collector bin: ​This is a simulation from earlier this week showing transport between the bins but weird stuff is going on in the upper and lower collector bin. Fixed the mass conservation and now I get ​this. Upper collector looks good but again I think destruction is progressing too quickly and obviously something is very wrong with the lower collector bin. Need to check whether that's a dust issue or a boundary/inflow issue. Also working on a plottings script to show the evolution of the total dust mass per bin across the entire domain to check whether the mass is indeed conserved or not.
• Gas Feedback: Gas feedback for sputtering is done but I didn't run a simulation of it as it heavily relies on the lower collector bin to work and since that isn't the case yet there's little point in feeding faulty number densities back into the gas. But the basic principle is this: It is assumed that the material in the lower collector bin is dust sputtered to such a small grain size that the material can be considered gaseous. I then calculate the total mass in the collector bin, translate that to a number density of gas atoms and feed that back into the hydro gas while setting the content of the lower collector bin to 0 (or rather the minimum dust density).

Up Next:

• Next three weeks: Probably won't get much coding done since I'm on break for a couple of days, then I need to finish the last stage of this year's HPC proposal and then I need to finish a thesis chapter. Might try to work on the collisions though if I get the chance.
• More debugging for gas drag & sputtering + getting started with debugging the grain-grain collision routines
• Gas feedback for gas drag + grain-grain collisions

Objectives

• Debug the processing routines (focus last week: plasma drag)
• Continue with the grain-grain collisions

Progress:

• Still testing stuff but plasma drag might be working now as well (see short description ​here)
• Here are some screenshots (couldn't make videos, something's wrong with my cluster atm): ​Collisional Drag only (from last week), ​Plasma Drag only (new), ​Plasma and Collisional Drag (new)

Up Next:

• Debug sputtering routines (probably this week and next week)
• Continue implementation of collisions (probably in Sept. after the 2nd round for the HPC proposal is done)

Objectives

• Debug the processing routines
• Continue with the grain-grain collisions

Progress:

I'm somewhat confident that the gas drag is now working correctly. There were a couple of smaller issues in the routines (e.g. wrong units) and I had accidentally counted the ion impacts twice instead of counting ions impacts and electron impacts once each. With the updates, it now looks like the drag might be ok or at least I can't see anything obviously wrong with it:

• Here is a video of what it used to look like without coupling the dust to the gas: ​No Dust-Gas Coupling
• Here is the video of what it now looks like with dust and gas coupled via (collisional) drag forces: ​Drag Dust-Gas Coupling (more videos: ​momentum, ​number density (original limits))

A couple of things to note here:

• no matter how long I ran the simulations or with which initial dust grain size distributions, it seems the bug where the dust number density just explodes at some point is gone, as is the NAN error
• the RK routines also don't throw any timestep related errors indicating that the equations used so far are not too stiff which is good!
• there's some ​pile-up happening in the centre cell row towards later timesteps; probably happening because there's not dust pressure as discussed previously…
• there are some ​numerical issues at the start in regions where the dust and gas have not yet interacted with each other significantly; probably also because of lack of dust pressure so no y-velocity yet in those dust fields…

Up Next:

• Do the same sort of study with the plasma drag
• Debug sputtering routines
• Implement collisions (see theory notes from ​last week)

Objectives

• Research
• Code development (debugging + grain-grain collisions)

Progress:

Spent last week on theory research again, summaries are below

• Distributions: See pdf for the power-law distribution ​here. (Currently working on a couple more here + some user-defined options to make the code more flexible, will add once that's done.)
• Grain-Grain Collisions: See pdf ​here

Up Next:

• Back to debugging and implementing the things gathered over the past two weeks.

Objectives

• Continue testing & debugging
• Work on Grain-Grain Collisions

Progress:

• HPC Proposal: This is finished for now, the next phase will be in September.
• Debugging: Ongoing benchmarking using external processing.
• Grain-Grain Collisions: See the pdf ​here!

Up Next:

• Need to look some more at the fragmentation case
• Continue with benchmarking

(I'm on leave next week so I won't be at this meeting!)

Objectives

• Debugging (now moved on to sputtering as well)
• Refine dust routines

Progress:

• HPC Proposal: I need to write one this month so progress on other topics might be slow towards the end of the month.
• MHD Clouds Crushing Simulations: Did a ton of simulation for our next paper to show the effect of different magnet field orientations on the Cloud Crushing Setup. I will make another blog post with plots and videos at some point but my computer decided to transition to laptop heaven and while I'm waiting on a new one I can't properly use visit.
• Debugging: Made progress here and have also started looking at the sputtering but I don't trust it yet and want to invest more time comparing my results with the dust post-processing results. Currently, simulations look a bit "grainy" which seems off considering that the gas looks smooth.

Up Next:

• Haven't forgotten about the ​equations for the dust grain size distributions! I will add them once I've checked something and made sure I'm actually implementing what I want. Looks a bit strange at the moment.
• Continue debugging of sputtering + check ion trapping & gas feedback is implemented correctly.
• Setting up Sputtering comparison with the dust post-processing code.
• I'm aiming to do production runs for drag and sputtering in August and while that's going on I'll add the rest of the grain-grain collisions code.

Objectives

• Debugging
• Refine dust routines

Progress:

• HPC Proposal: Still really busy with this, need to get through 1m core hours by the end of June…
• Dust Grain Size Distribution: Previously the code was designed to read in a dust.data file which contained the dust grain size distribution information. I realised that might not be so great cause it would be hard to reproduce a simulation if that file were to be lost + it should scale with the volume of the cell. After some edits, the dust grain size distribution is now set via input parameters in physics.data and then calculated during setup initialisation.

(I'll try and add the equations for the different distribution types when I get a chance; a lot of unexpected meetings today so I might not manage before the meeting tonight.)

Up Next:

• Continue Debugging

Objectives

• Debugging
• Refine dust routines

Progress:

• HPC Proposal Had to do some work there for the current & next round. We want to do different setups next year to look at different regions of dust destruction
• Simulations We're heading towards the end of Q1 so I'm currently really busy setting up some final simulations (post-processing)
• Debugging Made some progress regarding gas drag. Now the gas in the ambient medium is moving (see ​screenshot) but not in the clump. Not sure whether that is still density related or not, will need to check. I'm also adding a dust distribution routine so that those can now be calculated inside the code instead of with an external Python script beforehand (makes it more compatible with the post-processing).

Up Next:

• More debugging

Objectives

• Debugging
• Refine dust routines

Progress:

(A lot of research clusters in the UK have been down for almost three weeks now due to some security issue. I've moved my stuff to a working cluster but queue times have been long so progress has been slow.)

Dust Drag Progress & Issues

• Comparison with Dust Post-Processing: I've got the necessary simulations to compare my routines with the results from the post-processing
• Debugging: I fixed some stuff this morning but it must have broken something because suddenly I'm somehow getting an inflow from the lower x boundary about halfway through the simulation with passive advection: ​Start of Simulation → ​Middle of Simulation. For active advection, the velocity is suddenly completely gone ​Gas and Dust Momentum (FIXED! Forgot to remove some debugging stuff from last week and that messed stuff up there)

Update

• Debugging: The inflow was a bug in the dust wind initialisation. Some time ago I modified the wind so that the shock could come with or without dust. Turns out that causes an issue if the conditions inside the domain allow the bow shock to hit the lower boundary. I went back and changed that case to allow dust inflow and I found that for dust density 1e-16 the dust in the wind has no noticeable effect on the dust in the domain so effectively that makes the wind dustfree. ​Dust Clump evolution with negligible dust in the wind and passive advection
• In other news: My velocity values are now only 1e4 off down from 1e8… getting there

Up Next:

• More debugging

Objectives

• Debugging

Progress:

Dust Source Routines

Not much to report this week. Have been busy debugging and progress is ok but probably still a lot to do for each processing & drag type. Here's a comparison between ​last week and ​this week. Dust Velocities are better now. Need to set up a simulation with realistic values now as so far I've just been using random numbers which in the case shown here resulted in very heavy dust grains.

Last week, I noticed some ​weird behaviour of the dust x-momentum in the first couple of steps but I wasn't able to reproduce that. It now looks ​fine so that might have been an issue with my local cluster.

Up Next:

• More debugging

Objectives

• Dust Source Terms for Dust Advection

Progress:

Issue: Code crashes with dust source routines

Previously, the dust source terms routines would always lead to the code crashing. I have no integrated the dust source terms with the new Riemann solver for the dust and the code no longer crashes. A simulation can be found ​here and here are two screencaps: ​Fig 1 and ​Fig 2. There are a few issues:

• The dust momentum is incredibly small (~ 1e-50) which is not realistic. I suspect units/conversion issues which is not surprising and that's been on my debugging list for a while. Because the dust momentum is so small, the dust grains are not accelerated enough to noticeably move which is why the dust clump remains at rest.
• The non-zero dust momentum structure in Fig 1 is strange and I'm not sure yet what's happening there. Might be a setup issue…

Issue: RK2/45 and the dust source terms

The drag calculations we use calculate the total change of the velocity of the dust over one timestep meaning dv (summary is here: ​PDF). However, the RK routines don't expect the new velocity but the change in velocity over time. For a simple integration scheme, I can simply declare dvdt = dv / dt where dt is the current timestep. This works if dt is constant throughout the integration but for RK45 for example, dt is not constant and the integration also features the RK h factors. I could rewrite the RK45 routine to adapt but that would be messy. I've been looking through different source routine now to see how other routines handle it, in particular, PointGravity and I'm a bit confused about the following:

• Why does moving PointGravity from SrcDerivs to SourceCell improve the momentum conservation?
• Does PointGravity simply change the q-vector before the integration? If so, how does that not become unstable?
• From what I understand, the issue would be solved if I move my routines into SourceCell and instead of altering dqdt, I should be altering q directly?

Up Next:

• Continue Debugging

Objectives

• Grain-grain collisions routine (draft)
• Source term integration (debugging)

Progress:

• Sod Shock Tube: As suggested, I ran the dusty Sod Shock Tube problem setup with different interpolations methods: ​PCM, ​PLM, and ​PPM.
• Dust Processing: Made progress regarding grain-grain collisions & source term integration but not enough to run simulations right now

Up Next:

• Continue source term integration & draft grain-grain collisions

Objectives

• Properly integrate dust source routines
• Continue with grain-grain collisions
• MHD production runs for post-processing dust treatment

Progress:

• Sod Shock Tube: Finished (for now). The issue before was that the gas movement is triggered via initial values of the energy and as the dust has now energy component it remained at rest. Thus active advection + Sod Shock Tube is currently not possible (it would require writing a different initialisation routine and I've noted it down as something to be done in the future but I don't feel it's necessary right now), passive advection + Sod Shock Tube, on the other hand, works: ​Setup | ​Gas Density (​Video)| ​Gas Velocity in x (​Video) | ​Dust Density (​Video)
• Dust Post-Processing: We were awarded HPC time so I'm currently setting up MHD (instead of just HD as in the previous papers) runs for the dust processing done by a colleague. (Not sure about the best way to visualise field lines…)

Up Next:

• Source terms!! plus further development of dust routines but source terms now take priority

Objectives

• Set up Sod Shock Tube Problem for dust
• Properly integrate dust source routines
• Continue with grain-grain collisions

Progress:

• Grain-Grain Collisions: Made some progress here while I didn't have internet access would like to pause it now though to focus on the source routines from next week onwards.
• Sod Shock Tube: The default Sod Shock Tube problem is incompatible with the dust modules so I'm trying to integrate the dust by setting up my own Sod Shock Tube using split region. It did run but due to lack of internet, I couldn't check whether it was doing the right thing and unfortunately, the starting setup looked like ​this. I will try and rectify that next and maybe reduce the problem to 1D as I don't really need 2D here.

Update:

• Sod Shock Tube: Turns out the issue above was a) due to using 2D in global.data when problem.f90 was for 1D and b) colourmap choice. Using 1D and a different colour map I get ​this for the gas which looks alright (density cut ​here). However: While the dust is initialised it's not given any momentum at the beginning so it just remains stationary. Need to find out why the dust velocity initialisation is not triggered.

Up Next:

• Figure out dust velocity issue for Shock Tube problem
• I really want to get the source routine integration working next so that will be my priority over the next weeks.

Comment:

I wasn't able to make blog posts for the last couple of weeks so this is a summary of the progress made since then.

Objectives:

• Finish active advection routines
• Extend interpolation options to dust
• Continue dust processing routines

Progress:

Interpolation

• I've extended the interpolation scheme to include dust. So PCM, PLM, and PPM are now also available for the dust.

Dust Advection

• I reduced the Cloud Crushing setup to just a clump moving through gas as I wasn't sure whether there's an issue with the boundaries, the wind, the fluxes, etc.
• Found out that way that the x- and y-fluxes kept getting messed up → fixed in F_prim and included and added an HLL solver for the dust to both the HLL and HLLC gas solver.

Below are videos of the tests I ran to check the advection:

Simulations with the HLL Gas Solver (and HLL Dust Solver)

• Moving Clump (initial velocity in x): ​https://youtu.be/RcKyjCDi2HM
• Moving Clump (initial velocity in y): ​https://youtu.be/eQp_NtMpx0M
• Moving Clump (initial velocity in x and y): ​https://youtu.be/Z0rndg-28ug
• Moving Clump (initial velocity in x and y with the dust velocity equal to 1.5 times the gas velocity): ​https://youtu.be/_rNcDg2Ba64

Simulations with the HLLC Gas Solver (and HLL Dust Solver)

• Cloud Crushing Setup (constant dust distribution): ​https://youtu.be/uE0pXMotopQ
• Cloud Crushing with HLLC and PCM: ​https://youtu.be/2W6eQVvJEgY
• Cloud Crushing with HLLC and PLM: ​https://youtu.be/D48wUrkqpwg
• Cloud Crushing with HLLC and PPM: ​https://youtu.be/rhRvjBOLxjQ
• Cloud Crushing with HLLC and PCM and reflective boundary along the y-axis: ​https://youtu.be/Zi-2GXu0oQc
• Cloud Crushing with HLLC and PCM and dust grain size distribution across the bins: ​https://youtu.be/0O9YaNauS1o
• Cloud Crushing with HLLC and PCM in 3D (Slice): ​https://youtu.be/GJWClB9gUEA
• Cloud Crushing with HLLC and PCM in 3D (Contour): ​https://youtu.be/PbXodsg2PZ4

This looks like it's working alright considering the dust is implemented as a pressure-less fluid. Once the dust and the gas are coupled via drag, the dust should start getting a velocity in y as well.

Next steps:

• Now working on turning on the drag routines I have in the code (they still need to be properly integrated into the source routine)

Objectives

• Debug & test advection & sputtering routines
• Continue implementation of grain-grain collisions

Progress:

• Debugging: Still debugging, found one big bug and fix it last week so going well. Found an issue with the interpolation routines for the dust. Not quite sure yet what's happening here.
• Grain Growth: I think the ion trapping routine is more or less bug-free now, working on feedback effects on the gas since we couldn't consider that for dust post-processing.
• Grain-Grain Collisions: I'd say this is about 30% done now, aiming to finish by the new year
• New Problem Setup: Finished tweaking that. A sample simulation can be found here: ​Video.

Issues & Questions:

• Performance: Is there a way to find out how much time the code spends in certain routines?

Objectives

• Debug & test advection & sputtering routines
• Continue implementation of grain-grain collisions

Progress:

• Debugging: Currently cleaning up the routines. Fixed some stuff but still ironing out some issues.
• Grain-Grain Collisions: Continued the collision routines a bit but as mentioned this will take some time as I want to prioritise sputtering and drag.
• New Problem Setup: Set up a problem for supernova ejecta interaction with a wind-driven shell.​Illustration here. (Data for time-dependent wind can be taken from ​Schaller et al. 1992 and ​Vink et al. 2001)

Issues:

• Dust Processing: Still debugging but going ok.
• New Problem Setup: Working but is there an easier way?

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)

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.

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

Current Objective

• Debugging the dust drag advection routines (collisional drag)
• Finishing the dust sputtering routine
• Finishing the total mass calculation routine

Current Status

Drag Advection

(See also report from last week ​https://www.pas.rochester.edu/astrobear/blog/fschmidt10142019)

• dqdt for the dust velocity is way too large (in the order of 1e42) which (apart from being quite obviously incorrect) causes an issue with the RK integrators
• I found a couple of bugs already (e.g. there was an issue with the transformation into the dust grain rest frame) but it's still giving me values that are too big. The frame of reference is changed multiple times throughout the calculation so I suspect more errors in there somewhere. Currently still debugging.

Sputtering Routine

• Started implementing the sputtering routines based on Dwek & Arendt (1992), Tielens et al. (1994), and Nozawa et al. (2006). This is a work in progress.

Total Mass Calculation

• Did some more work on this internal processing routine which calculates the total mass of dust in each bin and across all the bins for each timestep (not just each output timestep) and stores them in a file. This additional dust-related output can be turned on or off and will mainly be used for 3D simulations (in case we don't get a new HPC award next year; our local facilities can't handle many AMR 3D Chombo output files. So using this output option would allow us to cut down on the Chombo output while still getting dust mass information with high temporal resolution). Still a work in progress but should be done soon-ish.

Questions

Feedback on Gas

Accelerating the dust should remove energy from the system and growing or destroying grains should either generate or remove gas mass from the system. I've been thinking about how best to implement that:

• For the grain destruction or grain growth, the most obvious choice would be to simply add or remove material to/from the gas during the source step. So

\sum_{i=1}^{nDustBins+2} \left( \bigtriangleup m_{\textrm{dust evaporation}, i} - \bigtriangleup m_{\textrm{dust growth}, i} \right)= \bigtriangleup m_{\textrm{gas}}

• For the gas drag the defining quantity that's calculated is the force F_{\textrm{dust}} the dust grains experience. F_{\textrm{dust}} is then used to calculate the acceleration that each grain experiences. My current approach is to set F_{\textrm{gas}} = - F_{\textrm{dust}}, where F_{\textrm{gas}} is the total force experienced by the gas particles in one cell. The force which a single gas particles experiences is then F_{\textrm{gas}}/n and that then results in a declaration. All of that is currently calculated in the source term routine alongside the dust acceleration. Questions: Will this result in energy conservation problems?

Stiffness of Equations

I'm very certain that the current issue with the RK integrator error is down to a bug and not down to an intrinsic issue with the stiffness of the drag equations (since by definition the acceleration of the dust grains cannot exceed the acceleration of the gas particles). However, it is, in theory, possible that the destruction routines are very stiff under certain conditions (very high densities, very high velocities etc.). Question: I was wondering what could be done in that case?

(Can't make it to today's meeting due to a time conflict with a seminar I am required to attend. Back next week!)

Current Objective

• Developing the dust drag advection routines (collisional drag)

Current Status

• Collisional drag subroutine is implemented though simulations crash with the following error message: "Src Error::RKOrder45 — timestep insignificant. Halting." I'm currently checking whether this is a bug in the source routine I've set up resulting in an unrealistically rapid dust acceleration or whether the system is too stiff.

Supporting Material

Gas Drag Equations (for Collisional Drag)

The new dust velocity after one time-setp is given by:

\vec{v}_{\textrm{dust}}(t + \bigtriangleup t) = \vec{v}_{\textrm{dust}} (t) + \bigtriangleup \vec{v}_{\textrm{dust}}(t) \, ,

where

\bigtriangleup \vec{v}_{\textrm{dust}} (t) = \frac{\vec{F}_{\textrm{drag}} (t)}{m} \cdot \bigtriangleup t\, .

The drag force \vec{F}_{\textrm{drag}} is given by

\vec{F}_{\textrm{drag}} = 2 \sqrt{\pi} k_B T_{\textrm{gas}} a^2 \sum_j n_{\textrm{gas,j}} \mathcal{F}_{\textrm{col,j}}\,.

The sum over j covers all particles species in the gas (we are using an average value for an oxygen dominated SNR). \mathcal{F}_{\textrm{col,j}} is given by:

\mathcal{F}_{\textrm{col,j}} = \left( S_j + \frac{1}{2 \cdot S_j} \right) \cdot e^{-S_j^2} + \sqrt{\pi} \left(S_j^2 + 1 - \frac{1}{4 \cdot S_j^2} \right) \cdot erf{S_j} \, ,

where S_j is

S_j = \sqrt{\frac{m_{\textrm{gas,j}}}{2 \cdot k_B \cdot T_{\textrm{gas}}}} \cdot v_{\textrm{rel}} .

The remaining variables are (all in cgs):

• m: mass of a dust grain
• n_j: number density of gas species j
• T_{\textrm{gas}}: gas temperatire
• a: grain radius
• m_{\textrm{gas,j}}: mass of a gas particle of species j
• v_{\textrm{rel}}: absolute value of the relative velocity between gas and dust
• \bigtriangleup t: timestep

The resulting drag force is then multiplied with a scaled unit vector to extract the drag forces in each dimension.

(Refernces: Kirchschlager et al. 2019, Baines et al. 1965, Draine & Salpeter 1979)

Passive Advection of Dust

Last Status

• passive advection was working if the gas-to-dust density ratio was 1. For any other ratios, the flux calculations were way off resulting in a faulty advection and eventually in a termination of the simulation run.

Current Status

• passive advection is now working for all gas-to-dust density ratios and bin numbers. The bug was in the conserved ↔ primitive conversion routines of EOS.f90. For passive advection the conversion routine should divide by or multiply with the gas density not the dust density.

Test Simulation: ​https://youtu.be/KJW2fLnDyG0

Next Steps

• Processing subroutine to store total dust mass in domain for each bin → Put in Sweep After Step in sweep_scheme.f90?
• Add other advection schemes (gas drag & plasma drag, starting with gas drag) → Put call to advection routines in Sweep Before Step in sweep_scheme.f90 but the actual calculations will be in dust.f90

Equqations

• Conservation of Gas Momentum:

\frac{\partial}{\partial t} (\rho_g \cdot u_g) + \frac{\partial}{\partial x} (\rho_g \cdot u_g^2 + P) = S_{\textrm{GD}_g}

where

S_{\textrm{GD}_g} = - \rho_g \cdot u_g \cdot \kappa_g + \rho_d \cdot u_d \cdot \kappa_d

• Conservation of dust momentum:

\frac{\partial}{\partial t} (\rho_d \cdot u_d) + \frac{\partial}{\partial x} (\rho_d \cdot u_d^2) = S_{\textrm{GD}_d}

where

S_{\textrm{GD}_d} = \rho_g \cdot u_g \cdot \kappa_g - \rho_d \cdot u_d \cdot \kappa_d = - S_{\textrm{GD}_g}