Mass unbinding

​Notes comparing unbinding using different criteria for "unbound"

Energy conservation and effective core profile in simulation

​Notes on algorithm that conserves energy and effective core profile in simulation

Energy paper erratum

​Draft of erratum for Paper II to correct unbound definition and initial PE profile figure

Revision

Revised paper is ready to submit again. ​The last version

Jet project - plan for this week

Runs:

1. Rerun the super jet model, fix tracer, add tracer for ambient
2. Rerun the no-jet fiducial, M2 = 1 Msun, use my jet modual, this is like CE run 152 (or 143, but don't change softening radius), use the initial orbit conditions for the jet runs, add tracer for ambient
3. Half-solar mass companian, fiducial, no jet run
4. (Maybe) Restart, add jet from middle of the fiducials

Analysis:

• Orbits, average separation plot. see the eccentricity formula form Luke's blog post
• Mass, bind vs. unbind, find a best way to present this.
  • maybe something like Fig. 8 in the AGB paper, where more mass gets elevated but still not unbind
• Energy, use data from the totals, do post-process if necessary
  • maybe analyze drag force as well, but orbit plot also shows the effect of drag force
• Present figures similar to the AGB paper
  • Fig. 6, 7, 8
• Update the separation plot
  • use CE run 152 for high mass fiducial, which doesn't change the softening radius

Calculations:

• how much energy do we pump into the simulation (jet)?
• Does the jet material fall back? (tracer plot will also help)

Extending the CEPN simulations

​Link to a summary document

Starting the Rad Hydro

​my plan to start

​Charge Exchange Paper

Sent out. Reviewing together Thursday 1 PM.

MHD tests

Something weird happened to the run log on BlueHive. It only has 4 of the ~24 frames that were actually completed by the last run. It's looking much better than the previous couple of tests, though:

This Test (DIFF_ALPHA2 = 0.01)

No Diffusion

High Diffusion (DIFF_ALPHA2 = 0.1)

Acceptable? Does the diffusion need to be even lower? Is it safe to turn off the diffusion, or give it another 10-15 frames?

GJ 436b

Got a bunch of test frames for a short early-time test (for the previously-mentioned pressure protections). The last frame copied so far shows the planet getting "hairy," for lack of a better word.

Spring code

Nearly written a full test suite. A few left in shapes, output tests, and then it'll be good to merge into the master branch.

A new way to calculate EOS

• Using another routine in MESA
• Compare MESA-r11701 EOS and MESA-r11701 Profile: ​https://www.pas.rochester.edu/~yishengtu/research_files/CEE_gp_meeting_05262020/my_profile_eos_compare.png
• Compare MESA-r11701 EOS and MESA-(older, from Luke) Profile: ​https://www.pas.rochester.edu/~yishengtu/research_files/CEE_gp_meeting_05262020/LC_profile_eos_compare.png

Inverting new EOS

• new EOS is not fully monotonic
• Numerical difficulties outside table range and dragon region
• Compare inverted EOS and original EOS, plotted percentage difference
• Energy: ​https://www.pas.rochester.edu/~yishengtu/research_files/CEE_gp_meeting_05262020/EDT-TDE.png
• Pressure: ​https://www.pas.rochester.edu/~yishengtu/research_files/CEE_gp_meeting_05262020/PDT-TDP.png

Directly compare new EOS and old EOS

• P(DT) ​https://www.pas.rochester.edu/~yishengtu/research_files/CEE_gp_meeting_05262020/old-newPDT.png
• E(DT) ​https://www.pas.rochester.edu/~yishengtu/research_files/CEE_gp_meeting_05262020/old-newEDT.png

Convergence Plot of Uniform Resistivity Simulation

Plotted the rate of change of Amplification as a function of Max field in each frame for all cases with uniform resisitivity structure and 2 winds (1000 #/cc and 300 #/cc).

We can see that those with high resisitivity start off with Amplification lower than the theoretical predicted max amplification limit and then increase towards their final steady state, while the ones with high conductivity start at high amplification and then decrease to steady state. (Kinks are due to restart with change in boundary conditions and buffer zones.)

New Lunar Proposal

We are trying to submit a proposal for study of lunar samples. Need to simulate the amplification for the present day moon, moon at 3.8 Ga and moon with thin atmosphere at 3.8 Ga.

Most recent experimental data of moon resistivity profile (from Apollo 12):

Approximation of above data into models from recent data (Shimizu et al. 2013):

Current Profiles running on BlueHive (very slow convergence):

Miscellaenous

Will soon be on my 3rd laptop in as many months. In general interest: ​https://xkcd.com/2083/

HEDSA Talk

HEDSA is organizing a series of talk about the future of Frontier Plasma Science. There will be a talk given by Carolyn Kuranz next week, Monday June 1 at 4pm EDT and I (Prof. Pierre Gourdain) was hoping that you could join the discussion. It is also very important to get students involved because we will be discussing something directly impact you. So, I would encourage you to reach out to all graduate students you know (get that Rolodex smoking), so they can ask questions using the Google form link ​https://forms.gle/4DupvDRRJKBn8p537

(Email me for link to the Zoom meeting.)

AGB paper

The current version of the AGB paper (with latest changes in red)

1. ​Corrected paper

Some notes to illustrate the reason for the changes

1. ​Notes on potential energy and potential energy density

RLOF and CE sims in general

Textbook chapter on the two-body problem showing that eccentricity can be computed from total energy and angular momentum

1. ​Notes on two-body problem including formula for eccentricity

These notes can be found at ​https://tmg-web.lehman.edu/faculty/anchordoqui/chapter25.pdf
I got the idea from this preprint, which is of interest: ​Kuruwita, Federrath & Haugbolle 2020

Energy conservation in the code

Some notes (now wrong I realize because Eq. 1 is wrong!) which prompted me to correct the ½ factor — extending energy conservation routine to include particles should be trivial

1. ​Notes on algorithm that conserves energy

Relevant paper (just putting it here for easy access)

1. ​Jiang+2013

​Charge Exchange Paper

Postprocessing finishes tomorrow or Wednesday. Still plan to have draft done by end of week.

GJ 436b

Having some issues with pressure protections. I suspect it's related to the planet-ambient interface, so trying a test run now.

Spring code

Started writing tests (using googletest framework). Going well so far. Found plenty of bugs.

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

Effect of Resolution, Interpolation and Resistivity Smoothing

The resistivity of the asteroid was being raised to high values from near-zero ambient value. Adding a exponential buffer made the transition of maximum field (in the zone just outside the surface) less step.

The default interpolation scheme (constant) was being used which was causing steep changes in the field across AMR levels. Changing it to 1 (minmod) makes the slope gentler.

(Yet to run Lv 4 simulations for comparison.)

Finalized Movie and Figure for Paper

External Image editing had to be involved for scale bars and labels. Draft was submitted today.

​Charge Exchange Paper

Postprocessing is done for the high-wind sim. I believe the low-wind postprocessing should be done by the beginning of next week, so I'm hoping to have a complete first draft of the paper to send around to everyone by the end of next week. See ​last week's post for a table comparing the major charge exchange papers. I didn't include Ekenbäck et al 2010, since they use a baked-in magnetosphere - don't think that would be a very useful comparison.

Exo3

I'm moderating. Once charge exchange paper is sent around I'll start running magnetic field tests again, so we can hopefully have some results for the poster.

Spring Code

I've integrated OpenMP parallel for loops everywhere I think they'll be useful. I realized this weekend that although Rebound has MPI built in, I'm not sure if the spring code add-ons will work with it. I need to dig a bit to figure out how the particles are distributed across processors.

Also on my list are to add tests, make sure all of the functions present are actually used and delete the unused, and write up some documentation.

Clump

• Tried ideal gas gamma = 1.1 (wind speed 1e8, sound speed ~1.1e7)
• upper: table eos, lower: ideal gas with gamma = 1.1
• ​https://www.pas.rochester.edu/~yishengtu/research_files/CEE_gp_meeting_05042020/Template/out_combined/combined.gif

Convective instability

• stable if |dT/dP| < |dT/dP|_ad
• If blue line above red line, the star in convectively stable. Otherwise unstable
• ​https://www.pas.rochester.edu/~yishengtu/research_files/CEE_gp_meeting_05042020/StarOnGrid/convective_instability.png
• Notebook of math ​https://www.pas.rochester.edu/~yishengtu/research_files/CEE_gp_meeting_05042020/StarOnGrid/CEE_convection_notebook.pdf

StarOnGrid simulation

• A simulation with just a star sitting there (left: ideal gas, right: table eos)
• Radial velocity (v_r), red is radially outward and blue is inward ​https://www.pas.rochester.edu/~yishengtu/research_files/CEE_gp_meeting_05042020/StarOnGrid/combined_vr_figs/combined.gif
• There seems to an infall, leading to a wave propagating outward
• Slow motion of what happened between 0th and 1st frame ​https://www.pas.rochester.edu/~yishengtu/research_files/CEE_gp_meeting_05042020/StarOnGrid/ideal_slow_figs/movie.gif
• Maybe an change in resolution issue?

Future step

• email MESA (later this week or next week)
• check the resolution issue
• perhaps we are ready to run CE?

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