​Plots for the new jet runs

Charge Exchange

Here's a movie of the frames so far. Few more frames still until it hopefully speeds up and we can throw more cores at it.

MHD

Tried a run with the planet with outflow-only hydro boundaries and extrapolating B field. No dice. So many pressure protections that I had to delete the log.

Think this is confirmation that we need to change tack to getting spherical line transfer up and running so we can use global simulations (or one of the other global simulation methods we've discussed). This will require some discussion.

Job search

Making some good progress. Created a short resume - may need a little bit more modification, depending on the job. Need to get around to writing a cover letter for at least one application, so I have somewhere to start.

Also got a template for thesis from Erica. How does this outline​ look?

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

Coupled EBM

​Trajectories

Fermi

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)

Charge Exchange

Shock still working its way across the grid. 5(?) frames until we should hopefully get some speedup.

MHD tests

Extrapolating runs work well. Here's a comparison of the strong and weak cases of the corotating with stellar gravity runs:

Strong

Weak

Also tried adding magnetic field to outflow-only boundaries, which fixes the initial pressure protections, but kills the run later on (only get 3 frames, tons of protections after 1.5).

EoS stuff

Work done

1. Explored effect on recombination energy of including temperature dependence of partition function.
2. Made pdf with all the initial profile plots, for easy reference in the future.
3. Went back and did a bit more reading on the EoS, particularly the OPAL EoS which is used everywhere in the initial RGB profile ​Rogers+ 1996, ​Rogers+Nayfonov 2002.
4. Some reading, including this preprint, which is relevant for our EoS CE project: ​Hirai+ 2020.
5. Figured out that we were missing the recombination energy from recombination of two hydrogen atoms to produce a hydrogen molecule, which is included in the MESA EoS and discussed in Hirai+ 2020. So included that contribution in the plots.

Results

1. The temperature dependence of the partition function cannot explain the decrement in the internal energy near the surface (see last blog post). Therefore it must be caused by other physics included in the EoS but not captured by the contributions included (thermal+thermal radiation+recombination) or, another possibility is that the Saha equation does not accurately predict the recombination energy near the surface, perhaps because local thermodynamic equilibrium cannot be assumed. In any case, the profile gets modified near the surface (r >~ 45 Rsun up to the outer radius of 48.1 Rsun) because of the matching onto the ambient, so we are not too concerned about this region. However, during the simulation there will also be regions where it doesn't quite work…is this a problem for what we want to do?—update: I now realize that the decrement is caused by the neglect of the latent recombination energy due to recombination of H atoms into molecular H2! See the last item of Results and discussion below…
2. Here is the pdf with all the various profile plots for easy reference: ​eos_figs.pdf.
3. The Hirai paper mentioned above also has the same idea as me of including the recombination energy of H or not, of He or not, etc. in the simulation…but they do not use a tabular EoS. Rather, they use an analytic solution to approximate the tabular EoS.
5. Although the number density of H2 molecules is completely negligible, this component of the recombination energy H+H —> H2 is not negligible, and in fact dominates at temperatures below 10⁴ K (see Fig 17 in the above pdf). When we include all the contributions, excluding recombination energy of metals, the results for the internal energy density agree to within 0.6% everywhere. If we include recombination energy of metals as well, we do a bit worse, especially at radius below 12 Rsun, where the percent difference becomes as high as 1.4% (over predicts the internal energy by up to 1.4%).

Discussion

1. We have confirmed that by adding the thermal energy contributions of gas and of radiation and the recombination energy from hydrogen and helium, one gets back the MESA internal energy profile (the discrepancy is on the order of half a percent at most).
2. We can now say that it is reasonable to estimate the recombination energy contribution using the sum of hydrogen and helium recombination energies computed from the Saha equation.
3. We need not worry about the temperature dependence of the partition functions, as ignoring it produces results in good agreement with MESA.
4. Recombination energy of metals can safely be neglected, it seems, but recombination energy of H+H —> H2 cannot be neglected, though it was not considered in ​Reichardt+ 2020.
5. As mentioned in the last blog post, one can see from the above pdf that the latest version of MESA (12778) gives results that are completely consistent with the version used (8845) (the slight differences apparently come from the small mismatch in the times of the snapshots). One could run the simulation with higher time resolution and better match the snapshot, and then use the snapshot from the latest version of MESA. Or it could be simply mentioned that the results are completely consistent and just use the original MESA 8845 RGB profile for consistency with earlier papers. It does not make a difference so might as well do the latter.
6. I also ran into a problem when doing test runs, which is that one cannot add much refinment around the secondary at t=0 because by doing this one asks for more resolution than contained in the initial profile. It is easy to increase the resolution in the ambient, but the extra refinement region (i.e. particle buffer) around the secondary could overlap the envelope if the separation between the secondary and RGB radius is smaller than the particle buffer size. To increase the resolution in the envelope profile, one needs to rerun MESA at higher spatial resolution. Since we know from the AGB paper (see Fig B1) that the orbital evolution is insensitive to the resolution around the secondary before the first periastron passage, it is unnecessary to refine around the secondary at a high level right from t=0. So I did not bother rerunning MESA to get higher spatial resolution in the profile. However, if we, for example, want to add even one level of AMR in the outer envelope from t=0 in the future, we would need to rerun MESA with higher spatial resolution and reproduce the initial condition. Something to keep in mind.

Next steps

• Reproduce EoS tables using mean atomic weight of metals as 17.35 instead of 18 (differences will be minor but might as well).
• Migrate new code to Stampede and test (produce 1 frame for run 207 and make sure it agrees with bluehive result).
• Prepare the new fiducial run with gamma=5/3 EoS (increased resolution, etc.) and perform test runs (This was named Run R1 in the ​Plan for EoS project presentation).
• Prepare new EoS run where recombination energy is completely thermalized (Run R2) and perform test runs.
• Prepare new EoS run where recombination energy is completely radiated away (Run R3—this will involve a bit of coding in astrobear—how to proceed?) and perform test runs.

This run is a repeat of one of my previous runs (two identical heavy jets), except with more controlled refinement:

most of the run only goes to level 2, but there is a sphere of radius 0.5 in the center that goes to level 4 and a larger sphere of radius 1.0 that goes to level 3. This allows us to see a higher resolution in the region of collision without wasting computational resources where they are not needed.

This is a run for a 1D Jet with cooling

Labels on the Axes will be added soon. For now, red is log of temperature, blue is log of density, and green is velocity divided by 100

I also read "Numerical Analysis of the Dynamic Stability of Radiative Shocks" (Strickland and Blondin 1995)

Update 8/16: Labels have been added for axes, but still need to add colour key

Charge Exchange

Postprocessing done on the collisional ionization runs. Solar wind run is going well:

MHD

+----------------------+------------+--------------+------------------------------------------------------------------------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|                      |            |              | Strong Field (4.5x10^1^ G)                                                                                 | Weak Field (4.5x10^-4^ G)                                                                                                                                                      |
+======================+============+==============+============================================================================================================+================================================================================================================================================================================+
|                      |            | No wind      | Immediate pressure protections, but they slowly decrease. By the time the run ended (frame 10), none left. | Runs well. Slight acceleration, moving at max of ~5 km/s at 21 frames. Possibly from boundary condition?                                                                       |
|                      | No gravity +--------------+------------------------------------------------------------------------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|                      |            | Stellar wind |                                                                                                            |                                                                                                                                                                                |
| Non-corotating frame +------------+--------------+------------------------------------------------------------------------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|                      |            | No wind      | Pressure protections all over the place after a few frames, but they also disappear just before frame 10.  | Runs well. Accelerates toward star, bending field lines. 20 frames.                                                                                                            |
|                      | Gravity    +--------------+------------------------------------------------------------------------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|                      |            | Stellar wind |                                                                                                            |                                                                                                                                                                                |
+----------------------+------------+--------------+------------------------------------------------------------------------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|                      |            | No wind      | Pressure protections after a few steps, still many at frame 9.8 (furthest it got).                         | Runs well. Accelerates as expected, fixed at left & top because of boundary condition. 20 frames.                                                                              |
|                      | No gravity +--------------+------------------------------------------------------------------------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|                      |            | Stellar wind |                                                                                                            |                                                                                                                                                                                |
| Corotating frame     +------------+--------------+------------------------------------------------------------------------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|                      |            | No wind      | Pressure protections after a few steps. Only 1 frame.                                                      | Runs well. Accelerates outward (I think as expected? Orbital velocity of planet is higher than Keplerian velocity of ambient). 4 frames (for some reason only ran on 4 cores). |
|                      | Gravity    +--------------+------------------------------------------------------------------------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|                      |            | Stellar wind |                                                                                                            |                                                                                                                                                                                |
+----------------------+------------+--------------+------------------------------------------------------------------------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+

There still appears to be a problem with the stellar wind being injected. I'll check it without fields to make sure I can get it working properly. Also, the strong wind runs are looking suspiciously similar. I may run clean runs of everything for comparison purposes.

Run 1

Here are the last pressure protections from the simplest run:

q=      0.73E-04    -0.11E-01     0.37E-03     0.12E-01     0.39E+01     0.79E+00    -0.25E-01    -0.18E+01     0.00E+00     0.36E-04
Info%q=      0.73E-04    -0.11E-01     0.37E-03     0.12E-01     0.39E+01     0.79E+00    -0.25E-01    -0.18E+01     0.00E+00     0.36E-04
pos=       0.51875000E+02      0.11250000E+01      0.11125000E+02
i,j,k,mx,my,mz    20    45    25    20    53    20
Info%level     0
tnow=       0.10220561E+01
 Temp is   -17241.5756049894     
q=      0.73E-04    -0.11E-01     0.37E-03     0.12E-01     0.39E+01     0.79E+00    -0.25E-01    -0.18E+01     0.00E+00     0.36E-04
Info%q=      0.73E-04    -0.11E-01     0.37E-03     0.12E-01     0.39E+01     0.79E+00    -0.25E-01    -0.18E+01     0.00E+00     0.36E-04
pos=       0.51875000E+02      0.11250000E+01      0.11125000E+02
i,j,k,mx,my,mz     0    45    25    20    53    20
Info%level     0
tnow=       0.10220561E+01
 Temp is   -17241.5756071120     
 Pressure protection encountered at           20          45          25
 Pressure protection encountered at            0          45          25

EoS stuff

Work done

1. More plots and analysis of MESA RGB
2. Downloaded and installed latest version of MESA (12778 instead of 8845). Ran "same" MESA simulation with this version and compared to old sim.

Results

1. Initial conditions
   1. ​RGB energy profiles.
   2. ​mass fractions (direct MESA output).
   3. ​number densities of ions, neutrals and electrons. Oxygen shown as an example but can also plot C, N, Ne, Mg.
   4. ​neutral fractions (direct MESA output).
   5. ​mean atomic weight.

Discussion

• Gas thermal energy + recombination energy (+Prad which is negligible) almost equals to internal energy but a small decrement that grows with radius, becoming sizable right near surface.
• Went to trouble of including metals but rather negligible!
• Used Saha equation with T—>0 partition functions, which are just the degeneracies of states (so positive integers). In reality the partition functions can be much larger at the temperatures found in the star >~10⁴ K. And their ratios, which is what is important for the Saha equation, can change appreciably. One can estimate these from atomic databases (NIST, ViZieR) but numbers can also depend on electron pressure…complicated. What effect would this have on the recombination energy?
• Found that the average atomic mass of metals is 17.35. This is an input into Yisheng's code that produces the tabuler EoS because it is a parameter for the ideal EoS part of the table. He used 18 and anyway, the ideal EoS part of the table is rarely utilized.
• Found that new MESA version 12778 gives very similar results. Should we go with new MESA RGB model or stick with the old one?

Next steps

Initial conditions

• Estimate recombination energy using partition functions from NIST and/or ViZier databases and see what difference this makes.
• Learn more about the EoS to understand what other physics contributes near the surface (optically thin region).
• Put all the MESA profile figures in a pdf, with figure captions, for reference (to be used throughout the EoS CE project).
• Write a short section on EoS for the paper while it's all fresh, with help from Yisheng and Jonathan.

Simulation

• Reproduce tables using mean atomic weight of metals as 17.35 instead of 18 (differences will be minor but might as well).
• Steps mentioned in last blog post.

Idea for paper

• I had the idea of subtracting out the recombination energy (physically it is escaping by radiating away) and seeing what effect this has. Results will be at least somewhat different from ideal gamma=5/3 EoS. But it would also be interesting to do a run where we subtract ONLY the H recombination energy, say, to see what happens if recombination energy of hydrogen leaks out, but that of helium (and metals) is retained, to constrain the importance of H rec energy vs He rec energy.

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)

I did a run which repeated one of my previous runs but at a higher resolution. The patches for different levels of refinement are shown here

However I need to learn how control refinement better (right now I only know how to use levels of refinement)

Here is the result of the run

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.

Coupled EBM

​I Made a Poster for the TechnoClimes Meeting this Week

Fermi Project

​Click Here for The Most Recent Update