GJ436b

Charge Exchange

Solar wind preview postprocessing is done. Here are the strong and intermediate winds with hot and cold split:

Pathways Proposal

Deadline was extended to this Friday (I get the impression that they feel Frontera is underutilized so far), so taking the opportunity for a few quick questions:

Fields:

120 - Astronomical Sciences (Primary) 510 - Atmospheric Sciences (Secondary) 130 - Physics (Secondary) 122 - Planetary Astronomy (Secondary)

Jonathan: University Research Staff or Faculty?

Grants: So far I've included the DoE grant, but it's MHD-focused, so maybe just barely fits the criteria: "Is the proposed work supported by a grant or grants that have undergone review for intellectual merit and/or broader impact (if relevant to the allocation)? If so, is the allocation request consistent with the objectives of the supporting grants and is the scale of resource use commensurate with the level and purpose of support?" and "the stated scientific objectives must match or be sub-goals of those described in the listed funding award(s)". Could perhaps describe as a sub-goal. Don't think we can include the NSF grant at all, based on this.

Charge exchange

Three more frames with startup allocation, but only have 103 node-hours left, so definitely can't get any more frames with it. May try to do some profiling if I have time. Here's the current state:

Note of fun - the simulation is surprisingly robust to a large empty cube suddenly appearing in the center of the planet.

Frontera Proposal

Strong scaling test is done for 8-128 nodes. It looks like the proposal should be a single document, so I've combined the scaling document and the previous supplement main document. Still needs a little work, but I'll try to send it around by end of day for review.

One sticking point: the steady-state run is about 12 times faster than the current state of the simulation. How to best request the appropriate number of hours? Pathways allocations are meant to help with scaling up to make better use of Frontera (see description below) - perhaps frame single run as test case for larger request that combines both charge exchange and radiation pressure (and B fields?)? And possibly profiling time? Would still like to see if the line transfer can be improved.

Pathways

Experience suggests that not all user teams can seamlessly make the transition from more general-use resources, such as Stampede2, to a system that emphasizes larger scale computing. It is here that latent scaling bugs (such as race conditions, hard-coded array bounds, and so on) are often encountered for the first time. The “Pathways” allocations track is for projects that believe they are ready to begin scaling up, but have not yet fully tested their applications at scale. Less stringent scaling data will need to be provided to receive a pathways allocation, and the total award size will be smaller, with a focus on preparing to reach the LRAC scale.

Here's the statement of requirements, for reference:

Proposal Format

LRAC, Pathways, and LSCP Submissions Requests

When writing your LRAC request, be clear and concise. We strive to have domain experts review every request, but they may not have deep expertise in your specific subdomain. Someone outside of your area should be able to understand the scientific objectives and understand why the chosen technique is preferred over another.

Renewal requests require less documentation, as described later (jump link to section).

The documents required of a new LRAC request ensures reviewers can effectively determine how each request satisfies the Frontera allocation Review Criteria. LRAC requests are limited to 10 pages.

Research requests must include a well-documented resource-use plan that describes how the requested allocations are necessary and sufficient to accomplish the project’s research objectives. An effective resource-use plan must address the Frontera allocation Review Criteria and, in particular, must include the following elements:

Scientific Background Research Objectives and specific questions to be pursued Resource Usage Plan to achieve the research objectives Justification of the allocation amounts for all resources and resource types Access to other CI resources and why those resources are not available or sufficient for the work proposed in this request

Scientific Background and Support

Succinctly state the scientific objectives that will be facilitated by the allocation. The existing merit-reviewed supporting grants should be listed and briefly described; please include the end date for each grant listed. Requests with merit-reviewed supporting grants will not be subject to further scientific review by the LRAC; however, the stated scientific objectives must match or be sub-goals of those described in the listed funding award(s). The description of the objectives should be sufficient to allow the reviewers to assess the resource usage plan.

Research Questions

Identify the specific research questions that are covered by the allocation request. Within the context of the scientific background and supporting grants, these objectives and questions should be stated so that the reviewers can understand how elements of the resource plan will contribute to relevant answers. If you consider the allocation request in terms of computational experiments, the research objectives and questions define the experiments to be conducted.

Resource Usage Plan

The bulk of the document should focus on the resource usage plan and allocation request. Inadequate justification for requested resources is the primary reason for most reduced or denied allocations. Once again, the PI should keep in mind the Frontera allocation Review Criteria when describing and justifying the choice of resources and allocation amounts.

Justifying Allocation Amounts

This section should contain code performance timings, resource usage details, and scaling information to support the calculation of the resource request. Ideally, the code performance data should be provided from benchmark runs on the resource requested, using the model configuration(s) needed for the computational plan proposed and demonstrating the scaling efficiency to the job size(s) planned. If applicable, the panel may accept well-justified performance data from architecturally equivalent systems. PI’s should contact the Frontera Allocations Coordinator to obtain a small allocation for benchmarking or to discuss non-Frontera performance data to be used in their submission.

The justification of allocation should take the quantitative parameters from the resource usage plan and combine them with benchmark cost information to calculate the allocation needed. Where possible for computational experiments, the justification should tabulate and calculate the costs of conducting each experiment for each resource and resource type.

Disclosure of Access to Other Compute Resources

Because of the high demand for NSF-funded resources, Frontera and the LRAC closely consider whether the PI has access to other, less constrained resources on which the proposed work could be conducted. All PIs must describe their local resources and other non-Frontera resources available to them, including both local and other national resources.

Page limit: 10 pages

References

[Note: probably unnecessary]

A PI may use this OPTIONAL document to separate a lengthy bibliography of cited work to take full advantage of the page limits in the Main Document. This bibliography is not the publications resulting from prior Frontera support, which should be entered into the Frontera publications database, but rather other citations referenced in describing or supporting the intellectual merit of the proposed work or the appropriateness of the proposed approach for addressing the research objectives.

Page limit: No limit.

Document Formatting

While readability is of greatest importance, documents must satisfy the following minimum requirements. Documents that conform to NSF proposal format guidelines will satisfy these guidelines.

Margins: Documents must have 2.5-cm (1-inch) margins at the top, bottom, and sides. Fonts and Spacing: The type size used throughout the documents must conform to the following three requirements:

Use one of the following typefaces identified below:

Arial 11, Courier New, or Palatino Linotype at a font size of 10 points or larger; Times New Roman at a font size of 11 points or larger; or Computer Modern family of fonts at a font size of 11 points or larger.

A font size of less than 10 points may be used for mathematical formulas or equations, figures, table or diagram captions and when using a Symbol font to insert Greek letters or special characters. PIs are cautioned, however, that the text must still be readable.

Page Numbering: Page numbers should be included by the submitter.

Frontera

Have the planet problem compiling, but running into runtime errors that appear to be related to the new Intel Fortran runtime library version. Specifically,

forrtl: severe (257): formatted I/O to unit open for unformatted transfers, unit 76, file /work/05515/adebrech/frontera/ChargeExchange_high/global.data
Image              PC                Routine            Line        Source             
libHYPRE-2.18.2.s  00002B3589B43B86  for__io_return        Unknown  Unknown
astrobear          00000000019B4FBE  Unknown               Unknown  Unknown
astrobear          00000000019704BE  amrcontrol_mp_rea         218  amr_control.cpp.f90
astrobear          0000000001960424  amrcontrol_mp_amr          72  amr_control.cpp.f90
astrobear          00000000019837E8  scrambler_IP_astr          63  scrambler.cpp.f90
astrobear          00000000019835A6  MAIN__                     51  scrambler.cpp.f90
astrobear          000000000040E592  Unknown               Unknown  Unknown
libc-2.17.so       00002B358BF26555  __libc_start_main     Unknown  Unknown
astrobear          000000000040E4A9  Unknown               Unknown  Unknown

Relevant code:

! Open data file (quit if not successful).
OPEN(UNIT=GLOBAL_DATA_HANDLE,FILE=GLOBAL_DATA_FILE,FORM='formatted') ! I added the FORM specifier to see if that would fix it

! Read in global data values.
READ(GLOBAL_DATA_HANDLE,NML=GlobalData)

Haven't gotten past this one yet, but I suspect something similar will happen with future read attempts.

When I attempt to INQUIRE about GLOBAL_DATA_HANDLE I get an even stranger message:

inquire(UNIT=GLOBAL_DATA_HANDLE, access=acc, name=filename, form=formatType, exist=exists)

forrtl: severe (257): formatted I/O to unit open for unformatted transfers, unit -1, file /proc/265613/fd/1
Image              PC                Routine            Line        Source             
libHYPRE-2.18.2.s  00002B69A5BBBB86  for__io_return        Unknown  Unknown
astrobear          00000000019C6A58  Unknown               Unknown  Unknown
astrobear          000000000196FF02  amrcontrol_mp_rea         219  amr_control.cpp.f90
astrobear          0000000001960424  amrcontrol_mp_amr          72  amr_control.cpp.f90
astrobear          0000000001983AD0  scrambler_IP_astr          63  scrambler.cpp.f90
astrobear          000000000198388E  MAIN__                     51  scrambler.cpp.f90
astrobear          000000000040E592  Unknown               Unknown  Unknown
libc-2.17.so       00002B69A7F9E555  __libc_start_main     Unknown  Unknown
astrobear          000000000040E4A9  Unknown               Unknown  Unknown

More investigation leads me to conclude that there's something wrong with I/O in general. I've tried re-chowning and re-chmoding all of the .data files, but no change.

print *, 'open status is ', iErr

forrtl: severe (257): formatted I/O to unit open for unformatted transfers, unit -1, file /proc/271281/fd/1
Image              PC                Routine            Line        Source             
libHYPRE-2.18.2.s  00002AD3EA073B86  for__io_return        Unknown  Unknown
astrobear          00000000019C6828  Unknown               Unknown  Unknown
astrobear          000000000196FE81  amrcontrol_mp_rea         218  amr_control.cpp.f90
astrobear          0000000001960424  amrcontrol_mp_amr          72  amr_control.cpp.f90
astrobear          00000000019838A4  scrambler_IP_astr          63  scrambler.cpp.f90
astrobear          0000000001983662  MAIN__                     51  scrambler.cpp.f90
astrobear          000000000040E592  Unknown               Unknown  Unknown
libc-2.17.so       00002AD3EC456555  __libc_start_main     Unknown  Unknown
astrobear          000000000040E4A9  Unknown               Unknown  Unknown

If no one else has thoughts, I'll shoot TACC support another message.

Charge Exchange

Used up the allocated 20% on Stampede2. Gotten 27 frames so far; don't think we're at steady state, but haven't been able to check yet. Have 5000 SUs on Frontera now to get to somewhere around 35-40, hopefully.

Upcoming Conferences

​Exoplanet Demographics, free?, 11/9-12

​AAS 237, $150, 1/11-15

​Numerical Modeling in MHD + Plasma, $100, now

Charge Exchange

4 more frames. Second most recent frame:

MHD

line_transfer branch is up to date for both the PlanetaryOutflow (radiation-driven) and OutflowWind (Parker-ish wind) modules, though not everything implemented in the first is implemented in the second (various non-radiation rampings and field-related things).

To-Do

• Test B field in OutflowWind with no planet
• Implement and test embedded line transfer
• Implement and test planet import

Misc

• Metastable HeII line
• Profile line transfer on stampede
• GJ 436b

Need more computing time for those last two

Charge Exchange

Making better progress. Still not quite at steady state after 20 frames.

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?

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).

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

Exo3

Posters. Mine is supposed to be MHD, though it can still be edited. I think we may have gotten enough to make an interesting poster, though, especially with the last frame or two of the most recent run. Current plan is to discuss Eric's theory about the field guiding the wind in a more LOS direction in the context of original parameter space run and new results.

MHD tests

Truly uniform, constant magnetic fields. Initial conditions are the same in each pair of strong/weak field runs, barring the magnetic field strength. All simulations are using outflow-only boundary conditions.

Haven't run the simulations including a stellar wind yet, a little more difficult to double check the implementation of those.

+----------------------+------------+--------------+--------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|                      |            |              | Strong Field (4.5x10^1^ G)                 | Weak Field (4.5x10^-4^ G)                                                                                                                                                      |
+======================+============+==============+============================================+================================================================================================================================================================================+
|                      |            | No wind      | Immediate pressure protections. No frames. | 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      | Immediate pressure protections. No frames. | Runs well. Accelerates toward star, bending field lines. 20 frames.                                                                                                            |
|                      | Gravity    +--------------+--------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|                      |            | Stellar wind |                                            |                                                                                                                                                                                |
+----------------------+------------+--------------+--------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|                      |            | No wind      | Immediate pressure protections. No frames. | Runs well. Accelerates as expected, fixed at left & top because of boundary condition. 20 frames.                                                                              |
|                      | No gravity +--------------+--------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|                      |            | Stellar wind |                                            |                                                                                                                                                                                |
| Corotating frame     +------------+--------------+--------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|                      |            | No wind      | Immediate pressure protections. No frames. | 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 |                                            |                                                                                                                                                                                |
+----------------------+------------+--------------+--------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+

Images of runs, numbered in order of appearance in table above:

Run 2

Run 6

Run 10

Run 14

Charge Exchange

Solar-type wind crashed with MPI issues before completing a single frame. I've halved the number of cores to reduce the messaging, which should fix that.

Completed 3 (of 5 planned) frames of strong wind with collisional ionization:

MHD

Still having issues with both the restart and the current fresh run.

The restart is starting properly now, but still has immediate pressure protections (not unexpected) and very high CFL to start, so it almost immediately freezes.

I'm trying to restart the fresh run from earlier frames with a higher diffusion coefficient. Running into the same issue with high CFL conditions, but this time it seems linked to the higher diffusion coefficient. Going back to an earlier frame reduces the equilibrium predicted wall time, but still trying to figure out where I need to go back to in order to get a reasonable value. Removing the maxsolverspeed helps with, but does not entirely solve, this problem. The pressure protections here are mostly in ghost zones, and pretty much entirely at various grid intersections.

Charge exchange

Solar-type wind going slow, pretty much as expected. Should speed up as the shock relaxes. But appears to be running well.

Charge Exchange

Collisional ionization appears to be working. Definitely a strong upper bound without it…

Frame 300, old run

Frame 299, new run

I'll queue the GJ 436b run (as soon as I figure out why the slurm script isn't working…) and the new HD 209458b run.

MHD

Officially stuck, so I'll turn the diffusion up some more. See ​last week's post for latest frame.

I've also set up a run to try restarting from the final frame of the original run with an extremely high field. Need to figure out how to reload a hydro run with MHD, and add to total energy when I add the field.

MHD

Only got one frame last run, still getting tons of pressure protections (despite not being close to MinTemp anywhere, as far as I can tell). Perhaps diffusion needs to be turned up a bit? But things are starting to look very interesting. Also, realized the stellar field wasn't being applied correctly at the left boundary, so I've fixed that.

GJ 436b

Slurm script seems messed up. Need to fix. Otherwise working well.

​Charge Exchange Paper

Made most of the changes we discussed at the meeting. Waiting on results of collisional ionization test to start the 1.2x solar mass loss rate run, with more solar-wind-like properties.

These runs are currently using the parameters from Run 2 of the ​parameter space paper, to wit:

The initial B field excludes the center of the planet.

Here are some relatively current movies:

Fixed left boundary

B field

rho, temp, P

beta, no ram

beta, with ram

​Charge Exchange Paper

Haven't gotten any new comments yet. Still working on the collisional ionization run on Stampede.

MHD

Still lots of pressure protections, but running reasonably ok (> 1 frame per day). Looks like the high-temp vacuum is starting to smooth out, so hopefully we're over the hump.

GJ 436b

Still no issues at equivalent of frame 4 (where it crashed with pressure protections before frame 1 last time). If it makes it to t=1 (frame 10), I'll just rename the files and have it restart from there.

Jobs, etc.

Started updating my (long long form) CV based on suggestions from AAS, also finishing up the requirements for the university's teaching recognition. This includes a first draft of a teaching philosophy, which I think is sufficient for the requirements but needs some revision.

Also probably going to apply to that job Laura sent around - may want people to start too early, but at least some gov't contractors have very long hiring processes, and worst case it's good practice and I might (maybe) get some useful feedback.

Spring code

Lots of tests fail with OpenMP. Need to determine which fail because of bugs and which fail because of pseudorandom ordering. On the backburner for the moment.

​Charge Exchange Paper

Made the changes discussed, responded to Luca's comments. Have a job queued on Stampede to test changes with collisional ionization.

MHD

Pressure protections started immediately, so there must have been some in the missing portion last time, too. But it still seems to be running reasonably ok. Here's the current frame:

GJ 436b

Haven't hit a problem yet, it seems, so I'm continuing the current run to see where it does. If we do run into an issue, I'll run some tests with the outer boundary out farther, see if I can find a better balance between keeping the atmosphere hydrostatic and still having a decently-low ambient temperature.

​Spring code

Pretty much completely done with current to-do list. Merged into master branch on GitHub. Would still like to figure out how particles are distributed using MPI (how is an essential tree implemented?) and possibly adapt the spring portion to be a hybrid code, but OpenMP works well enough for now (though I'm segfaulting in the tests while using OpenMP, so need to fix that).

​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.

​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.

​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.

Charge Exchange

Here's the split of absorption between hot and cold neutrals:

Shows pretty much what the theoretical calculation does - the hot neutral density just isn't high enough to get significant absorption. You can see a lot more broadening in the hot case, the total absorption is just very low compared to the cold.

Papers

Summary: After reviewing all of the papers, it's starting to look like everyone can fit GJ 436 but not HD 209458.

(Not sure why only the first attachment link is working…. Open post to download attachments.)

Spring code

Basically done. A few more problems to update, some minor refactoring and upgrades for the spring part (particularly, need to figure out the correct places to parallelize with OpenMP).

Charge Exchange

With a theoretical optical depth as a function of frequency of

\tau = \sigma_{\nu_0} n_{HIhot} d \frac{c}{\nu_0} \sqrt{\frac{m_H}{2 \pi k_B T}} e^{{-\frac{(\nu - \nu_0)}2 c^{2 m_H}{2 \nu_0}2 k_B T}}

0.125704 e^{{-8.92598\times10}{-25} (\nu - 2.47\times10¹⁵⁾2)

Which gives an optical depth at \nu_0 of .126, an optical depth at \pm 100 km/s of 0.061.

I've attached a Mathematica notebook with my work.

I've also made it possible to output both the unsplit and split (as separate files) Lyman-alpha absorption from postprocessing. Should be able to get some preliminary results in a day or two.

Spring Code

Split into 10 modules now, approximately 6ish left to refactor. May be able to finish by the end of the week. Almost certainly by the end of the month.

​Charge Exchange Paper

Figures and captions (except for postprocessing) done. Also basically done with methods. Need to talk about other papers in intro still (which? of Khodachenko19, Khodachenko17, Shaikhislamov16, Christie16, Tremblin13, Ekenback10, Bourrier16).

Postprocessing: Completed up to frame 276 for the split-species high wind. Completed up to 226 for the split-species low wind. High wind is taking 7-8 hours per frame, so should be done late next week (~10-12 days). Low wind takes closer to 10 hours per frame, but hopefully can be done by the end of the month.

MHD

Test started recently. Will do some visualization checks soon, but appears to be running (slowly, but well). Will have to stop to finish the charge exchange postprocessing.

Charge Exchange

Low wind finished. Here are movies of both (missing last few frames of low wind):

High Wind

Low Wind

Spring Code

Converted Alice's code to C++ to make some of the math significantly easier to express. Still working on sorting and refactoring the majority of the code. There also appear to be dozens of unused routines - need to check on those.

Misc

Contacted by Research Outreach, paid magazine publication service. I'm assuming we're not interested, but wanted to make sure.

MHD

Diffusion coefficient still appears too high. However, when set at what I think it should be, I run into min timestep issues (Atma appears to be having a similar issue). I first believed this to be related to the addition of a DIFF_ALPHA2 term for the maxspeed:

maxspeed(Info%level) = max(maxsolverspeed(Info%level), max(maxspeed(Info%level), 10d0*DIFF_ALPHA2/levels(Info%level)%dx))

But looking again, the lower DIFF_ALPHA2, the lower that term should be, so I don't believe that's the problem.

Charge Exchange

Low-wind run finished over the weekend. Working on the postprocessing now, along with the new separate-temperature postprocessing for both.

AMR Line Transfer

Buffered sends with a stack of rays is a significant improvement over blocking and looping. Tests on 120 cores:

I've merged the AMR and projected versions of line transfer into one branch, so it's possible to choose between them with a single switch. Default is non-AMR. During these tests, noticed that AMR version is ~2x faster than non-AMR on 24 cores. Hasn't affected execution time for non-AMR version.

Still probably a way to improve the ray processing order to minimize dependencies, but will take significant thought.

Charge exchange

To get a reasonable time for averaging, need about 9 more frames. Almost there…

MHD

Testing a diffusion (DIFF_ALPHA2) run. Still having pressure protections, but entirely in ghost zones now. Not clear yet if it will improve runtimes if these protections are solved, or if it will behave well (seems like too much diffusion right now).

DIFF_ALPHA = 0

DIFF_ALPHA2 = 0.1

With Diffusion, Frame 39

Without Diffusion, Frame 27

Charge Exchange

Weak wind should be done this week. Here are the current figures for the high wind:

I find the differences between the initial state and the high wind interesting. The absorption is definitely more spread out, but I'm not sure enough for the significantly less-deep line center absorption. Need to compare total optical depths.

Weak wind

MHD

Results from fixed left boundary field. Left is original run, right is with fixed boundary conditions.

Radiation Pressure Paper

Proofs received and returned last week. Should be complete now.

Charge Exchange

Current frame of low stellar wind:

Think we've reached a steady state. Running 15 more frames, and double-checking the charge exchange rate.

High wind postprocessing should be done mid-week.

Charge Exchange

Preliminary figures for the strong stellar wind case:

There's a bug. I'm checking the impact.

Weak stellar wind:

At frame 242, out of room on stampede again. Transferring more radiation pressure runs to afds3 as soon as possible.

Radiation Pressure Paper

Accepted, finally. I'll fill out the paperwork and update the ArXiv version this afternoon.

Charge Exchange

High wind

Here's the final movie. Waiting for postprocessing job to start on BlueHive.

Low wind

MHD

Here's a movie up to the current frame:

We're getting some significant slowdown due to pressure protections along the boundaries of the wind. It's actually sped up over the last couple of frames — down to 2.5 months instead of 4 months.

What I believe is happening is that a vacuum is trying to form along some of the wind-ambient boundaries. Which leads to a couple of thoughts: first, have we ever tried to build vacuum handling into AstroBEAR? I know it's possible with the exact solver, so it may be something where we could have a failsafe trigger to use the exact solver at pressure protections instead of averaging or otherwise modifying the pressure. Second, would lowering the ambient density initially help with this? I.e., if it's already closer to vacuum, would it improve runtime?

Misc

Still waiting on a bunch of MHD and AMR line transfer test runs on BlueHive.

Radiation pressure paper

Revision submitted. Hopefully it'll be accepted soon and I can move the data off BlueHive.

MHD

Ramp vs No-Ramp

Slight differences. But the non-ramping version is working, so that's probably the better way to go. Working on some resolution tests now.

Charge Exchange

I have 10 more frames of this, but no room on BlueHive currently. The low-pressure wind is working (have 5 frames or so), but I believe I've pretty much used up my allotted portion of the allocation on Stampede2.

I will start the postprocessing on BlueHive, so that we can start looking at synthetic observations.

AMR line transfer

This may become useful for the higher field cases. In order to decrease the discontinuity between the planet and the ambient, we need higher resolution around the planet - but we can force the rest of the simulation to stay at e.g. level 4 maximum while refining the planet to a much higher level. However, this would be practically impossible to run in the "projected to max level" version of the line transfer routine. I'm working on updating and merging this into all of my branches.

Currently waiting on a timing test for the buffered nonblocking send with ray stack version, to see if that's an improvement over the last looping version (which may or may not have used IBsend during the test - if this isn't an improvement, probably have to check again).

At some point I would like to make a graph of the dependencies and see if there's a best way to prioritize the rays.

HPC summer school

Statements are due midnight tonight, AoE.

Exoplanets III

Talk abstract for MHD submitted. Deadline is this Friday, I assume we'll hear back Feb-Mar. Working on flight and hotel.

Misc

I've updated the development branch to use the MPI_F08 module, which allows type and bounds checking, makes iErr optional, and in general is both better and the only way to use MPI that's completely compatible with the Fortran standard. Sent a merge request on Gitlab.

I now have my own professional domain: ​alexcomputes.com (was tired of using my much less professional free domain on my CV).

MHD

See my page for ​movies of moderate-field runs.

Next steps: Fix field lines in left (& right?) boundary. Test starting from complete field. Attempt higher fields (will likely require high resolution at planet boundary - means need to spend some time integrating AMR line transfer).

Charge Exchange

Successfully replaced ambient for low-wind run, but running into pressure protections - trying to work through them now.

Misc

Exoplanets 3 registration is open. AAS dress code?

MHD

The bug was tracked down last week. The low-field test case now works (known working up to a stellar beta of 10^-2), so I'm experimenting with some values of the ramping to get the higher-field cases working.

Charge Exchange

Here's a movie up to frame 30:

And here's the latest frame:

I've implemented a check to make sure the stellar wind has a greater total pressure than the ambient, and if not to replace the ambient (or material at a density threshold approximately equal to the ambient density) with the stellar wind (which will be a much more significant shock to the simulation, I think, but probably can't be helped). Still testing this.

Scaling Test

Still fighting with the KNL nodes. I'm talking with Baowei to try to figure out why.

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.

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.

MHD

We've found the bug, in upwinded_emf. I'll be checking the implementation from ​Gardiner & Stone 2005, but for now it runs.

We now have a more physical issue, with infall from around the AMR boundaries which I believe is due to (magnetic) pressure gradients in the ambient material. Setting the total (thermal+magnetic) pressure to a constant equal to the far-field value creates more problems (pressure protections in the initialization), but I have one more idea to try (in addition to trying a different interpolation method, to see if the boundaries will smooth out that way). Jonathan pointed out there's also magnetic tension - I'm hoping that countering the pressure will be sufficient to make things stable long enough to start the wind.

Charge Exchange

Also started a low-mass-loss-rate stellar wind run. Will check and see how that one's running soon.

MHD

I think a CFL issue has been successfully ruled out (Jonathan, if you think a direct wave speed comparison would be useful, I can do that - but based on the tests already done, don't think they're necessary).

Is the prolongation of the aux fields still the most likely place for a bug? It may be time to sit down for a bug hunt.

Charge exchange

Still moving right along. Hoping we reach a steady state soon. Here's the latest frame. (At some point I'll get a movie, but waiting to hear back about the paper so we can hopefully archive the rad pressure run and make room….)

I meant to turn off the lower y boundary, turned off the upper one instead - been corrected, and shouldn't affect the simulation as it progresses (but I've queued the restart from an earlier frame anyway).

AMR line transfer

Set up a stack, still need to change ray tracing to use it. Think it should be a stack of rays rather than ray positions, which means the main ray tracing loop needs to be rejiggered again.

Once the MHD bug is worked out I should have time to go a little deeper into this, see if there's a better way to approximate the graph that represents the dependencies.

MHD

Believe current problems are related to prolongation of the aux fields. Because we lose information when going from face-centered to cell-centered B fields, we can't just reconstruct the "correct" (uncorrected) aux fields from q. I believe this is the relevant piece of code, in ProlongateParentData:

DO i=1,nDim
  DO j=0,r**nDim-1
    IF (MOD(j/r**(i-1),r)==1) CYCLE
    DO n=1,nDim
      l(n)=MOD(j/r**(n-1),r)  
    END DO
    ic(1:nDim,1)=1+l(1:nDim)-rmbc; ic(1:nDim,2)=Info%mX(1:nDim)-r+1+l(1:nDim)+rmbc;ic(i,2)=ic(i,2)+2
    ip(1:nDim,1)=mB(1:nDim,1)-mbc; ip(1:nDim,2)=mB(1:nDim,2)+mbc;ip(i,2)=ip(i,2)+1
    Info%aux(ic(1,1):ic(1,2):r,ic(2,1):ic(2,2):r,ic(3,1):ic(3,2):r,i) =     &
    Parent%auxChild(ip(1,1):ip(1,2)  ,ip(2,1):ip(2,2)  ,ip(3,1):ip(3,2)  ,i)
  END DO
END DO

Haven't looked to see if there are any obvious errors yet. Outer ghost zones were the source of the error in last test.

AMR line transfer

Thinking some more about prioritization. Setting up a stack isn't a bad idea - may be the next thing I try, since it should be easy.

In the truly ideal case, we'd know which processors each ray depended on, and we'd prioritize based on number of dependencies. There may be a better way to approximate this than a stack. I'll continue thinking.

Charge Exchange

Still progressing. Here's the latest frame:

AAS

Current version of the abstract:

The photoionization-driven evaporation of planetary atmospheres has emerged as a potentially fundamental process for planets on short period orbits. Using AstroBEAR, an AMR multiphysics code, in a co-rotating frame centered on the planet, we model in three dimensions the transfer of ionizing photons into the atmosphere, the subsequent launch of the wind and the wind's large scale evolution subject to tidal and non-inertial forces. We run simulations for planets of 0.263 and 0.07 Jupiter masses and stellar fluxes of 2x10¹³ and 2x10¹⁴ photons/cm²/s. These simulations reveal new, potentially observable planetary wind flow patterns, including the development, in some cases, of an extended neutral tail lagging behind the planet in its orbit. In addition, the role of radiation pressure in shaping exoplanet photoevaporation remains a topic of contention. Radiation pressure from the exoplanet's host star has been proposed as a mechanism to drive the escaping atmosphere into a "cometary" tail and explain the high velocities observed in systems where mass loss is occurring. Using simulations of a planet modeled after HD 209458b and subject to both ionizing and Lyman-α radiation, we demonstrate that for the Lyman-α flux expected for HD 209458, radiation pressure is unlikely to significantly affect photoevaporative winds or to explain the high velocities at which wind material is observed. Charge exchange between the stellar and planetary winds has also been suggested as a method for creating the observed Lyman-α absorption signature. We present results of simulations that explore the effect of charge exchange on our synthetic observations. In addition, we present simulations of the effect of stellar and planetary magnetic fields on our self-consistently launched wind.

MHD

Pressure protections around the planet. Some very small B fields in the planet - may need to fix the field to zero inside the planet at every step?

Also realized that I set the field to zero inside R_p only, not inside R_ob. That may help some. Perhaps it's even worth fixing the field to zero out to within a small but significant distance from R_ob?

Also plan to try fixed-boundary global sim when I get a chance.

Charge Exchange

Still running fairly slowly. I expect it should reach steady state soon, though (assuming we don't enter a cycle with the stellar wind - I am using John's intermediate case).

The wind is turning inward as it approaches the outside boundary - I think this is expected due to Coriolis/centrifugal force/differential rotation? But haven't really had time to investigate yet.

AMR profiling

Actual profiling should finally run successfully on Stampede. Also on my to-do list is to make sure we're not waiting longer than necessary for local/MPI rays.

Radiation Pressure Paper

Resubmitted. Hopefully not too long until we hear back.

Ray tracing

Almost working perfectly. Getting a few unexpected zero-length rays that I need to investigate, and a few strange rays toward the end.

Charge Exchange

Appears to be running well. Trying to make a movie of first six frames.

Magnetic Fields

Still having both pressure and divergence issues in the initialization. May need a brainstorming session for this one. See Mathematica notebook​ with equations to satisfy.

AMR profiling

Still waiting on Stampede.

AAS

Booked hotel for 3rd through 9th. No longer sure where I'm flying out of, but should hopefully know soon.

Radiation Pressure

Submitting Friday. Haven't gotten Ruth's comments yet.

Charge Exchange

See ​last week's blog post for test results.

Test comparison

Production run with John's "moderate" stellar wind is queued on Stampede2. Also queued non-rotating parameter space runs on BlueHive.

Magnetic fields

Added option for initial stellar B field that excludes planet (as if it were a perfect conductor). Getting divergence warning immediately after initialization. I believe the field should be divergence-free, but will this be a problem?

AMR profiling

New test still queued on Stampede.

Ray tracing

Tested patch output first. That's working well. Now just need to figure out best way to read mixed data types in Matlab, and should be able to plot rays without too much issue. Using 8-core, 3-level, 2-D 16x256 base resolution test mesh, which should leave enough room for each ray to be distinct even on the finest level.

Here are the first three steps of charge exchange for a simulation where the densities of stellar and planetary material differ by five orders of magnitude.

 q(iHIIhot) =    55.0433622655540     
 q(iHcold) =    5504281.18319314     
 q(iHIIcold) =   0.000000000000000E+000
 q(iHhot) =   0.000000000000000E+000
 After step:
 q(iHIIhot) =    54.7826014982936     
 q(iHcold) =    5478205.36722786     
 q(iHIIcold) =   0.000000000000000E+000
 q(iHhot) =   0.000000000000000E+000
 Charge exchange coefficient is   3.455999999999997E-003
 Before iHIIhot protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 After iHIIhot protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 Before iHhot protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 After iHhot protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 Before iHIIcold protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 After iHIIcold protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 Before iHcold protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 After iHcold protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 Advanced level  0 to tnext= 0.6000E-08 with dt= 0.6000E-08 CFL= 0.5764E-07 max speed= 0.1921E+02 radiation sub-cycles=       1
 Info allocations    =   148.7 kb  148.7 kb
 message allocations =   ------    ------  
 sweep allocations   =   ------    122.5 kb
 filling fractions   = 
 Current efficiency  =  75%  16%  91% 
 Cell updates/second =        211       209 100%
 Wall Time Remaining =     9.1 day at frame    0.0 of     10
 AMR Speed-Up Factor =       0.2004E+03
 Before step:
 q(iHIIhot) =    54.7752475699149     
 q(iHcold) =    5478092.25506669     
 q(iHIIcold) =   6.222857684820517E-003
 q(iHhot) =   6.222857684820517E-003
 After step:
 q(iHIIhot) =    54.7751651392209     
 q(iHcold) =    5478084.01114326     
 q(iHIIcold) =   6.222848320107427E-003
 q(iHhot) =   6.222848320107427E-003
 Charge exchange coefficient is   3.455999999999997E-003
 Before iHIIhot protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 After iHIIhot protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 Before iHhot protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 After iHhot protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 Before iHIIcold protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 After iHIIcold protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 Before iHcold protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 After iHcold protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 Advanced level  0 to tnext= 0.3124E-04 with dt= 0.3124E-04 CFL= 0.3001E-03 max speed= 0.1921E+02 radiation sub-cycles=       1
 Info allocations    =   148.7 kb  148.7 kb
 message allocations =   ------    ------  
 sweep allocations   =   ------    122.5 kb
 filling fractions   = 
 Current efficiency  =  77%  14%  90% 
 Cell updates/second =        303       302 100%
 Wall Time Remaining =     1.8 min at frame    0.0 of     10
 AMR Speed-Up Factor =       0.2879E+03
 Before step:
 q(iHIIhot) =    22.3818513598937     
 q(iHcold) =    5478051.58190027     
 q(iHIIcold) =    32.3995362683480     
 q(iHhot) =    32.3995362683480     

----

 After step:
 q(iHIIhot) =    22.3818513491939     
 q(iHcold) =    5478051.57928144     
 q(iHIIcold) =    32.3995362528591     
 q(iHhot) =    32.3995362528591     
 Charge exchange coefficient is   3.455999999999997E-003
 Before iHIIhot protection:
 q(iHIIhot) =   -4.07465558395285     
 q(iHcold) =    5478025.12276545     
 q(iHIIcold) =    58.8560431859153     
 q(iHhot) =    58.8560431859153</span>

----

 After iHIIhot protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624     

----

 Before iHhot protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624     
 After iHhot protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624     
 Before iHIIcold protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624     
 After iHIIcold protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624     
 Before iHcold protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624     
 After iHcold protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624
 Advanced level  0 to tnext= 0.9368E-04 with dt= 0.6244E-04...

Jonathan was right on the money with the timestepping being too long. I've broken out the relevant blocks above. For the charge exchange step that follows the first line, we get an overall \Delta of

\Delta = \beta(HIIhot \cdot Hcold - Hhot \cdot HIIcold)dt = 3.456\times10^{-3}(22.382*5478052-32.4^2)6.244\times10^{-5} \approx 26

Since HIIhot is only ~22, this gives us HIIhot of ~-4 after charge exchange is complete.

This run is from after I implemented protections for any species less than 0. The section that follows the first line shows that 4 is added back to HIIhot and Hcold, and removed from HIIcold and Hhot, reversing the excess charge exchange that brought HIIhot negative.

This method of protection can result in large oscillations arount equilibrium. Better would be to prevent any one species from changing more than e.g. 10%, as we do with the other quantities, but this is more computationally expensive thanks to the amount of subcycling that would be required in some places (namely those with large \Delta in comparison to any one species).

AAS

Registered. Have flight info ready, will get hotel once I have flight confirmed.

Abstract: Charge exchange? Should have results by January. Poster, iPoster, or iPoster-Plus?

Radiation Pressure

Ruth said minor comments are on the way. Haven't heard anything from Eric.

Charge Exchange

Last week's problem with no hot neutrals appearing was exactly as expected. That's been fixed. But in the second frame I'm still ending up with stellar material right by the planet, for some reason (which is what the changes that caused the previous issue were prompted by…). I've added some debug statements - if that doesn't work, may be time for another small test problem, though it seems odd that it happens only near the planet - suggests to me that it's a problem with the simulation, or possibly some other multispecies routine (species protections? These should preserve 0s, though), rather than the charge exchange routines. In particular, the fact that there's no problem for an entire frame suggests it's something about the changes caused by the stellar wind.

First frame

Second frame

Freezes on BlueHive when compiled with optimization level 3, before reading in from restart file. Runs fine when compiled with optimization level 2. Should probably figure out why at some point, but for now I'm letting them run.

AMR Profiling

Trying to run on Stampede. Minor problem with profiler, trying suggested workarounds now.

Ray Tracing

Queued on BlueHive. One job should start on Friday; also requested job with shorter time allocation, to see if it might start earlier (don't need long).

AAS

Need to register by 9/25, abstracts by 10/8 (if we want to submit something). I'm not currently an AAS member - not sure what the turn-around time is to become one, but it's a significant difference in registration fee.

Flights are ~$1200, similar for (AAS) hotel, though looks like about 2/3 price is available.

Radiation Pressure

Sent around what should be the nearly final version of the ​radiation pressure paper along with the response to the referee. I'll plan to submit as soon as I have a response from everyone (hopefully end of the week or beginning of next).

Charge Exchange

I've removed everything that I believe can eliminate the added hot neutrals from charge exchange (ionizing radiation, reverse charge exchange, collisional ionization [not turned on originally]), and checked that the charge exchange rate is being calculated correctly in the actual line transfer step (rather than in diagnostics). We have rates up to order 10 in some places, so with a dt of 0.1 over a frame, we should end up with some significant amount of hot neutrals after a frame now. Still getting none, however. Feel like I have to be missing something stupid now.

AMR profiling

Stampede has the Intel profiler available. Haven't had a chance to use it yet, but as soon as the charge exchange is working I'll check it out.

Non-rotating charge exchange

Feel like this should wait until we think we have charge exchange working, unless there's some reason it might give a clue to what's going wrong? I can't think of any, though. Ready and waiting, though.

Radiation Pressure Paper

Believe almost every comment has been addressed. Two things remaining:

• Neutral speed plot: better version below?

• Timescale for bulging of wind. It may actually be a box effect (cycle is very close to 2 crossing times), though I'm not yet sure physically how reaching the edge of the box might increase the effect of radiation pressure.

Charge Exchange

See ​my previous blog post for the latest thoughts on charge exchange. It's recently started again on Stampede and is moving fairly quickly (~6 days left as of this morning). Here are the latest frames:

Species proportions

Charge exchange potential

Non-rotating Charge Exchange

Using the terminal velocity and mass loss rate from ​this thesis to approximate the wind of an O6 star at the orbit of Jupiter, I get a wind density of ~2x10^-18 g/cm³, which is within an order of magnitude of our current values of 4.5-6.5x10^-19 g/cm³. Given this, I've queued the rotating tests with the same wind as HD209458b (though technically they should be corrected for a different planetary and stellar mass).

AMR line transfer

Mixed news on this front. Bad news: The (current on the repo) version of the AMR line transfer is taking about double the amount of time on BlueHive as the previous. Good news: Based on one frame, it appears to be working for a full problem (rather than the simplified tests).

I'd like to run this in a profiler, but I haven't yet gotten it working (getting errors when attempting to run through slurm). I may try a run on an interactive session on a single node, though that will be far fewer cores, which may affect the hotspots…

For reference, here's the updated table:

Charge Exchange

A few more frames. So far, charge exchange is occurring, presumably creating a hot neutral population which is then immediately ionized (on the next line transfer step [not subcycle], by both photoionization and possibly ionization by "reverse" charge exchange). This suggests that the amount of neutral stellar material should be alternating between zero and some small value every step, which means at some point a frame may be output with nonzero amounts of Hhot….

Charge exchange should be creating at most ~10^-4 Hhot (in CN). The highest-rate locations are in the shadow of the planet. HIIcold reacts with Hhot, and we have large amounts (~10^-2-10⁰) of HIIcold. This reverse charge exchange is already taken into account in the rate calculation, but if it's oscillating we may not have seen the opposite side of the step yet.

All this to say I'm not yet sure where the Hhot is going, and I haven't figured out the appropriate diagnostics to look at to make sure everything is happening as it should be, even if it's not quite the expected result (everything being zero when Hhot is zero, and not having seen a frame with nonzero Hhot yet).

Here are some figures for the latest frame:

rho, 1-X, T

Proportion of species

Clockwise from top left: Ionized cold, ionized hot, neutral hot, neutral cold

Amount of species

"Charge exchange potential"

min(Hcold, HIIhot)/max(Hcold,HIIhot) - so proportion of material that can undergo an exchange, disregarding the reverse interaction.

Charge exchange rate (CU)

Multiply by dt to get delta Hhot

AMR line transfer ==

Still hasn't started on BlueHive.

Non-rotating charge exchange

Still waiting on AMR line transfer.

Paper

Good discussion Thursday. Haven't had much chance to start making figures or new changes yet.

Charge Exchange

Working on stampede. Here's the first frame:

I believe the neutral planetary material top-left is from the higher density (higher recombination rate, possibly also from higher temperature) in the stellar-planetary wind shock. Note that no charge exchange has occurred yet (but there's little contact between neutral planetary and ionized stellar material so far).

Trying to make a VisIt expression for charge exchange potential, but the following gives a "cannot divide by zero" error:

if(gt(Hcold,1e-2), HIIhot/Hcold, 0)

Note that the division should never happen if there could be a zero in the denominator, so I'm not sure why VisIt is complaining here…

AMR line transfer

Fixed the missing condition on the displaced child mask update. Still waiting on newest test on BlueHive.

Non-rotating charge exchange timing tests

Have the steady states ready. Need to calculate the stellar winds - use original distance (vs. larger orbital distance e.g. O-type)?

Radiation Pressure Paper

Got referee response. Sent around preliminary thoughts, will try to start addressing them this week. To recap my email:

• Collisional nature of the wind (not hard-body collisional but electrostatically collisional)
• Clearer distinction between our focus on the cometary tail vs. the observed velocities
• Discussing more thoroughly the absorption of momentum by ionized hydrogen

Charge Exchange

As usual, Stampede is continually frustrating (related: Jonathan, you mentioned a plan to migrate from gitlab.circ to git.circ—any ETA on that?). Still making progress, though.

AMR timing

Baseline test complete (single frame of HD209458b steady state). What we believe to be the most efficient version is queued on BlueHive. Here are the relevant stats from the baseline run:

Other tests

• Charge exchange timing for non-rotating parameter space planets
• AMR ray tracking

Misc

Posted my suggested progress through Toro on a ​blog post. Feel like there should be a better place to put this - perhaps an "Introduction for New Students"-type page would be useful?

• Read Chs 1-3 on the Euler equations and solutions to hyperbolic equations in general. Some of this will likely be familiar.
• Read Ch 4 and write the exact Euler solver described.
• Read Ch 5 and write the Godunov upwind scheme, for the linear systems described.
• Read Ch 6 and update the Godunov upwind scheme to use the exact Euler solver.
• Skip 7&8
• Read Ch 9 and write the approximate-state solvers (good place to find out if you're following good programming practices).
• Read/skim Chs 10-12 and write at least the HLLC solver (which is the one used in AstroBEAR). Note that there have been some significant improvements to the HLLC chapter in the third edition. Write the rest if you're interested, though I wouldn't recommend the Osher solver (it's interesting mathematically, but doesn't seem to offer much over the other solvers).
• Read Ch 16 and write a 2D/3D solver (test problem in Ch 17). Another place where you'll really notice if you've been following good practice.
• For second-order solvers: Read Ch 14 (skim 13 for an introduction to the same concepts for scalar equations) and write the WAF and MUSCL schemes (esp. MUSCL-Hancock, used in AstroBEAR).

Note that the ends of the code-related chapters have example code. I'd strongly suggest looking at some modern Fortran style guides rather than following the old style given in these, not least because it's much easier to follow your own code if it's written with meaning.

Charge Exchange

Modified to track stellar and planetary ionization rates separately. I believe this is the most efficient solution to dividing the ionization and recombination properly among the species, since a single rate can't account for the difference in proportion between the recombination and ionization rates.

rates(i,j,k,iPlanetIoniz) = f*(1d0-emtau)/dx*(q(iHcold)/q(iH))
rates(i,j,k,iStellarIoniz) = f*(1d0-emtau)/dx*(q(iHhot)/q(iH))

rates(i,j,k,iPlanetIoniz) = rates(i,j,k,iPlanetIoniz) - RecombinationRate(T, ne*(1-ambFrac), q(iHIIcold)*(1-ambFrac))
rates(i,j,k,iStellarIoniz) = rates(i,j,k,iStellarIoniz) - RecombinationRate(T, ne*(1-ambFrac), q(iHIIhot)*(1-ambFrac))

Non-rotating timing on BlueHive

Need to get the data from Baowei. Then will do this after AMR tests are done.

AMR line transfer

The first timing test (projected-grid, for comparison) is queued. I'm using a single frame (220→221) of the steady state of HD209458b, no radiation pressure.

Tracking rays

Lowest testing priority right now, but ready once other tests are complete. Contemplating how best to implement plotter for output in Mathematica/Matlab/python.

Misc

Working through pre-req to fall course. ​My website is up.

Radiation Pressure Paper

Updated version is on the arXiv today. I've emailed it out to a few people.

Charge Exchange

Waiting for job to start on BlueHive to get diagnostics.

AMR line transfer

Pushed two commits with separated local and MPI rays, and keeping an integer pointer to the next free ray and ray to be processed (instead of looping through entire array each time). Not tested yet, but weren't complicated changes, so (fingers crossed) they should work immediately. Now just need to compare timings to see if it was worth it for plane-parallel.

Also pushed commit with output for following ray propagation across the grid in time order (synchronized time, ray initial position, and processor ID). This may require some testing, particularly to determine the best way to process it for visualization.

Misc.

Set up ​GitHub repository of my Euler code. Working on personal website, as well.

Radiation Pressure Paper

Submitted to MNRAS and the arXiv. Also looked at renormalized absorption:

Will update paper on Overleaf soon.

Charge exchange

Ionized and neutral planetary material appear to be interacting and becoming stellar material. Think the problem must be here:

if (q(iHhot) /= 0) then
  prop(i,j,k,lt_iHhot) = q(iHhot)/(q(iHcold) + q(iHhot))
else
  prop(i,j,k,lt_iHhot) = 0d0
end if

if (q(iHcold) /= 0) then
  prop(i,j,k,lt_iHcold) = q(iHcold)/(q(iHcold) + q(iHhot))
else
  prop(i,j,k,lt_iHcold) = 0d0
end if

if (q(iHIIhot) /= 0) then
  prop(i,j,k,lt_iHIIhot) = q(iHIIhot)/(q(iHIIcold) + q(iHIIhot))
else
  prop(i,j,k,lt_iHIIhot) = 0d0
end if

if (q(iHIIcold) /= 0) then
  prop(i,j,k,lt_iHIIcold) = q(iHIIcold)/(q(iHIIcold) + q(iHIIhot))
else
  prop(i,j,k,lt_iHIIcold) = 0d0
end if

data(:,:,:,lt_iHhot)=data(:,:,:,lt_iHhot)-rates(:,:,:,iIoniz)*dt*prop(:,:,:,lt_iHhot) + rates(:,:,:,iChargeExchange)*dt
data(:,:,:,lt_iHcold)=data(:,:,:,lt_iHcold)-rates(:,:,:,iIoniz)*dt*prop(:,:,:,lt_iHcold) - rates(:,:,:,iChargeExchange)*dt
data(:,:,:,lt_iHIIhot)=data(:,:,:,lt_iHIIhot)+rates(:,:,:,iIoniz)*dt*prop(:,:,:,lt_iHIIhot) - rates(:,:,:,iChargeExchange)*dt
data(:,:,:,lt_iHIIcold)=data(:,:,:,lt_iHIIcold)+rates(:,:,:,iIoniz)*dt*prop(:,:,:,lt_iHIIcold) + rates(:,:,:,iChargeExchange)*dt

rates(i,j,k,iChargeExchange) = chargeExchangeRate*(data(i,j,k,lt_iHIIhot)*data(i,j,k,lt_iHcold) - data(i,j,k,lt_iHhot)*data(i,j,k,lt_iHIIcold))

Charge exchange rate should be zero whenever there is no planetary or stellar material. Proportion of stellar material should be zero when there is no stellar material. So it should be impossible for q(iHIIhot/cold) to become nonzero when they're both zero at the start of the line transfer step. Clearly they are, however - will put some diagnostic variables in to see what exactly is nonzero when it shouldn't be.

Also need to check timing for charge exchange runs on non-rotating case of parameter space runs on BlueHive (will need to get the last frames).

AMR line transfer

Will try to run timing this week (after charge exchange diagnostics are complete), in addition to working on a new commit with optimizations.

Upcoming/Possible Projects

In approximate order of importance (though not necessarily priority).

• Literature review
• Spherical line transfer
• Metastable helium line
• Hot Neptune
• Flares

Fall 2019 registration

Last semester to take courses for credit. Would like to take CSC 482 (designing efficient algorithms), I think - will talk with Laura as well.

MHD proposal figs

Also, a test of the streamlines. I like the LICs better, since the streamlines have a habit of not ending.

​Radiation Pressure Paper

Found issue with postprocessing. Reprocessing synthetic observations. Here's a first look at the difference. The difference is significant, but doesn't change our conclusions.

Charge Exchange

Current frame

We also need to decide on the high-med-low for the stellar wind. Here are the current parameters (taken from John's simulation, except the mass loss rate, which is the solar mass loss rate - which corresponds to a 4x higher wind density than their strongest wind):

AMR line transfer

Timing run is next up (for comparison with non-AMR and new optimizations), along with work on optimizations.

Also working on a way to output propagation of rays across grid for a step (just for a nice visualization):

Need: ray position, synchronized time

Ray position is easy.

Synchronized time pseudocode:

call MPI_barrier     ! To synchronize time - unless MPI_wtime_is_global, in which case all processors are synchronized already and this becomes much easier
initTime = MPI_wtime
call PropagateRay

subroutine PropagateRay

time = MPI_wtime - initTime
MPI_file_write(time,position)    ! Or might be easier to save an array of time, position and write it out after all rays are processed, in processor order (so we don't have to calculate the buffer)

​Radiation Pressure Paper

All changes received so far made. I'm ready to submit.

Also, we still need to upgrade the Overleaf account to enable sharing again.

Charge Exchange

Running successfully on Stampede. On 20 nodes (960 cores), predicted time to completion is currently 1.8 months (got 2 frames in 48 hours). I expect this will speed up significantly once the stellar wind is more completely introduced and we reach closer to steady state.

Cost to allocation would be 1.8 mo x 30 days/mo x 24 hr/day x 20 nodes = 26000 node-hours. Much higher than I'd like.

Current frame

We also need to decide on the high-med-low for the stellar wind. Here are the current parameters (taken from John's simulation, except the mass loss rate, which is the solar mass loss rate - which corresponds to a 4x higher wind density than their strongest wind):

AMR line transfer

Final functional fixes made (except a couple relating to subcycling). Some optimizations remaining, but could start using it for runs. Still want to do timing comparisons between old and new routines, and would be good to get a baseline for comparison of various optimizations.

Radiation Pressure Paper

Moving this to the top of my to-do list. Should have first final look done in next couple of days.

Charge Exchange

A few more Stampede problems. Made it up to line transfer step with last run, so I'm hopeful it will just run now.

AMR line transfer

First clump test (radiation pressure only) still partially successful (See ​code page, though I have a few more frames). Local ray communication works, and fixed race condition (though still may switch to buffered sends). Currently tracking down source of unmatched sends and/or receives.

Charge Exchange

Square appears to be from a combination of diffusion and advection (advective ionization is significantly negative). Chombo conversion appears essentially correct (slight discrepancies between max and min of rho?), so attempting to run on stampede.

AMR line transfer

First clump test (radiation pressure only) partially successful (See ​code page. I'm running into what appears to be a race condition in the MPI library - multiple recvs are being used to accept the same sends. I've put debugging this on hold while developing the local (non-MPI) communication.

Currently testing the MPI-less implementation.

AMR Line Transfer Tests

Running the same test seen ​here, which is IonizationFront module. Parameters are given in Mathematica notebook attached to that page, but for thoroughness:

Dependent on test (initial test parameters listed)

Working my way up the test hierarchy.

* Single ray (1D), single processor. Ran for very long simulation time, to give solid point of comparison.

* Multiple rays (2D, GmX(2) = 32, GxBounds(5) = 2.5d-1), single processor (final_time = 1d1, final_frame = 100), refined from x = 0.25 to 0.75.

* Multiple rays (2D), 2 processors (neighbor patch).

* Single ray (1D), single processor, 2 levels in middle-ish (child patch, MaxLevel = 1).

Result of refinement for optical depth

* Multiple rays (2D), single processor, 2 levels (child patch).

* Multiple rays (2D), multiple processors, 2 levels (neighbor(? - may need to force somehow) and child patches).

* Multiple rays (2D), multiple processors, 2 levels (neighbor and child patches), with domain boundary refined.

* Single ray (1D), single processor, for radiation pressure.

* Single ray (1D), single processor, for radiation pressure.

* Compare runs from ​combination test.

Comments: y and z directions are treated identically, so it's unnecessary to test both "2D" and "3D" cases (note that nDim is always set to 3 in global.data, resolution just set to 1 and size very small in "1D" and "2D" cases). Believe the above set should test all cases in each if statement.

Some working changes have been pushed.

AMR line transfer

Jonathan and I worked out the remaining logic last week. Short summary of what we discussed ​here. Also includes some thoughts on spherical line transfer.

Current to-do list is on AMR_line_transfer_2 branch in linetransfer_control.f90. Will need to clean up git when it's ready.

Other projects

• Radiation pressure
  • Needs comments, punching-up of how we address Cherenkov paper
• Charge exchange
  • On debugging to-do list, for weird recombination(?) in bottom right hand corner

Radiation Pressure

Final measurements made. Here are related plots:

During bubble-out portion

During blowback portion

Charge Exchange

Compiles, runs, and does something. Need to test it — possibly determine steady-state mixing ratio for some density? Also still seeing that excess neutral square behind the ablated planet material. It almost looks like the recombination shadow we see in the other simulations, but since we don't have any ionizing radiation I don't see what would maintain the completely ionized state outside that region (perhaps it's not actually complete?). Could also be from some mixing, or more time spent there by the stellar wind… thoughts?

Uses HD209458b parameters, as given in ​this post.

Initial state

Final state

AMR line transfer

Made good progress last Wednesday. Latest commit just has a short to-do list, mostly cleanup. Primary functional code left is implementing an MPI derived type for sends and receives of rays, and a custom scatter and gather (because of differing array sizes).

Radiation Pressure Paper

Believe it's nearly done. Could use some help with the results comparison section and the punchiness of the conclusions, and there's one more measurement I need to make (will likely have to wait until after my qual).

Next Steps

• Finish AMR line transfer
  • Believe I have a working algorithm, just need to finish implementing it
• Charge exchange
  • Will run just stellar wind, no ionizing flux, to see effects of mixing only
  • May be able to run later this week
  • Think it would be useful to run HD209458b with just stellar wind as well as charge exchange, both for comparison with Cherenkov et al. and to separate the effects of each (under the assumption that John's results may not be representative of the results for the significantly different planet)
• What next? Spherical line transfer?

Qual

• What should I be reviewing?
• What should I focus on in the brief (10 journal-style pages)?
  • Should I try to work in all of the papers, tell the story so far in limited detail? Focus on just one?

​Radiation pressure paper

• Sent an email around with information related to our previous discussion. I've begun incorporating some of the information where appropriate. Body of email copied below for reference.

• Could use some brainstorming on alternate ways to make the same point as these measurements:
  • Comparing the average speed of a particle being significantly affected by radiation pressure to the escape speed of the system
  • F_{L-\alpha}, F_{th}, F_{ram} comparison (could also compare pressures, if I can figure out a good length scale to convert volumetric radiation force to radiation pressure)

Even using the high-memory nodes on the visual partition, VisIt either crashes (at end of trying to load Chombo) or hangs (while processing) when trying to make the necessary plot. Possible to load all of the VTKs created when saving 3D plot in python, maybe?

• Cherenkov et al 2017 ionization procedure:
  • This is the entirety of what they say about their ionization procedure: “In the grey approximation, the last term in equation (6), the rate of photoionization, can be written in the form , where h is the Planck constant, σUV = 6.3 × 10−18 cm2 is the photoionization cross-section at the Lyman limit frequency ν0, FUV = F0 (apl2/l2) is the intensity of ionizing radiation at distance l from the star centre, F0 = 884 erg cm−2 s−1 is the radiation intensity at the planet orbital distance, apl, (Schneiter et al. 2016) and τ is the optical depth.”
  • Sections 2.3 and 2.4 talk about the radiation transfer model they use, but only in the context of the radiation pressure (Lyman-alpha line). My assumption would be that they intend for this to extend to the ionizing radiation, but they never make this explicit.
  • None of the previous studies that could reasonably be considered to be by their group use self-consistent ionization
• Launch conditions from Bourrier et al.
  • Bourrier et al 2013 launch metaparticles of neutral hydrogen only, and stop tracking ionized hydrogen—this may have a significant effect on the dynamics, since much of what is being pushed on (indirectly) is ionized
  • Their highest density is at launch radius, and is ~107 /cm3. Around the same location (Mach surface/Hill sphere at ~2.8 Rp), we have a total density of ~6x106 /cm3, with a neutral density of ~8.5x104 /cm3 (for no-flux simulation). So their optical depth to Lyman-alpha radiation is significantly higher.
  • They say “the number of atoms in a meta-particle is calculated to keep a low value of the corresponding optical depth dτ in front of a fraction of the occulted stellar surface.” Unclear to me what exactly this means.
  • Their metaparticles are launched randomly in every direction with thermal speed of 11000 K (much higher than our temperature at this point).
• HD 209458 vs. other stars' Lyman-alpha fluxes
  • From Wood 2005, Lyman-alpha fluxes for other G0 V stars normalized to HD209458b orbital distance, in erg/s/cm2: 18344, 18955, 14374, 23836; HD209458 = 6409
  • So it’s on the low side, by a factor of ~3
• Optical depth in wings (near planet)
  • At 3 Rp in +y, 50 km/s ± 2km/s bin, optical depth of 1.66 (in high flux case)
• What happens as we transition between ionizing flux regimes?
  • From Sec 5.3 of John's paper: While in E-limited, neutral fraction constant, density higher, so total neutral density (and therefore optical depth) higher. Crossing into R-limited, neutral fraction reduces as 1/F1/2, and density increases at the inverse rate, so optical depth becomes constant. At some point, wind may start turning around to create the swirly arms – at this point, optical depth should ~double from previous.
• Local ionization equilibrium
  • Essentially maintained at 0.03 events/sec throughout wind, outside of small region - see attached VisIt plots (​high side, ​high top, ​med side, ​med top, which plot the absolute value of RecombinationRate-IonizationRate (so they're a measure of the distance from ionization equilibrium)
• Effect of averaging over bubbling in high-flux case
  • Very small effect, mostly deeper center transit in the phase where the wind is expanding out of the Roche lobe. See ​attenuation and ​observation comparisons.
• Plots of velocity, neutral density
  • ​Final attached plot is neutral mass density with velocity contours plotted over. Contours are: blue, 1x10^{5 cm/s; red, 5x10}5 cm/s; cyan, 1x10^{6 cm/s; yellow, 5x10}6 cm/s; magenta (none present), 1x10^{7 cm/s}

AMR line transfer

Making good progress. I believe I have the implementation almost completely defined, just remains to write some of the organizational code. The algorithm I have isn't as efficient as I'd like, but I haven't yet figured out a way around the limitations. Once I have a working prototype (end of the week?) it may be worth reviewing to see if there are any speedups or reduced memory usage possible.

In some searching, I also found what appears to be a ​decent spherical algorithm for ray tracing. It's an approximation method combining the long-characteristic and short-characteristic methods that was implemented for FLASH. I definitely think it's something we should consider for AstroBEAR (I haven't tried yet, of course, but I believe it's possible to implement it in our code).

Radiation pressure paper

Have looked at Luca's comments. The biggest is ​a paper by Cherenkov et al., who have done something similar but included the Doppler shift in the absorption coefficient. I can't tell how they calculate the optical depth for ionizing radiation. They also include the stellar wind and start with a much higher-temperature (~2x) and higher-density (~100x) planet. Their resolution is significantly smaller, though it changes in a way that's unclear (I believe they just use unevenly-spaced points) from the outer portion to the inner portion of the grid.

There are also a couple of comments on clarity that I could use another set of eyes on.

Finally, the measurements that I'm attempting to make keep crashing VisIt. I'm going to try yet another method this week.

Project Statuses

• Radiation Pressure
  • ​Paper sent around to everyone (astrobear Overleaf account still on free plan)
  • Luca's looked at it and made some small edits and some comments — haven't looked at any of them yet
  • Two measurements still to make (I've commented at the location)
    • Radiation pressure vs ram pressure vs thermal pressure at edge of wind (need to think about this one some more)
    • Average speed of material absorbing significant amounts of flux (currently using between 100% and 1% of maximum radiation pressure as boundary of "significant")
• AMR line transfer
  • Sketch of code written, significant modifications required
  • Main focus for next couple of weeks
• Charge exchange
  • Will run just stellar wind, no ionizing flux, to see effects of just mixing (can run this week)
• HD209458b
  • Need to discuss at some point whether we want to do a paper on just the planetary wind
  • High ionizing flux - effect on mass loss rate/wind structure

Charge Exchange

​Radiation Pressure Paper

Figures pretty much done. Need to put the medium flux and high flux average results side by side, and most could use some minor tweaking for presentation. Currently only using the streamline image for the high flux case, but may want to add at least radiation pressure (especially if it turns out that a lot of the absorption is happening in a fast-moving shell).

High flux results are as expected on average, but the results for the blowback phase are worrisome - could use a second opinion.

Bubbleout (​245-290)

Full average (​245-315)

Blowback (​290-315)

Will have a first draft shared with everyone by end of day Friday.

Radiation Pressure

Another 3 frames done. I have it running still, but if this is sufficient, I can stop it and start postprocessing so we can get the remaining plots for the paper.

Parameter space

Spent a good portion of last week retrieving data and making movies. Checking the proofs and sending back today - should be completely done by tomorrow evening.

Parameter Space Paper

Submitted. Need to get password for one of macs to get data for movies to upload to YouTube. Also need to make list of people to send to.

Radiation Pressure

7 additional frames finished. VisIt is having problems drawing the Mach contour the first time (works on redraw) — any thoughts?

There may be a slight bulge starting at the planet again, but it's still too early to say definitively.

AMR line transfer

Jonathan and I had a productive meeting two Wednesdays ago. I should be able to sit down and get a functional version this week.

Project Statuses

• Radiation Pressure
  • Outline of ​paper
  • Need to make a few more plots for high Lyman flux case, once it's completed (early December?)
  • Will address Eric's email, as well - measure speeds of blowoff at a variety of points along side and compare to energy requirements
• AMR line transfer
  • Sketch of code written, significant modifications required
• Charge exchange
  • Next will be stellar wind, no ionizing flux, to see effects of just mixing
• HD209458b
  • High ionizing flux - effect on mass loss rate/wind structure

Project Statuses

• ​Parameter space paper
  • Ready to resubmit
  • Need to retrieve data for some movies, then to be posted on YouTube
  • Waiting on John to submit to arXiv
• Radiation Pressure
  • Outline of ​paper
  • Need to make plots as listed
  • Need a few more frames of high flux - keep running out of memory, currently trying 8 GB/CPU
• HD209458b
  • High ionizing flux got 7 frames in 5 days, so could probably see changes in ~1 month
• Charge exchange
  • Believe I have stellar wind on the grid successfully
  • Need to test charge exchange now
• AMR line transfer
  • Sketch of code written (see ​previous blog post)
  • Jonathan and I are going to discuss this week sometime

HD209458b

Here's frame 7 of the high ionizing flux run. You can see the hotter wind material bubbling out into the Hill sphere.

Charge Exchange

Stellar wind appears to be successfully on the grid.

Initial:

Final (link to movie):

​Parameter Space Paper

Just a few comments left. Copied from my email:

Places/comments that could still use some attention:

• General qualitativeness and informality (anything in particular stand out?)
• Where to discuss the appropriateness of the case B recombination assumption - in methods where we state the assumption, or later in the results/conclusions somewhere?
• I've added a few citations to John's paper. Further suggestions for places to add a reference, or more discussion (particularly in section 4.2), are welcome.

Response-only:

• How to respond to single-frequency suggestion (Just say we've considered it?)
• Square features in streak images from lack of velocity to convolve with noise

In addition to reviewer comments, I've updated the figures with new colorbar ranges that should highlight some of the features more clearly.

HD209458b synthetic observations

Still processing. Should have results for low and medium radiation pressure soon. High radiation pressure run is complete. Average observations over burp, non-burp times? I'll look to see what John did, as well.

HD209458b mass loss rate

High ionizing flux crashed with MPI_WAITANY error, at level 4 after a few time steps. First time this has happened (code is unchanged), so I've requeued it, but I'll also go see if I can find the cause. Unfortunately, the error message doesn't give a location of the call. Fortunately, we only have a few calls to MPI_WAITANY.

AMR line transfer

Have a maybe halfway there commit on new branch AMR_line_transfer. Here's the relevant code:

linetransfer_declarations

Ray type and linked list of rays. Ray may need a few more data structures inside.

  type RayDef
    integer :: sendRequest, recvRequest
    real(kind=qPrec) :: currentIonizingTau, currentIonizingFlux, currentLymanTau, currentLymanFlux
    type(InfoDef), pointer :: info
    integer :: localCoord(nDim), level
    type(RayDef), pointer :: next
  end type

  type RayList
    sequence
    type(RayDef), pointer :: self
    type(RayList), pointer :: next
  end type

linetransfer_control

Truncated code for integration along ray, only change is we don't have to hack q from the layout data any more.

    type(infodef), pointer :: info
    type(RayList) :: rayList
    type(RayDef), pointer :: ray

    !!! Still only want to do the calculation once, after each hydro update is complete
    IF (n < FinestLevel) then 
      linetransfer_iters = 0
      SubcyclingReason = 'N'
      RETURN
    end if

    call StartTimer(iLineTransferTimer,n)

    ! Convert total to internal energy on all levels
    CALL ConvertLevelsTotalToInternalEnergy

    dt=levels(n)%dt
    tnow=levels(n)%tnow
    tfinal=tnow+dt

    linetransfer_iters = 0
    DO WHILE (tnow < tfinal)

      rates = 0d0
      ! Set info to default values
      if (lt_iF /= 0) info%q(:,:,:,iFlux) = 0d0
      if (lt_iD /= 0) info%q(:,:,:,iDepth) = 230d0
      if (lt_iFa /= 0) info%q(:,:,:,iFlux_a) = 0d0
      if (lt_iDa /= 0) info%q(:,:,:,iDepth_a) = 230d0

      !!!!!!! Things to think about with sends/receives:
      !!!! Should tag each ray individually (related to position on left boundary; corrections for split and joined rays?)
      !!!! When crossing child patch, can't continue on parent until we've received the ray back from child
      !!!! How and when can we collect the rays from finer grids?
      !!!! When can we distribute the rays to coarser grids?

      !!! Nodelist is local, right? And if a processor has a parent and child, can it communicate the ray without MPI?
      ray=>rayList%self   ! Head of linked list for this processor
      DO n=0,MaxLevel
        nodelist=>Nodes(n)%p
        DO WHILE (ASSOCIATED(nodelist))
          info=>nodelist%self%info    ! info points to current patch info
          ! Loop over left boundary of current patch
          do j = 1,info%mx(2,2)
            do k = 1,info%mx(3,2)
              if (n == 0 .and. left boundary == xmin) then    ! If we're on domain boundary at root level, set ray to boundary conditions
                ray%level = n
                ray%currentIonizingTau = 0
                ray%currentIonizingFlux = IonizingFlux
                ray%info=>nodelist%self%info    !!! ray%info now points permanently to the correct patch info?
                ray%localCoord = (/1,j,k/)    ! zone index for local patch
              else
                if (info%childmask(0,j,k) == -1) then   ! If we're a finer grid, receive ray from parent and divide
                  call mpi_irecv(receivedRay, 1, mpi_defined_type, parent, linetransfer_tag, mpi_comm_world, recvRequests(parent))   ! Receive ray from parent   !!! Define a ray type for mpi communication? Only need tau and flux; where is processor of parent stored?

                  !!! but this can only actually be done once ray is received... so move it to processing loop, and need to process rays based on origin
                  rayNum = 1
                  ! Divide ray from parent across appropriate cells
                  do while (rayNum <= 4)
                    ray%level = n
                    ray%currentIonizingTau = receivedRay%currentIonizingTau
                    ray%currentIonizingFlux = receivedRay%currentIonizingFlux/4d0
                    ray%info=>nodelist%self%info
                    ray%localCoord = (/1/)    !!! not sure how to set this one yet
                    if (rayNum < 4) rayList=>rayList%next
                    rayNum = rayNum + 1
                  end do
                else if (info%childmask(0,j,k) > 0 .or. info%childmask(0,j,k) == neighborchild) then        ! if we're on a coarser grid not at domain boundary or patch boundary, do nothing   !!! How is neighborchild set? It's not a parameter...
                else
                  call mpi_irecv(ray, 1, mpi_defined_type, neighbor, linetransfer_tag, mpi_comm_world, recvRequests(neighbor))   ! Receive ray from neighbor   !!! Define a ray type for mpi communication? Also probably need an array of requests, rather than storing them in the individual rays, so we don't have to iterate through the list to check them all...
                  ray%level = n
                  ray%info=>nodelist%self%info
                  ray%localCoord = (/1/)    !!! not sure how to set this one yet
                end if
              end if
              rayList=>rayList%next
            end do
          end do
          nodelist=>nodelist%next
        end do
      end do

      !!! This will return to head of list, right?
      ray=>rayList%self
      do while (associated(rayList))
        if (ray%currentIonizingFlux /= IonizingFlux) call mpi_waitany(recvCount, recvRequests, indx, mpi_status_ignore)   !!! wait for new ray if we're not on the left boundary of domain
        info=>ray%info
        do i = 1,info%mx(1,2)
          j = ray%localCoord(2)
          k = ray%localCoord(3)
          if (i == info%mx(1,2)) then   ! If we're at right boundary of current patch
            if (info%childmask(i,j,k) == 0 .or. info%childmask(i,j,k) == neighborchild) then
              call mpi_isend(ray, 1, mpi_defined_type, neighbor, linetransfer_tag, mpi_comm_world, sendRequests(neighbor))    ! Send to neighbor
            else
              !!! Need to sum appropriate rays here - mpi_reduce? Except collectives are blocking, so we'd never get to another ray to do the sum... so need a conditional to test when we're ready to send a ray to parent
              call mpi_isend(ray, 1, mpi_defined_type, parent, linetransfer_tag, mpi_comm_world, sendRequests(parent))    ! Send to parent
            end if
          else if (info%childmask(i,j,k) > 0) then
            call mpi_isend(ray, 1, mpi_defined_type, child, linetransfer_tag, mpi_comm_world, sendRequests(child))      ! Send to child
            do i = i,info%mx(1,2)   !!! I think this will do what I intend... but not certain yet
              if (info%childmask(i,j,k) < 0 .and. .not. neighborchild) then
                call mpi_irecv(ray, 1, mpi_defined_type, child, linetransfer_tag, mpi_comm_world, sendRequests(child))     ! Receive from child
                exit      ! Return to outer loop at new location      !!! Need to wait to receive summed rays from child before we can continue along this ray...
              end if
            end do
          else

            !!! can break this back into subroutine(s)

            ! Goal is to take input flux and integrate equations across grid to find ionization rate and heating rate
            ! Flux has units of number per area per time in computational units
            ! Data has computational units of the density of neutral hydrogen, ionized hydrogen, ionization rate, and heating rate
            q=>info%q(i,j,k,:)
            ! Calculate radiative transfer of ionizing flux
            .....
            ! Calculate radiative transfer of Lyman-alpha flux
            .....
            !!! Could do all of this as array operations instead of zone by zone - probably more efficient?
            ! Adjust heating and cooling rates if they exist
            .....
              ! Recombination
              .....
              ! Recombination cooling
              .....
              ! Lyman alpha cooling
              .....
            ! Charge exchange
            .....
            !Update minimum time step by looking at rate of energy, neutral, and ionized species
            .....
          end if
        end if
      end do
      rayList=>rayList%next
    end do

#if defined _MPI
    CALL MPI_ALLREDUCE(MPI_IN_PLACE, dt, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, iErr)
#endif     

    dt=min(dt,tfinal-tnow)
    IF (dt == 0d0) THEN
      WRITE(*,*) 'dt = ', dt, tnow, tfinal
      STOP
    END IF

    linetransfer_iters=linetransfer_iters + 1
    tnow=tnow+dt

    call mpi_waitall(sendCount, sendRequests, mpi_statuses_ignore)   !!! wait for all sends after we've finished the rays that we needed to receive
        
  END DO

Parameter space paper

Spent most of last week working on this. Fixed the issue where the colormap wasn't clipping at the lowest value (interpolation gives NaNs beyond range). Need to figure out the best new range still, but here's a preview:

Also made a number of edits, and a few places still to work on. Would like to have it resubmitted by end of next week (proposals are all due this week, I believe?).

HD209458b

Synthetic observations aren't promising here. We're using up to about half the flux reported in the literature, which would probably up the mass loss rate by a factor of ~1.9 (to about 10^{10, still a little on the low side).}

Radiation Pressure

High Flux

We get the same sort of burping John sees in his medium-wind case.

Med Flux (Actual HD209458b Lyman-alpha Flux)

Very subtle widening of the arms, due to some material being pushed radially outward.

Low Flux

Even more subtle widening, very similar to medium flux.

Medium and low flux appear to be pretty much steady-state. Probably don't need to get more frames - doesn't look like there'll be much change to average out. All three still queued on Stampede, but we're pretty low on time now. Satisfactory to make synthetic observations?

Radiation Pressure

Runs on Stampede2:

Running on 128 nodes (6144 cores), we currently (frame 211) have an efficiency of 13% hydro, 38% line transfer (including layout communications), 52% total. Estimated completion 2.9 days. This is down from 23%, 54%, 76% at frame 203 with estimated completion in 1.5 days, where the current run started. Reduced estimated wall time by a factor of ~10. Produced 8 frames in the last 6 hours.

Max speeds are only on the order of 10s, and haven't changed significantly in current run (nor has dt). Here are some comparisons of the mesh:

Frame 180 (pre-radiation pressure):

Frame 203:

Frame 211:

So it's definitely related to increased resolution due to disturbance by the radiation pressure, which should settle down eventually.

Pseudocode for AMR line transfer:

Start line transfer
if (processor at -x boundary)
  calculate radiative transfer (to child boundaries?)    ! only actually needs to be calculated on finest grid
  if (hasChildren)
    prolongate at boundaries and send to child           ! somehow need to be traveling along the ray while we do this, so that we encounter children in order (so that prolongation and restriction are correct)
  end if
  send ray to neighbor
else
  receive from neighbor
  calculate radiative transfer
  if (hasChildren)
    prolongate at boundaries
  end if
  if (not at +x boundary)
    send to neighbor
  end if
end if

Related: Have we ever thought about using the MPI distributed graph routines to keep track of neighbors?

HD209458b

Discovered that there's something incorrect about the stellar size output by the code (because it's still set by the stellar lambda, since it's not needed on the grid). Rerunning postprocessing on BlueHive, should have some new results tomorrow morning.

Radiation Pressure

Stampede is working now, even without increasing request descriptors, because it's running much more slowly, as Baowei was finding (~1 month to get to 300 frames). I think it's time to tackle AMR line transfer.

Also, we could theoretically run on BlueHive. Based on the new speed on Stampede, though, it would take months to complete a run.

HD209458b

Some very interesting ​results from the synthetic observations. Very little variation over time, as expected from looking at the simulation.

Here's the highlight:

We see significant symmetric absorption in the blue and red wings during transit (but not much before or after). Perhaps worth investigating a finer time series of these.

Radiation pressure

Communicating with TACC about missing tmi fabric. While loop appears to be working (should be blocking, but it's hard to tell when it immediately errors…), so should be working once the execution problem is straightened out.

BlueHive predicts 2.3 years, before being killed (due to long timestep), so it's not going to be effective to run even with a large reservation.

HD209458b short paper

Postprocessing of HD209458b run in progress. I'll finish some of the analysis in the next few days. Things to look at right now:

• Tau/xi terms (ratios of dimensionless Coriolis, tidal, wind accelerations)
• Regimes (should be energy-limited; probably type II from Matsakos?)
• Velocities and densities along sub-/anti-stellar ray, with markers for sonic and Hill radius (Coriolis?)
• Synthetic observations - high v?
• Compare all of these with most similar (high mass, low flux) planet from parameter space paper

Radiation Pressure

No real progress to report, other than that stampede is beginning to test my patience. But I believe I've solved all of the problems (except the MPI issue - more below), so we should at least start getting some frames.

MPI issues

Located all of our calls to IRECV and ISEND, and all matching WAITs are there. No MPI errors are going ignored. Is there anything else we can do to be convinced that the problem's not on our end? (Here's ​one example of a case where irecv is called by RECV, which could potentially cause the issues we're seeing but not be anything we can affect).

I've found the environment variable that controls the request descriptors (not sure why I couldn't find it before, but I distinctly remember not being able to set it), and upped it to 10¹⁴, so hopefully that will be sufficient to allow us to run for the full two days we're allowed.

HD209458b ionizing flux only

Got some analysis done here. Mass loss rates are about an order of magnitude lower, which is to be expected given the tighter binding of the atmosphere.

Mass loss rate: ~5.3x10⁹ g/s

It is (almost?) just Roche lobe overflow. Here's the Hill sphere (in red) overlaid on the wind, with contours of the dimensionless Coriolis length.

MHD planets

Here's the global steady state.

Global simulations are actually much worse than local. The run stalls out after 6 frames. Running for 2 more days, we get three more frames; here's the last one output:

Stellar wind implementation

I have a Matlab script that's solving for the (radial) velocity, pressure, and density of a pseudo-stellar wind with adiabatic gamma. I have noticed a couple of odd jumps that seem to be artifacts of the numerical solver, however. Also, I can't solve for any points inside the reference radius. This seems like a Matlab (or, specifically, vpasolve) limitation, but I haven't looked into other numerical methods in Matlab or otherwise. I've attached a sample solution file​.

Current list

I'll just keep this as a running tally, to keep us all organized a bit:

To Do

1. Continue radiation pressure runs
2. MHD planet runs (use the larger local runs as good-enough?)
3. Invert subcycling for line transfer
4. Add MHD to line transfer run
5. Test John’s stellar wind for charge exchange
6. Test charge exchange
7. Run charge exchange
8. Development for AMR line transfer/point source line transfer
9. Development for metastable He

Complete

Debug MPI issue on stampede (ish)
Set up global MHD run (w/ stellar wind)

Project/To-Do List

• Radiation pressure
  • Run high/med/low flux
  • Debug MPI issue (using too many descriptors)
• MHD planets
  • Fill parameter space
  • Add fields to parameter space run
• Charge exchange
  • Get initial conditions working (stellar wind without charge exchange populations)
  • Test charge exchange (have an attempt to compare to Christie et al [with fixed boundary code, though] - any other ways?)
  • Run (with suite of stellar wind strengths as in John's paper, perhaps? Should maybe look at what's predicted to affect the efficacy of charge exchange)

MHD planets

Refer to ​Jonathan's page for simulation parameters.

Using box of effective resolution 512³ for size (1 CU)³, with sphere of radius 0.25 CU forced to be resolved to level 1 (so it's equivalent to Jonathan's previous runs). Currently running corotating sims (do we want non-rotating as well?). Beta of 0.1 corresponds to surface magnetic field of 0.23 G.

Corotating

Have hydro steady state:

Also have a late state of the beta=0.1 run. May actually be a couple of frames too long, with some of the field lines being distorted down-orbit by stellar wind/planetary wind interactions? Running time is increasing, but only 4 days remaining when I restart (worth trying to get a few more frames?).

Beta = 10 run hit a wall at frame 224 - 1.9 months remaining.

Radiation Pressure

Run with no radiation pressure appears to have reached steady state. Actually not too different from frame 180 where radiation pressure is introduced in new runs.

Have a couple of frames from the high flux and medium flux cases. Just had stampede queue access restored, so will be restarting them.

High

Med

Debugging

Fixed makefile issue for debugging on stampede. Now that I have queue access again, should be able to use DDT to see if there's a leaky thread somewhere.

Current calculation

​Here I'm attempting to add the effects of ablation (more energy is absorbed by the layers being blown off than is required to unbind them), as well as the nonzero width of the bound torus (the gas at the inner edge of the torus will be more tightly bound than the gas at the orbital distance).

See ​these ​movies for some physical intuition. Note the steady state reached by the first - I've approximated this as two nested parabolas for the calculation.

Other approaches

Ruth's 2009 paper shows that radiation pressure alone would be insufficient to produce the mass loss rate determined by observations. Since this calculation is concerned primarily with ablation of the planet itself, and assumes a priori that radiation pressure can accelerate planetary hydrogen to the orbital velocity, I'm not sure it's relevant to our current situation.

Shaikhislamov et al. used the reconstructed Lyman-\alpha flux of ~10¹⁴ phot/cm²/s to calculate the pressure applied to the planetary wind: P = F h \nu / c \approx 0.5\times10^{-7} barye, which is about 2 orders of magnitude less than the pressure from their stellar wind (so they argue it's appropriate to neglect). Our torus is supported by a thermal pressure of ~8x10^-7 barye, and a ram pressure (in the x direction) of ~10^-11 barye (so insignificant compared to the thermal pressure).

Slowly-ramping (full flux around frame 80)

Quickly-ramping (full flux around frame 20)

As a reminder, here are the parameters for the planet:

Simulation domain: 20 R_p cube (0.0146 AU)

rho_p = 1.625d-15 g/cm³

T_p = 3000 K (not really important as long as wind gets hotter - T_wind should be >5300 K for this planet)

M_p = 0.73 MJ

R_p = 1.529 RJ

a = 65 R_p = 0.0475 AU

M_s = 1.23 MS (for tidal forces)

F_UV = 2d13 phot/cm²/s (Sun-like for orbital radius)

F_Ly-alpha = 4.1d14 phot/cm²/s, 4.1d15 phot/cm²/s, 4.1d13 phot/cm²/s (HD209458b flux, 10x, 1/10x - would it be better to do 5x and 1/5x?) (ramping up in approximately 1.5 CT, or 1/6 orbit)

• Need to move out of 476 for a week or two. What should we do with computers/where should we store miscellany?

• ​Parameter space paper essentially done. If we want input from Ruth, we should set up a time to talk with her. I'll give it a final once-over, then submit.

• First cut of intro and methods for ​radiation pressure paper. Introduction may need to go a little deeper into a discussion of the wind itself - could use some input here. Also, I've pulled most of the methods from the parameter space paper.

• HD209458b status: These two frames compare almost the same amount of physical time:

The wind is building much more slowly. Based on Owen & Alvarez 2015, we're well into the cooling-limited regime here. This seems relatively borne out by the heating and cooling happening:

Doing a calculation similar to that for the radiation pressure blowoff (but I believe more accurate - fewer unknown factors here), without cooling, we would expect the mass between the ambient and R_p to be completely unbound in \frac{E_G}{P_\gamma} = \tau_{launch} = \frac{4 \pi G \rho_p (0.331) M_p}{3 e_\gamma F} = [12310,1231,123] CT ([min flux, med flux, max flux]), compared to [263,26,2.6] CT for our heavier parameter space planet (which is on the high side, but not too far off - plus, we're not launching all of the mass considered here). This suggests it will take ~50x longer to launch the wind with HD209458 There's no explicit dependence on planet radius, but changing the radius also changes the density (which is why these differ, in addition to the additional planet mass). Louden et al.​ reconstructed a lower total ionizing flux (a factor of ½) than we're using, but we could increase the flux to get a stronger wind.

• I believe charge exchange should be working now, with the stellar wind also implemented as required. I have the stellar wind set to maintain in the -x, \pm y boundaries (y because there's the rotational component of the interaction). Waiting for the most recent test to start (requires more RAM than available on debug partition).

To Do

• Still contemplating the ​blowoff threshold. Eric has suggested cooling as a possible mechanism; I'd also like to include some sort of ablation, since the outside layers can absorb a lot more energy than required to unbind them before they're no longer shielding the inner gas (this is partially because we assume the gas absorbs at line center, no matter its radial velocity - this may be more important than we were initially thinking).

• Finalize charge exchange. Will test against analytic 1D solution; anything else we can test against?

• Continue HD209458b radiation pressure.

Poster​

Also, dress code for conferences?

Radiation Pressure

Still looking ok at 16 frames.

I've attached a Mathematica file​ with calculations for the balance point in terms of various quantities. The most questionable approximation, I believe, is that of constant density (although the radiation pressure isn't actually pushing on all of the material, I expect the momentum to transfer between the leading edge and the rest of the wind as they collide). With current parameters, the balance point is calculated to be about 10¹² phot/cm²/s.

Flares

Clearly one hour at 100x steady flux is too strong of a flare. Luca's suggested 40 minutes might be a more reasonable time, and I can turn it down, as well.

HD209458b Lyman-alpha flux

Got a response from Luca - the most accurate number for the Lyman-alpha flux at HD209458b's distance is 4.1x10¹⁴ photons/cm²/s (6.7x10³ ergs/cm²/s) (​Wood et al., table 3). So our flux here is about a factor of 2 lower than what is experienced by HD209458b.

Flares

Flare test running on stampede. Using Run 2 with flare parameters of:

Should be about 16 hours to run once I have it working.

Paper

Going through and editing, putting in references. Should be done by the end of the day - finished with introduction and methods.

Archiving old data

afrank_lab full- Baowei, any progress on moving the tarballs?

WASP-12b paper accepted - will be posting to ArXiv today.

Radiation pressure:

\alpha = F_{rad}/F_{grav}

\beta = F_{rad}/F_{ram}

\delta = F_{rad}/F_{press}

All of the following are in computational units:

F_{grav} = \frac{G \rho M_p}{r^2}

F_{ram} = \nabla (\rho |v|^2) \cdot \hat{x}

F_{press} = (\nabla P) \cdot \hat{x}

For \beta and \delta, white is F_{rad} \approx F_{other}.

\alpha

\beta

\delta

Cross-cut

Ly-\alpha flux for HD209458b

From ​here (Bourrier & Lecavalier des Etangs 2013), flux peaks at J = ~28x10^-14 erg cm^-2 s^-1 AU^-1.

Converting to photons cm^-2 s^-1,

F = \frac{J (1.5 \times 10^{13} \mbox{ cm})}{1.63\times10^{-11} \mbox{ erg phot}^{-1}} = 2.5\times10^11 \mbox{ phot cm}^{-2} \mbox{s}^{-1}

Which is 3 orders of magnitude less than the maximum flux in our simulation. If instead it should be scaled to the orbital radius (0.04747 AU, not 1 AU), that reduces it to 1.2x10¹⁰.

Timescale for blowoff vs. replenshiment

Assume torus of axis a = 0.04747 AU, radius of 12 R_p = 12(1.5\times10^{10}) cm, and uniform density \rho = 4.11\times10^{-19} g/cm³. Given the mass loss rate \dot{M} of the planet, the timescale for replenishment is:

\tau_{rep} = \frac{M}{\dot{M}} = \frac{\rho V}{\dot {M}} = \frac{\rho \pi^2 (12 R_p)^2 a}{\dot{M}} = 6.2\times10^6 s

The timescale for blowoff is dependent upon the gravitational binding energy and the energy deposited by photons on the edge of the torus per unit time:

E_G = 2 G \rho \pi^2 (12 R_p)^2 a M_p = 8.7\times10^{27} erg

P_{phot} = F A e_\gamma = 1.45\times10^{17} erg/s

\tau_{blowoff} = E_G/P_{phot} = 6\times10^{10} s

\tau_{blowoff} \gg \tau_{rep}

Radiation Pressure

Side

Top

Shadow recombination

\frac{\frac{dX}{dt}}{v_y} = \frac{dX}{dy} \propto \frac{(1-X)X(1+X)^{5.4}}{v_y}

Recombination rate \approx 15000;

\frac{dX}{dt} = \alpha (\frac{T}{10^4})^{-0.7} n_e n_{HII} = 7.56\times10^{-13} (\frac{2175}{10^4})^{-0.7} 84214000^2 = 1.56\times10^4

Recombination timescale = \frac1{n_e \alpha} = \frac1{84214000 7.56\times10^{-13}} = 15699 \approx 2\times \frac{R_p}{v_y} = 7360

Paper

• Run5 done:

• Non-rotating results and rotating run 1 results done. Need to discuss 2 further and write about 5 and 6.

• Want to get a good start on the discussion this week. Any sections I'm missing from the current outline?

Radiation Pressure

• 2d14 almost done. 2d15 crashed - I suspect everything's just been violently blown away, but VisIt keeps crashing every time I try to make visualizations of these. I'll try to figure that out this afternoon.

Submitted WASP-12b revision.

​Parameter space paper coming along. What still needs work:

• Intro: Have a sketch, but still needs citations and probably needs to be filled out more.
• Methods: Regime (Energy vs cooling limited)
• Discussion: Have an outline. Anything else that stands out from the figures?

Run6 for paper stalled, Run5 almost done.

Radiation Pressure

The higher-flux run looks good. Switched from plotting X to plotting 1-X, so log plots now work well (also changed this in the figures for the paper).

2d14, top

2d14, side

2d14, force comparisons

2d15, top

2d15, side

2d15, force comparisons

Comparisons of force from radiation pressure and magnitude of gradient of thermal pressure, and magnitude of gradient of ram pressure (not actually a force, doesn't seem useful for comparison):

Comparisons of ambient pressures in frame 150 and frame 151. Here is where the problem lies - ambient pressure is increasing due to OUTFLOW_WIND boundary conditions being set in a different way than before. This disrupts the steady state and masks the effect of radiation pressure on the stellar side of the wind.

Neutral fraction, temperature, density, and \frac{d\tau}{dx}, which is dependent on these. At edge of wind, between ~4 and 5, the temperature drops and the neutral fraction increases enough that most of the Lyman-alpha radiation is absorbed near the edge.

Exoplanets II abstract:
Using AstroBEAR, a 3D MHD fluid dynamics simulation code, we have developed a new platform on which to study the dynamics of planetary winds, including a variety of microphysics such as ionization heating, radiation pressure, and charge exchange. We will present the results of a study of the effects of planet mass and stellar flux on the structure of hot Jupiter winds, as well as changes to this structure that result from introducing radiation pressure.

​WASP-12b paper sent out to authors for review of referee response.

​Parameter space paper coming along. First cut of methods done, with a start on results. How much repetition is appropriate for figure captions/commentary? Intro is next on list of things to write.

Related, running out of space on BlueHive again. My personal quotas are full on BlueStreak (could request more - only have 200 GB here), BlueHive, and Stampede, and afrank_lab is full on BlueHive. We could back up and remove all but the final frames of the four low-flux cases (do we need the low-flux, high-density, low-pressure test for anything?).

On the simulations, the low-pressure ambient test for Run5 is proceeding nicely. The low-pressure ambient tests are running significantly faster than the higher-pressure ones - something to keep in mind when designing future planet simulations.

LyAlphaFlux = 0

This seems to indicate a problem. Running a couple of tests indicates that there's no radiation pressure in the first dozen or so steps, and the flux is zero throughout the grid, so it may just be that we haven't reached a steady state yet in Run2 (on further consideration, though, this doesn't seem likely to be the case - Run2 looks very stable between frames ~120 and 150, at least).

Two possible tests spring to mind - let test run go further with debug output in the line transfer routines, or run Run2 further with iLineTransfer=1, see if the same thing occurs. Some minor parts (that shouldn't matter on a restart) were changed between the original run and the radiation pressure run, primarily the number of zones between the outer boundary and the zero radius.

Radiation Pressure, Run2, top

Radiation Pressure, Run2, side

Planet density comparisons

​Go here for gif/stills; test2_side and test2_top are the final frames acquired (frame 97)(the GIFs don't make it there).

2c max ram pressure @ 1.5 Rp, frame 65 = 1.74x10^-4
2c max ram pressure @ 1.5 Rp, frame 20 = 1.04x10^-4
2c initial ambient pressure = 3.9x10^-8

2 max ram pressure @ 1.5 Rp, tfinal = 8.375x10^-4
2 max ram pressure @ 1.5 Rp, frame 98 = 7.77x10^-4
2 max ram pressure @ 1.5 Rp, frame 65 = 8.47x10^-5
2 max ram pressure @ 1.5 Rp, frame 20 = 2.7x10^-5
2 initial ambient pressure = 1.7x10^-8

Note that the original Run2 didn't have any buffer zones (pressure matching was essentially right at the outer boundary), which is why the ambient pressure is lower.

Sample Figures:

Tiled vs. Single Window

Axes, scales (units)?

Streaklines

Planet density comparisons

Run5 side

Run5 top

Run6 side

Run6 top

Updated ​WASP-12b paper to address referee concerns. Should new figure be included in the paper, or just the referee response?

Updated ​summary post with run statuses.

Started ​parameter space paper.

Radiation pressure:

​Blog post summarizing everything working so far for the parameter space study.

x10 flux

2x10¹⁷ flux (high)

Instead of using theoretical calculation, used theoretical relationship for 1D run - so upped planet density by only 2 orders of magnitude. Looks pretty good.

Radiation Pressure

This was before implementing ramping, without setting up the steady state wind first. Currently working on running from frame 150 of Run2 with ramping of the Lyman-alpha flux.

​Blog post summarizing everything working so far for the parameter space study.

Fixed bug in WASP-12 run - needs to restart from frame 10.

Highest-flux case working, but not launching wind. I expect the problem is resolving the ionization front:

Run1 flux, for comparison

Run3 flux

Run3 Lyman-alpha cooling

Run3 ionization heating

Run1 Lyman-alpha cooling

Run1 ionization heating

Finally, streak images are inverted:

Just invert Y-axis in matlab?

Simulations

​Blog post summarizing everything working so far for the parameter space study.

Started test of radiation pressure run - first few steps are struggling (pressure protections and high CFL). Once VisIt is done with mass-loss rates, will be able to see what's up.

WASP-12 run is at frame 116. Need to look at output - console only has level 0, but when I checked previously the planet was resolved at level 2. Also, execution time shot up at frame 116.

Last but certainly not least, adding a tracer to the ambient and stopping recombination for ambient material seems to be working fairly well at solving that problem (first frame in short time, no recombination in ambient yet), but the execution time still shot up - only got 3.5 frames in 5 days. Again, once VisIt is free I'll be able to take a look at the problem.

Physical Parameters

Simulation Parameters

lScale = R_P
rhoScale = rho_P

2b is higher-density planet surface with low flux.

Movies

noRot, Run 1 (side view)

noRot, Run 2 (side view)

Rot, Run 1 (top view)

Rot, Run 1 (side view)

Rot, Run2 (top view)

Rot, Run2 (side view)

Mass Loss Rates

noRot, Run 2 mass loss

noRot, Run 1 mass loss

Rot, Run 1 mass loss

Streak Streamline Images

Run1

Run2

Planet Closeups (top views)

Run1

Run2

Simulation Status

Run 3 testing

Looks like we need to take a different tack. The ambient is still ionizing too quickly and creating a shield in front of the planet. Lowering the flux (and therefore the density requirements) is one option - to get an idea of how much, look at the recombination timescale:

\tau_{rec} = \frac1{\frac{\rho_H}{m_H}\alpha_{rec}}

For \tau_{rec} \approx 30000 \mbox{s} (one unit of time), need ~5x10^-15 g/cm³ of (fully ionized) hydrogen. With a flux of 2x10¹⁷, the ambient medium is at about 10^-14 g/cm³. So reducing just one order of magnitude, to a flux of 2x10¹⁶, may be sufficient to solve the ambient recombination problem, since it would increase the recombination timescale to about 10x the simulation timescale.

Another option is to consider other adjustments to the ambient profile. If we were able to start the exponential profile farther in and still get a sufficiently low density at around R_p, that seems like a good option, but I can't seem to fiddle the profile into such a shape.

Run2 side view

Run2 top view

Other analysis products?

1. Movies (as above)
2. JC’s streak images that show flow pattern?
3. Total Mass loss rates
4. Back flow rates?
5. Observations of the kind in Carroll-Nellenback ea 2016?

Simulation Status

Attempting to run Run4 until transients die down. After 20 hours, timestep is still about 1x10^-7 - doesn't look like it's going to be successful. It's almost certainly related to the high ambient temperature of ~10¹³ K (not sure how much the transient flows are important in the final analysis - sound speed is very high).

Could try lowering density ratio based on recombination timescale (make sure it's large). Could also up resolution so planet gets closer to zero density, but this wouldn't be particularly effective in saving computational time (because of line transfer step - planet itself is only small fraction of grid, so high resolution hydro there is not a problem - which brings an idea - limit line transfer to something higher than max level?).

Movies in progress for corotating Run2.

Radiation pressure

Merged charge exchange and radiation pressure into local line_transfer branch. Need to test to make sure merge was successful, then can queue radiation pressure version of Run2.

WASP-12b Hill vs. bow shock radii

Using Matlab script for planetary parameters, we would need a stellar wind density (at base) of about the same as the planetary wind density:

\rho_\star = 4.95\times 10^{-11} \text{g/cm}^3

with

\rho_P = 5.881\times 10^{-11} \text{g/cm}^3

Flux to ionize Jupiter

An O-type star would need to provide about 1.2x10³⁸ ionizing photons per second to a planet at Jupiter's orbit for the equivalent flux. An O star actually provides about 10⁵⁰ phot/s, resulting in a flux of 1.3x10²¹ phot/s/cm² (or about 4 orders of magnitude more than our high-flux case above).

Referee report, ​WASP-12 paper

Moderate revision - only one real concern, that being the effect of stellar rotation on the stellar wind-planet wind interaction. Eric's comment that "If the result depended on the distribution of mass around the planet, then the point of the referee would lead to very important differences" is accurate, I think, and I do believe that's what we were attempting to get across in that section - so we could either revise that section for clarity, or we could run a new simulation (already set up) where the stellar rotation rate isn't locked to the orbital rotation rate.

The other concern was Fig 2, the 3D rendering. I feel it adds something, but I'm not sure I can make a strong argument for its continued inclusion.

Simulation Status

I plan to prioritize finishing rot_Run2 over the weekend (Jonathan - change frame numbers before restart?).

Now that I've looked at a few frames of Run3, the problem is apparent. The ambient medium is too dense, so no flux is making it to the planet. I believe there's a quick fix for this (extend the planet profile a bit farther, so that we get closer to zero density at the pressure-matching condition), although it may require upping the resolution at the outer boundary as well.

I do have movies for rot_run1 and noRot_run1/run2:

Rot, Run 1 (side view)

Rot, Run 1 (top view)

noRot, Run 1 (side view)

noRot, Run 2 (side view)

Rad transfer module test

Re-running the ionization-only, radiation pressure-only, and combined tests, they look good. The secondary effects in the ionization-only test are pretty strong still, but in general it appears to act like we expect. In the combined test, the effects of ionization dominate radiation pressure at the current fluxes (2x10¹³ ionizing flux, 5x10¹³ Lyman-alpha flux), but you can still see the effects of radiation pressure (particularly after the clump has been blown away).

Fellowships

NESSF - not sure about my chances, but worth applying? Mostly not sure how strong of a proposal I can write.

Simulation Status

Still need to look at 3&4 - I suspect the excessive cooling is related to the higher density in those runs, though.

Rad transfer module test

I've fixed the bug causing the 1D ionization test to fail:

However, the bug was such that it wouldn't have affected the clump test - and in fact the tests are identical pre- and post-fix:

• Ionization-only test: Matches exactly the results from our original ionization tests (way back when) up to the final frame at 4ish hours (~0.205 time units). Since this simulation is running in 3D, it doesn't go as far as the potentially troublesome simulation, but I think it's safe to go back to 2D.

The long term behavior still appears to not be what we expect — given the behavior of the 1D simulation (see below), it may be due to a higher flux than expected being used.

1D test: With the HotJupiter branch line transfer module, we get a stable ionization equilibrium (did rerun to make sure none of the parameters had been changed). With the radiation pressure implemented, the ionization front moves across the simulation domain to the right, indicating something incorrect happening. However, the momentum is a constant 0, so it doesn't appear to be radiation pressure-related.

• Investigated neutral tail coming from the first planet. Velocity vectors show it's originating almost entirely well outside of the planet boundary.

• Have simulations queued for parameter space. For production runs, using total of 150 frames for 15 computational times (5.25 days, or ~30 planetary crossing times).

Was missing the TABLES from the non-rotating version, not sure why I wasn't getting any error output.

• From discussion on Wednesday, see this page.

• For planet density in the high-flux case (assuming recombination-limited flows), see previous blog post.

Runs

Currently running small planet, low flux case (still recombination-limited, but looks like thanks to tidal forces the \tau=1 surface is staying well above the reset radius). Modified parameters (from first link above):

Also queued high flux, same size planet case. Python script​ was modified so that the planet profile is calculated assuming the totally-neutral \tau=1 surface is at the planet radius (though this is absolutely not the case for the new density). Modified parameters:

Code stuff

Got another soft crash on the small planet - need to investigate (race condition seems likely - Jonathan suggested layout_comms).

Running ionization-only test with radiation pressure implemented again to confirm and investigate results (second link above).

Further thoughts

Do we want to think about including a stellar wind at all? My thought: maybe for later simulations, but not for the parameter study, and probably not for the radiation pressure.

We had discussed possibly increasing the distance in y to capture more of the up- and downwind arms — may be nice for production runs at least (I wouldn't expect it to break anything, but could be worth testing first).

Constants:

R_p = 1.5\times10^{10} \text{ cm}

\gamma = \frac53

M_p = 0.263 M_J

T_p = 3000 \text{ K}

Assume every photon entering a tube of dimensions (R₀-R_P, dx, dx) is absorbed by the time it reaches R_P, and all the hydrogen in the tube starts ionized (X=1). Using an adiabatic atmosphere, the density profile is given by

n(r) = n_p(R_{crit}(\frac1r - \frac1{R_0}))^{\frac1{\gamma-1}},

with R_{crit} = \frac{(\gamma-1)G M_p}{\gamma c_p^2 R_p}, R_0 = \frac{R_{crit}}{R_{crit}-1}, and n_p the number density at the surface of the planet (to be solved for).

Volumetric recombination rate is

\frac{dn_0}{dt} (X,r) = (X n_H(r))^2 \alpha,

\alpha = 2.59\times10^{-13} (at the wind temperature of 10⁴ K).

Total recombinations (per unit time) are then given by

dx^2 \int_{R_0}^{R_p} \frac{dn_0}{dt}(X,r) dr

Total photons entering the tube per unit time are given by

J_{phot} dx^2,

where J_{phot} is the photon number flux.

Solving the equality​ of these two quantities for n_p, find

n_p = 1.171\times10^{-17}\sqrt{J_{phot}}.

For the currently-running case, our flux of 2x10¹³ gives a mass density of \rho_p = 1.2\times10^{-11}, much larger than our current value of 1.32x10^-16 (but we're not recombination-limited in the current case). For the high-flux case, J_phot = 2x10¹⁷,

\rho_p = 5.23\times10^{-9}.

• Summarized important elements of radiation pressure papers (from Luca).

• Now that compute nodes are back up, can restart simulation (equilibrium state to restart from) from where it was killed.

• Registration for next semester — Geometric Methods in Fluids?

Radiation Pressure

• Static radiation pressure test looks good. Setting up some dynamic tests. (Normalize min vel to sound speed?)

• Radiation/recombination-limited simulation (same as previous, but with 2x10¹⁷ max flux) needs to have ramping adjusted - larger flux just obliterates the atmosphere. I'm thinking start from 10^-4 or so lower, and adjust the ramping speed and half time accordingly.

• Need to work out photon-limited parameters — mass of planet as given is too low for the radius and temperature and results in an infinitely extended atmosphere.

First Set of Parameters

See here for an enumeration of the possible parameter space with orbital distance, flux, and planet mass variable. What follows is the top left of the first table.

See here for other simulation parameters.

Top view

Side view

Radiation Pressure

Testing with slab of density 1 (in CU), all neutral, embedded in ambient of density 0.01, all ionized. Setting J = \frac{n_{H} g}{\sigma_H p_\gamma}, the whole thing should be in a stationary state (the domain is thin enough that \tau over the grid is close to 0). Some problems with initial test conditions (I don't think I've sufficiently isolated the effects of radiation pressure).

I've restarted test with no hydrodynamics and constant temperature rather than pressure (\sigma_H depends on T), and the entire domain neutral. Seems like having hydrodynamics turned off removes gravity, though. From comparing diagnostic variables, seems like scaling is off on the radiation pressure term.

Scaling:

…and fixed.

Timing

See previous post for most simulation parameters.

Timing estimates that we discussed at Monday's meeting are inaccurate - I think there may have been a soft crash.

Wind

Line transfer evolving steadily. Still looks good at about two crossing times further than last post.

Top

Side

Fixed boundary temperature profile is backwards, but I don't think I'd expect that to affect the evolution other than reflecting it over the y=x line.

Also still working on box size/flux ramp test, but now I'm getting dt=0 on BlueStreak. Will be looking into this.

1D

After quick discussion with Jonathan, added spherical geometric terms (dilution and cooling) to code. Queued on BlueHive.

Line transfer

Stellar gravity and Coriolis test went well.

Parameters are the same for linked line transfer run, except:

Side view

​mesh

Top view

​mesh

Contour key:

Next

Need to test bigger box on BlueStreak (and perhaps get OpenMP working?). Since line transfer step happens at maximum resolution, the resolution of the box plays only a small role until a large portion becomes resolved. Pace is ~50 frames in 20 hours.

Fixed boundary

The equilibrium state is too close to the corrected Parker solution for any evolution to be seen. At 10^-11 CT, the wind is the same as at 1 CT. Nothing more to be done here.

1D comparison

Physical Parameters

Simulation Parameters

lScale = R_P

Got a run that seems to have reached steady state. Stays well below supersonic, though. Still need to get lineouts to compare with Ruth's '09 paper, but the problem seems to be that as the wind launches, the \tau=1 surface moves out so that the wind is launching from what seems to be a very different planet.

Next

Comparison of the ionization diagnostic variables should show why the wind isn't staying sufficiently ionized for an insignificant optical depth. Results should lead toward either futzing with the code or altering the simulation parameters.

Physical Parameters

Simulation Parameters

Fixed boundary (old code)

\lambda_P set such that R_P is equal for both runs.

TimeScale = P_orb
lScale = a

Contour key:

Not enough frames to capture evolution (where I suspect the tidal forces may be apparent). Requeued with 1/10th time. However, top view shows clear Coriolis effect. Looks perhaps a bit stronger but otherwise consistent with \tau = 0.08 and \Xi_P = 0.76 given ​Fig. 6.

Line Transfer

lScale = R_P

Contour key:

Unable to start from steady-state run, planet in incorrect position. At frame 760, there still appears to be some subsonic material at the edges, but this may just be because it hasn't yet blown out the entire ambient. The flow profile seems to be more strongly entrained up- and downwind than the above run.

Next?

May be best to run a (slightly larger?) box to steady state with the planet located correctly (using HD209458b orbital parameters - appropriate?), then add stellar gravity and rotation.

• Ionizing planet: see previous post
  • Two more things to (possibly) examine:
    • Increase box size by R_p in each direction, to match size of box in ATHENA
    • Whether ramping time or box size has bigger effect on \tau=1 surface
      • Should just need one more run (high resolution, small box, long ramping time)
  • Also, doesn't seem to be resolving anywhere but where it's forced

• Testing planet with Coriolis and stellar gravity (point particle off grid - ok?): BlueHive problems, only at frame 60. Wind isn't far enough from planet to be affected significantly by Coriolis, but tidal forces are definitely working.

Side View

Top View

• Still working on 1D version of planet

• Charge exchange: Have a modified fixed boundary condition module that should be testable, but working through a strange issue with the code.

• Conferences: I sent an email with a few options, reproduced below. Ruth and John have not made any plans for next year yet - they're going to meet to discuss and get back to us.

​https://exoplanets.phy.cam.ac.uk/Meetings/exoplanets2 - Early July.

​http://rencontresduvietnam.org/conferences/2018/exoplanetary_science/ - Last week of February/first week of March.

​https://www.cospar-assembly.org/show_infopage.php?info=52 - Mid July.

​http://eas.unige.ch/EWASS2018/program.jsp - Early April.

• Queued two more: -5 R_p to 5 R_p in all directions (same domain as ATHENA), at medium resolution (32 cells/R_p - BlueHive) and high resolution (64 cells/R_p) - BlueStreak.

• Also would like to compare a quick ramping time for one set of good runs, to determine if it's the ramping time or box size that has the most impact.

1 unit of computational time is 2 crossing times (for the planet - from 0 to Rp at the sound speed at the planet surface).

No movies for first two, due to storage difficulties.

64 zones/R_p, 2x10¹³ flux, quick ramping time (half flux by frame 30)

(Ignore the velocity data in this frame - I suspect it's corrupted).

64 zones/R_p, 2x10¹¹ flux, quick ramping time (half flux by frame 30)

64 zones/R_p, 2x10¹³ flux (total execution time 3.56 days)

32 zones/R_p, 2x10¹³ flux

32 zones/R_p, 2x10¹¹ flux

The low-flux case doesn't appear to have reached a steady state yet, possibly because the flow stays subsonic even with the larger box.

16 zones/R_p, 2x10¹³ flux

Lyman-\alpha tests

2x10¹³

2x10¹³, no Lyman-\alpha

2x10¹³, X=0

Conclusion: Lyman-\alpha cooling is related to the problem (perhaps just exacerbating), but not due to the interior of the planet collapsing any more (since there's no cooling in the case X=0 inside the planet).

Bigger box with corrected? BCs

Some …interesting… grid-aligned effects, plus the jump in temperature and pressure of the ambient (HII and ne increase, but that would decrease rather than increase the temperature).

Fixed the temps in at least one of these. ​Frames for all. VisIt didn't behave as I expected, so missing the first half of 2x10¹¹, 2x10¹².

2x10¹¹

2x10¹¹, no Lyman \alpha

2x10¹²

2x10¹², no Lyman \alpha

2x10¹³

2x10¹³, no Lyman \alpha

2x10¹³, X=0

Simulations to Run

• Photon fluxes of 2x10¹¹, 2x10¹², 2x10¹³ (maximum), with and without Lyman \alpha
• X = 0 inside planet, 2x10¹³ flux with Lyman \alpha

Original Results

Line transfer routines crashing was a problem with density and pressure protections not recalculating ionization state. Fixed (see ProtectionOptions).

Longer run doesn't look any better, though. Note that the temperature may not be exact (may be off by up to a factor of 2).

After 4 days

Initial Conditions

Initial conditions for the planet are finally correct. The ambient media still differ - the planet ends slightly earlier, and the profile is very steep near the edge.

Density

Temperature

Pressure

Ionization Fraction

Radiative Transfer

• The optical depths are close, but not quite identical. At ~300 s,

Working backwards, I find \sigma_H = ~6.28x10^-18 for AstroBEAR, and ~6.14x10^-14 for ATHENA - tested, and the difference isn't significant for the optical depth:

(ATHENA in blue, AstroBEAR in orange)

• Ionization/recombination rates and heating/cooling rates are qualitatively similar, but still differ by 1-2 orders of magnitude.

Ionization

(Without advective terms)

Theory Comparisons

It appears that the actual data from ATHENA doesn't match the theoretical calculation (using the data for n, X, and \tau in the equations on the LineTransfer page).

Heating

(Without advective terms)

• Electron densities are correct:

Dynamics

For the first time, the dynamics are beginning to look sort of similar. Radial velocities are too large outside of the planet for AstroBEAR:

And the evolution of the temperature on the day side is correct:

• ATHENA and AstroBEAR comparisons are working nicely. There are still a lot of differences between the simulations, however. Dashed lines are ATHENA, solid are AstroBEAR.

​Initial profile:

Initial conditions

Flattening radius - ATHENA flattens at r_ib rather than r_mask
Electron densities - ATHENA goes to near zero inside planet, AstroBEAR is just 10^-10 of the H density (this is probably not significant in the end)
r_ob - ATHENA profile extends to r_zero, but uses values at r_ob for density and pressure at edge
Pressure profiles differ

These last two seem likely to cause many of the differences seen once things actually start happening.

Frame 1

Looks like this ATHENA run was actually run with ramping - John seemed to believe it was not, so either the flag he's using isn't successful or he grabbed the wrong simulation accidentally
dE/dt for PdV, Lyman-alpha, and ionization heating all differ fairly significantly
AstroBEAR recombination is incredibly small - likely not attributable to differences in profiles

Other differences at this point are almost certainly attributable to the difference in flux or the difference in outer boundaries.

​Adjusted outward to r_zero:

With the profile extended to r_zero (and the resolution increased by one level), the location of the outer boundary becomes closer, but the values no longer match. So something is going wrong either in the calculation or import of the profile, still.

• Gave Diyora edge-centered VTK output from VisIt volume plot to start with. I almost have the layout working, but can't seem to get the refined data (only the coarsest is showing in the chombo). Is only level 0 data loaded BeforeGlobalStep?

• Did testing for OpenMP version of code. Problem appears sometime between 5 tasks/4 threads and 10 tasks/12 threads (haven't determined exact cutoff yet).

• Working on presentation for eclipse at Ogden Farmers' Library.

Both with recombination/cooling, 2e13 flux, no Lyman alpha cooling

32³ + 2 levels AMR

​frames

128³ fixed grid

​frames

Lineouts

Circles are 32³ + 2 levels AMR, triangles are 128³ fixed grid. The 128³ fixed grid is 200 frames shorter - the last frame is repeated. Day side looks nearly identical - night side differs, probably due to lower resolution.

Density

Ionization fraction

Optical depth

Temperature

Heating/cooling time scales

Heating/cooling rates

• Add diagnostics for cooling and recombination. ✓
• Subcycling time for planet. ✓
• Add flag for Lyman cooling. ✓
• Check cooling terms for correctness. ✓
• Ramp up flux over time. ✓
• Check recombination on night side. ✓
• With stored optical depth, calculate ionization and ionization heating rate (can't be done straightforwardly with a diagnostic variable). ✓

• Jonathan will be adding optical depth variable. ✓

Runs

• Lyman cooling, recombination/cooling, 2e13 flux ✓

​frames

• Recombination/cooling, 2e13 flux ✓

​frames

• Recombination/cooling, no flux ✓

​frames

• Lyman cooling, recombination/cooling, 2e10 flux ✓

​frames

• Calculate new profile for planet unaffected by Lyman cooling ✓

​frames

• Lyman cooling, recombination/cooling, ramping flux 1.2e13 → 2.02e13 ✓

​frames

• Recombination/cooling, ramping flux 3.6e12 → 2e13 - running

• Fixed minor typo in scaling of hydrogen cross section and added floor to line transfer dt. Pushed to line_transfer branch.

• Now appears to be working as expected (as far as absorbing the flux, anyway). The dynamics become more interesting around frame 50 (~4 hours).

3D run, 3 levels AMR, 64x64x64, initial ionization fraction 10^-5 (0 was crashing, but I believe that's been fixed), 1 computational time = 1 day.

Lineouts at 6.144 hours:

Movie goes to 13 hours.

• Why We Need an Entropy Fix​

• Implemented multispecies gas for Roe solver in AstroBEAR, but need a solution for density or sound speed for van der Waals gas to implement it. Relevant code snippet:

     ! Entropy fix      !!! Implement a different entropy fix to use van der Waals gas - density equations have complicated/no solution
     lambdaLL = WL(ivx) - PrimSoundSpeed(WL)
     lambdaRR = WR(ivx) + PrimSoundSpeed(WR)
     if (lambdaLL < 0 .or. lambdaRR > 0) then
      call StarPU(WL((/ irho, iE, ivx /)), WR((/ irho, iE, ivx /)), PrimSoundSpeed(WL), PrimSoundSpeed(WR), pstar, ustar)
      ! Passively advected (tracers; y,z momentum)
      ExactL = WL
      ExactR = WR
      ! Solution in interaction region
      ExactL(iE) = pstar
      ExactL(ivx) = ustar
      ExactR(iE) = pstar
      ExactR(ivx) = ustar
      ! Densities differ
      if (mpi_id == 0 .and. iEOS == EOS_vanDerWaals) write (*,*) 'This solution for density will not be accurate. The attraction parameter is not yet included in these calculations.'
      ExactL(irho) = PrimExactDensity(WL,pstar)
      ExactR(irho) = PrimExactDensity(WR,pstar)
      
      lambdaLR = ustar - PrimSoundSpeed(ExactL)
      lambdaRL = ustar + PrimSoundSpeed(ExactR)
     end if
     if (lambdaLL < 0 .and. lambdaLR > 0) then    ! Check for transonic rarefaction on left side
      print *, 'Using entropy fix.'
      lambda(1) = lambdaLL*((lambdaLR - lambda(1))/(lambdaLR - lambdaLL))
     end if
     if (lambdaRL < 0 .and. lambdaRR > 0) then    ! check for transonic rarefaction on right side
      print *, 'Using entropy fix.'
      lambda(nDim+2) = lambdaRR*((lambda(nDim+2) - lambdaRL)/(lambdaRR - lambdaRL))
     end if

• Wrote WAF scheme for own solver. Working, but results are not as good as I expect based on Toro's. Not yet sure where the problem lies.

• Related: Transmissive boundary conditions for second set of ghost cells don't quite make sense yet. Why move back inward to get wave to flow out?

• Finished edits of ​WASP-12 paper; reprocessed absorption.

• Still getting explosion (but better than implosion) with adiabatic atmosphere.

Should be the same as Jonathan's result - not sure why it appears to not be working.

I've added some data and a Matlab script for creating images such as those below. Density is on a linear plot- looks pretty good overall, although there's a strange wiggle near the inner boundary. Pressure is on a semilog plot, and is orders of magnitude off (I believe these should be in the same (computational) units).

EDIT: After running for longer time, not as stable as it first appeared (although I'm not seeing the same values in the lineouts as before, either). Perhaps the outer boundary is not as balanced as it should be?

Density plot (same scale as below):

Lineouts:

Pressure lineout scale changes dramatically after first frame, but can't get VisIt to save window correctly.

• Why We Need An Entropy Fix​

• Roe solver implemented in AstroBEAR:

Still need to implement van der Waals gas and figure out how to implement multispecies (and find tests).

• Testing charge exchange:

​frames

Seems to be working pretty much as expected, although the proportions of cold HII to hot H seem a little wonky - likely because of the other physics present, though.

• ​WASP-12 paper nearly there. Got some comments back from Luca, so I'll look at those in the next few days.

Link to some thoughts.

• UCSC - Reimburse which parts from what?

• WASP-12 simulations are looking good so far.

• Synthetic observations close to working (need loadbov).

• Roe solver without entropy fix written.

• Wrote multidimensional Euler solver, found here. Some highlights of the tests:

Comparing HLL to HLLC with stationary CD:

2D explosion test:

Strips in x and y only seem to appear after wave hits boundary - something to worry about?

3D explosion test:

Can't see quite as much - can't get VisIt to make slices of pseudocolor plots. But still looks pretty good.

The remainder of the tests can be found ​here.

• Anything else to do with code? Haven't implemented RCM, FVS, Roe, Osher, higher-order (WAF, finite volume, etc.) schemes, source terms.

• WASP12:

(~5 days left)

• Wrote Godunov solver for linear advection, inviscid Burgers equations, as well as Godunov solver for Euler equations (see here). Need to write up and comment Euler solver still, though not much is different. Also looking at a quick parallelization - need to check shared and private variables in loops. See ​here for test results.

• Implemented (but haven't tested) charge exchange - see changelog​.

• Working on higher mass loss rate planet. Think I was trying to change too much at once, so upped the stellar wind density a bit (bow shock radius is now at ~0.3). Currently having trouble reloading grids on restart - MPI errors.

• Wrote exact Riemann solver.

• Finished correct WASP run - time extended by ~100000 years, so not terribly significant. Updating ​paper figures (still need to find good output format).

Graph​ Volume rendering​

• Working on stellar fallback in ​WASP-12.

• Also have some analysis of WASP wind asymptoting. Took ram and thermal pressures of wind and set equal to disk pressure at required density at various radii - see Mathematica notebook​. For hydrogen extinction, need \kappa_R \lesssim 100 for HI at T = 10^4 K, \rho \approx 10^{-14} \frac{g}{cm^3}. This temperature is right on the border of H- and ff/bf opacity, so hard to say with certainty, but given the relatively low density it should be fine. (​link)

• ​Paper is in decent shape. Still needs some editing and a sentence here or there (figures and a few references, as well), but meat is done.

• Working on (done with?) application for Blue Waters fellowship

• WASP-12b:

Ran for 10 more orbits - planetary mass loss rate not asymptoting yet, but it does appear to be becoming more variable (just numerical?).

Qualitatively, disk is identical:

Also running with \alpha=0.1

• MHD: Running with betas for comparison with Owen and Adams - time growing exponentially again for stronger B, so need to look at that.

Need 2x10¹⁷ /cm² column density of magnesium, so ~2x10²² of hydrogen. Currently too small - by calculating average density in a box around the disk and multiplying by length of box, I find ~2x10¹⁸ /cm² of hydrogen. Would need to increase density by 4 orders of magnitude - plotting average density in this box over time, estimate 10^{5 orbits, or about 300 years, to build up sufficient density.}

They estimate that a mass loss rate of 3x10¹⁰ g/s for WASP-12b will give an upper limit of 4x10²¹ /cm² for Mg column density.

At 3.3 R_p, mass loss rate of planet is around 8x10⁹ g/s (appears to be decreasing over time, though), which is about 3.75 times lower than the mass loss rate used to estimate the column density - linear in mass loss rate, so this still gives an estimate of 10²¹ /cm². A rough estimate of the mass loss rate of the disk gives 5.1x10⁷ g/s, so two orders of magnitude smaller than the mass loss rate of the planet - which correlates well with an increasing overall disk mass (using slope of disk mass over time, get ~2.4x10⁹ g/s as a very rough estimate - within an order of magnitude).

​Images

MHD

Anisotropic MHD simulation at a surface beta of 50 seems to have gone pretty well:

​frames

​frames

And top view:

​frames

​frames

The magnetic field goes a bit crazy at the end, but the flow still looks about what we'd expect, with outflow from the night side highly suppressed and material there concentrated near the equator by the field. Best, we don't get the strange fourfold outflow and fallback. Seems narrower than the suppression seen by Owen & Adams, but that may be because of a difference in temperature profiles (they use radiative transfer).

I'm currently running a simulation with a beta of 10, which may lead to more coherent field lines - I can next run a simulation with similar field strengths to Owen & Adams, if the anisotropic simulations appear to be good.

WASP

​frames

​frames

The planet has only existed for half an orbit, so there's not much to see for troubleshooting purposes yet. But so far it looks as expected, with lobes spiralling up-orbit and down-orbit. Side view shows strange wave patterns starting in stellar wind around frame 7 - perhaps I need to increase the size in z a bit more?

• MHD begins to throw pressure protections between 1e-5 and 1e-2. They are thrown, as far as I can tell, immediately upon creation of the planet (went down to final time of 1e-100, with 100 frames). Up to this point, the fields and contours all look identical to the zero-field case, but want to check out early times for small B:

(Attempting to narrow down where protections start to happen I did a stupid - running three jobs again.)

10 frames of output posted on alfalfa for B=1e-2, to 0.01 orbits: ​MHD data

• Ran a larger WASP-12 simulation, but it also crashed - what exactly does changing the start time do computationally-speaking (for data output purposes)? It also runs until the planet is created - connection? But up to that point, stellar wind evolution is uninteresting/unchanged.

• Removed extra refinement from MHD simulations and extended boundaries by 0.025 in each direction. Still have subsonic material leaving grid in x-y plane, as well as both subsonic and sub-alfvenic eventually leaving in the z direction. Also, alfven speed is becoming very large in a concentrated area around frame 24 and increasing computation time significantly - not sure if there's anything to do about this.

Definitely going to increase grid size more. Anything else to do about computation time and irregularity of contours?

• WASP-12b simulations appear to have gone extremely well, on the other hand. Ran a test with new parameters for density, bow shock radius, etc.:

Mass loss rates may be a bit low - can increase both densities to get higher, but sims look good:

​Frames

Still a slight spiral feature in the initial stellar wind, before planet is created, but doesn't seem to persist (may just be Coriolis effects - doesn't appear to originate from stellar surface). Ran 9 full orbits - reaches quasi-steady state after about 3 days:

​Movie

• Finished GRFP application on Friday. Waiting on two letters of recommendation (as of noon today).

• WASP-12b: According to matlab script, no bow shock. Plot for ram pressures is

Initialized stellar wind out too far, so that planet was created inside:

​WASP-12b global simulation

• MHD simulations of HD209458b: attempting to control refinement. On most recent attempt, segfaulted:

​Error message

Relevant code snippet:

SUBROUTINE ProblemSetErrFlag(pos, Info)
    TYPE(InfoDef) :: Info
    REAL(KIND=qPREC), DIMENSION(:), INTENT(IN) :: pos
    REAL(KIND=qPREC) :: x(3), r
    x=pos-PlanetPos
    r=sqrt(sum(x**2))
!   Force lowest resolution outside of 5*Rp
    IF (r > 5*planetRadius) THEN
        Info%ErrFlag = 0
    END IF
  END SUBROUTINE ProblemSetErrFlag

Next thought is to write a separate subroutine that can be called in ProblemSetErrFlag, similar to the subroutines being called in the before step - but is this too simplistic?

• Also time to register for next semester - high-energy astrophysics or compressible fluids?

• Finally, 2D Owen and Adams MHD simulations​.

• Ran simulations with small and zero magnetic fields, as well as turning off the planet boundary (Parker wind out to bow shock) and lowering the CFL variables to 0.75, 0.1, 0.3. Runtimes below:

Seems the magnetic field and the incorrect bow shock radius are both playing a role, but bow shock radius is the bigger problem (since we need magnetic field). Some side-by-side comparisons:

​B=0, small

​Boundary off, B small

​Boundary off, CFL low

And alone, without all of the mess:

​B=0

​B=1e-10

​CFL lower

​No planet boundary

Some higher-resolution simulations:

​2 levels AMR

​2 levels AMR, with B field

​Full time

​Full time with field

• WASP-12 presentation​.

• Finished drafts of GRFP statements.

• Began working through Toro.

• Debugging MHD simulation runs. Starting with B_p = 0.54 G, B_s = 13.5 G, cutoff radii of 0.125 and 1*10¹⁰ respectively. First run, time to completion approached infinity - fixed this, but simulations are still extremely slow. Currently running short time, high framerate. Also calculated Alfven speed (\frac{B}{\sqrt{4\pi \rho}}) and plotted contours - high (3.5*10⁷ cm/s) near poles of planet. This seems to be higher than I would expect for the maximum Alfven speed, which should be near 5*10⁶ cm/s.

• Read three papers on WASP-12b observations (​Absorbing Gas, ​Metals in Exosphere, ​NUV absorption. Observations of WASP-12 have found lines missing from the magnesium spectrum, particularly, which are extremely unlikely to not be radiating from the star. In order to get absorption to the observed levels, it requires a density of a least an order of magnitude more than the interstellar medium. It is therefore hypothesized that WASP-12b is creating a cloud of gas that is capable of creating the observed magnesium absorption lines.

• Determining parameters for magnetic fields of planet and star for HD209458b simulations:

Matsakos et al. use the plasma \beta = \frac{8 \pi P}{B^2} = \frac{8 \pi \rho c_s^2}{B^2} at the base of the outflows (surfaces) to define the magnetic field strength. With \beta_s = 5, the magnetic field at the surface is 11.5 G, or ~10 B_s.

The magnetic field of HD209458b was estimated to be ~10% that of Jupiter (​Ref), or ~0.5 G. This gives \beta_p = 0.25, which is right in the middle of the range used in Matsakos (0.002 - 400).

Currently, calculations of the bow shock radius use \sigma. Will need to convert this to using \beta by relating magnetic field:

r_{bow}=\mbox{solve} \left [\rho_s(a-r_{bow}) v_s(a-r_{bow})^2 + P_s(a-r_{bow} ) +\frac{B_s^2}{(a-r_{bow})^6} = \rho_p(r_{bow}) v_p(r_{bow})^2 +P_p(r_{bow}) + \frac{B_p^2}{r_{bow}^6} \right ]

Or, in terms of knowns,

\rho_p c_p^2 \left\{ \phi(\frac{\xi_{bow}}{\xi_p},\lambda_p) [\psi(\frac{\xi_{bow}}{\xi_p}, \lambda_p) + 1] + \frac{8 \pi \xi_p^6}{R_p^6 \beta_p \xi_{bow}^6}\right \}=\rho_s c_s^2 \left\{\phi(\frac{1-\xi_{bow}}{\xi_s},\lambda_s)[\psi(\frac{1-\xi_{bow}}{\xi_s}, \lambda_s) + 1] + \frac{8 \pi \xi_s^6}{R_s^6 \beta_s (1-\xi_{bow})^6}\right\}

With the \beta values above, this gives a bow shock distance of 0.155 (in computational units).

(What does \beta from PlanetaryAtmospheres represent? \neq to plasma \beta)

The other parameter that needs to be set is the cutoff radius for the fields. Currently planning on large (essentially infinity) for the star and near the edge of the box for the local simulations (0.125 orbital radii).

Once these are set, I'll be ready to start running simulations. Plan to begin with isotropic temp profile and no rotation, then jumping to anisotropic with rotation (& the global model) if it works - are we interested in the intermediate models still?

• Also beginning to look at simulating WASP-12b, to attempt to reproduce the observed absorption lines of Mg II and Ca II. Physically, this requires that a fairly dense (10⁹/cm³ of H) and stable torus be created by the capture of the planetary wind. Will possibly require inclusion of radiative transfer, depending on level of ionization of torus.

• Passed prelim.

• Began application for NSF GRFP.

• Ran simulation with lambda=5 for longer time, to get to steady state:

• Numerical integration technique for isothermal Parker wind - starting from differential equation in terms of u and r, parameterize in terms of s, so that

We can calculate the Jacobian for this dynamical system, and evaluate it at the critical point:

Jacobian = \begin{bmatrix} 0 & 2c^3 \\ \frac{GM}{2c^3} & 0 \end{bmatrix}

This gives the linear approximation to the system near the critical point; calculating the eigenvalues and eigenvectors, travel a small distance along the eigenvector with positive eigenvalue (corresponding to tangent to the stable manifold of the original system) and use numerical integrator of choice outward from the critical point.

• Also started determining B fields for planet and star - looking for sigma of both to be 0.1 initially (based on Jonathan's example parameter space).

• Looked again at the Parker wind differential equation:

If we assume spherically symmetric outflow in an isothermal corona, then by conservation of momentum, continuity, and the ideal gas law, we have

Continuity shows that r^2 \rho u = C, where C is a constant. If we substitute for p and rho into the momentum equation, we get a differential equation in u and r:

Define r_c = \frac{G M m_p}{4 T} so that u_c^2 = \frac{2 T}{m_p} (i.e., sonic radius); substitute and rearrange to find

Integrate both sides to find the Parker wind equation,

• Modified the domain on a couple of simulations to better resolve the planet/capture sonic radius. Also plotted contour of sonic radius.

For lambda=5, may not reach quasi-steady state.

T_p = 5x10⁴ K:

M_p = 0.25 M_J:

Sonic radii:

• Continuing to work my way through Zel'dovich.

• Solving parker wind equation,

psi − ln(psi) = 4 ln(xi) + 4/xi − 3, where psi = (v/v_s)² and xi = r/r_s.

There are 4 types of solutions - combinations of sub/supersonic at the surface and v → 0/v_f as xi → infinity. The only physical solution is subsonic at the surface and has a nonzero velocity infinitely far away, and passes through xi = 1, psi = 1.

Attempts at solutions using approximations and a single solve:

Alternate implementation of approximation from Jonathan's code; works for r > r_s, but not r < r_s.

Pertinent line is: yy(i)=vpasolve(y - log(y) == -3 + 4*log(xx(i))+4/xx(i),y,(xx(i) > 1) * xx(i) + (xx(i) ⇐ 1)*exp(3-4*log(xx(i))-4/xx(i)))

Interpreted as:

if xx(i) ≥ 1

yy(i)=vpasolve(y - log(y) == -3 + 4*log(xx(i))+4/xx(i));

else

yy(i)=vpasolve(y - log(y) == exp(3-4*log(xx(i))-4/xx(i)));

end

Correct & incorrect plots:

• Calculated mass loss rate in VisIt. Used a scalar expression for the radial mass (momentum) flux and summed over isosurfaces (can also use spherical slices) to get mass loss rate. Plots were unfortunately uniform in display; however, in experimenting, discovered the flux operator makes plots that appear qualitatively correct for the flux:

Mass loss rates for the various simulations from last time:

Aniso appears correct, based on plot from planet paper.

Working on calculating the sound speed to get the sonic radius.

• Started working through occasional prelim problems, as well.

• Ran planetary simulations with no rotation and an anisotropic temperature profile with various parameters. See ​my page for density movies. Changing lambda performs as expected - with lambda = 5, there is a marked increase in the strength of the wind, and with lambda = 15 there is very little to no wind. Changing the mass of the planet didn't result in quite a clear changes, but overall it doesn't appear to have affected the wind very much. Hotter planet appears to have a stronger outflow, while the cooler planet doesn't differ qualitatively from the original run. The last simulations I've performed so far were decreasing the ambient temperature, and this doesn't appear to affect the wind very much in terms of density. Momentum plots, however, appear to show outflows that are stronger by a factor of two.

• Also researched justification for using a fluid approximation for charge exchange. In ​Christie et al., they calculate the mean free path (sound speed over rate of reaction) for charge exchange and note that it is less than the planetary radius (in regions dominated by the planetary wind) to justify using the fluid approximation, with mixing and exchange primarily at turbulent boundaries. In ​Murray-Clay et al., they justify the fluid approximation in general by comparing the scale height H to the mean free path of a particle, with the fluid approximation holding for H > lambda_mfp at the sonic point (in other words, the exobase [where H = lambda_mfp] is above the sonic point). Thus, H > R_p as well, and therefore greater than mfp of charge exchange near the planet.

• Finally, working my way through Zel'dovich and Raizer for hydrodynamics. Beginning viscosity and heat conduction at the moment.

• Read ​Tremblin and Chiang for computational charge exchange. The paper is a followup to earlier studies of charge exchange between planetary and stellar winds, which used Monte Carlo simulations of particles. Here they use the hydrodynamic equations (no magnetism, Coriolis and centrifugal forces, or tidal gravity); the isothermal planetary wind was initialized as 80% ionized, following Murray-Clay et al. The planetary wind incorporated photoionization/recombination and advection. To incorporate charge exchange, the hydrodynamic code was augmented with chemical reaction solvers, where beta is the reaction rate.

These equations take reverse exchange into account, so as not to overestimate neutral hydrogen (still slightly overestimated). They are very similar, but not identical to, the equations used by ​Christie et al.

The simulations appear to reproduce the observed absorption curves well, with asymmetry between the two sides of the Doppler shift.

• Used Jonathan's Matlab code to examine change of bow shock radius with magnetic fields. If sigma_* and sigma_p are equal, the bow shock radius is unchanged with or without magnetic field - ratio of radius to orbital separation, chi_bow = 0.240468. With sigma_* = 1, sigma_p = 0.1, chi_bow = 0.148204; sigma_* = 0.5, sigma_p = 0.1, chi_bow = 0.187300; sigma_* = 0.1, sigma_p = 0.5, chi_bow = 0.302483; sigma_* = 0.1, sigma_p = 1, chi_bow = 0.363674; and with sigma_* = 0.5 and sigma_p = 1, chi_bow = 0.297793 ≈ chi_Coriolis.

• Also attempted to recreate isotropic planetary wind with no rotation (Run5 from planet directory). Copied .data files, but clearly didn't turn out correctly. Need to figure out why.