Using AstroBEAR for Lab Experiments
AstroBEAR needs more physics in order to make it more applicable to lab experiments. I've been toying with this idea, and I don't think I'm going to get very far because it is quite complicated. However, I have been somewhat successful at hard-coding some 0th order approximations into a separate problem module (not usable by the rest of the code). Even the act of attempting to implement these things has forced me to learn about new physics, yay!
Laser Deposition
A simple dump of some fraction of laser energy onto a target surface (resonance absorption). Target heats up and the thermal pressure gradient launches a plume of plasma.
One of the tricky things here was the directional nature of the laser heating. Once the laser hits some critical density, it is only absorbed to a certain depth. Thus far, my code only works with single proc and a fixed grid.
A more clever implementation could make use of multiple processors and AMR. And to do laser deposition more accurately, you'd need ray tracing.
Now, the fun MHD stuff…
Biermann Battery
In ideal MHD, we have the induction equation:
.
In resistive MHD, the induction equation becomes:
This is in some version of the code thanks to Shule. There are other MHD terms that can be important depending on the context of the problem. Here is the induction equation with the Biermann term:
The Biermann term makes self-generated B-fields possible. In other words, you don't need an initial B to generate more B. This is also in my new problem module, and could someday be implemented into AstroBEAR for general use. This formulation assumes the plasma is fully ionized (i.e.
).
Nernst Effect
Lastly, we have the Nernst term. Here's how the induction equation becomes even more complicated:
This
is the Nernst velocity, and it can be written as:
where
is the heat flux and is the thermal conductivity. is technically a tensor since it is different in different directions (anisotropic). We can simplify the tensor and write as follows:
where
is a unit vector in the direction of . The first two terms are implemented somewhere in the code for conduction, but the third term is missing. Furthermore, I believe the conduction implementation in the code is strictly for transporting heat; in other words, there is no magnetic field generation.This is something that I haven't implemented yet. There should be a way to add the Nernst velocity to the fluxes within the code, but I think I will first try to separate
in the induction equation and treat it as a source term.Mach Stems: Cooling Strength ==> Effective Gamma
Below is a new set of 2-D runs designed to explore how the strength of radiative cooling affects Mach stem formation and size. An approximation of the critical angle for Mach stem formation depends only on the adiabatic index gamma. We hypothesize that the cooling strength implies an "effective" gamma such that very weak cooling implies a gamma of 5/3 and very strong cooling implies a gamma of 1.
Furthermore, if a Mach stem forms, its size is limited to the smaller of the cooling length and the clump diameter. Thus, as cooling strength increases (decreasing cooling length), the size of the Mach stem should decrease.
Here are the important parameters:
vs = 50 km/s M = 5.2 T = 8322.56 K nclump/namb = 5000 tfinal = 100 yrs Rclump = 10 AU
dcool / Rclump | Ambient Density [cm-3] | d / Rclump | Cooling Run | Effective Gamma Run |
---|---|---|---|---|
30 | 3.18 | 5 | ![]() | 1.46 ![]() |
10 | 9.77 | 4 | ![]() | 1.25 ![]() |
3 | 33.95 | 3 | ![]() | 1.19 ![]() |
1 | 108.21 | 2.4 | ![]() | 1.14 ![]() |
0.3 | 388.55 | unstable | ![]() | ⇐ 1.10 ![]() |
0.1 | 1252.4 | unstable | ⇐ 1.10 |
Then, there are 5 more runs with M = 30 and T = 250.047 K.
dcool / Rclump | Ambient Density [cm-3] | d / Rclump | Cooling Run | Effective Gamma Run |
---|---|---|---|---|
30 | 3.634 | 5 | ![]() | 1.50 ![]() |
10 | 11.16 | 4 | ![]() | 1.30 ![]() |
3 | 38.76 | 3 | ![]() | 1.18 ![]() |
1 | 123.30 | unstable | ![]() | ⇐ 1.10 ![]() |
The effect of cooling is stronger than I had anticipated, and I'm not sure if we can use Pat's figure from his recent Mach stem paper to make quantitative comparisons. However, qualitatively, we are already proving our point, and more simulations with lowered gamma instead of cooling can lead to values for effective gamma.
Meeting Update 05/11/2016 - Eddie
- Finished thesis last week! woohoo!
- Clumpy paper and response to referee is written. Ready to resubmit?
- Cooling paper is written. Need more details on other codes. Pat will add some writing in 1-D radiative shocks section?
- Resumed 3-D pulsed jet sims. Hydro and beta = 1 runs are done. beta = 5 is 95% done, beta = 0.4 is 66% done and may take a few more weeks to finish on bluehive. Maybe only 2 weeks if I can get a reservation for more nodes?
- Finished poster for next week's HEDLA conference:
More 2-D Mach stem runs w/ cooling
Previous blog posts: ehansen04122016, ehansen04192016
Did a few runs with d/Rclump = 6 where the previous runs were closer together at d/Rclump = 4.
Rclump/Lcool | At 50 years | Mach stem size |
---|---|---|
2 | ![]() | regular reflection |
1 | ![]() | 0.3 Lcool |
0.5 | ![]() | 0.4 Lcool |
Also did a run where we are no longer in the frame of reference of the clump and shock; both ambient and clump are moving now, but relative velocity is the same as before. With d/Rclump = 4 and Lcool = 0.05 Rclump. It looked the same, i.e. no instability.
Meeting Update 04/19/2016 - Eddie
- Updated table on Mach stems and Lcool: ehansen04122016
- Checked to see if low resolution is cause of bow shock instability: doesn't appear to be, image below
- Checked if 3-D effects cause bow shock instability: plausible, image below
- Checked if even stronger cooling (Lcool = Rclump / 100) is the cause: also plausible, movie below.
- My thoughts: any small perturbation can trigger the instability whether it be numerical or physical. The stronger the cooling, the easier it is to perturb. I believe my simulations in the clumpy paper are more susceptible simply because they are 3-D; extra degree of freedom for motions and/or noise. Also, in the paper, we are not in the frame of reference of the fast clump.
Lcool = 0.1 Rclump, 2-D, low res
2-D, nice Mach stem, Lcool = 0.5 Rclump movie
2-D, unstable?, Lcool = 0.01 Rclump movie
- I should be able to finish up referee report and resubmit clumpy paper by end of week.
- Pat is working on cooling paper, and I will add some more details this week.
- Thesis is coming along nicely. Most of the writing is finished except for the conclusions section.
Meeting Update 04/12/2016 - Eddie
- HEDLA XI accepted my abstract, gave me a poster presentation
- While writing up the introduction for my thesis, I decided it would be nice to have a schematic of a disk-jet system so I drew one.
- In response to referee report on the 3-D clumpy paper, I am running some 2-D sims to look at Mach stem size when cooling is present. For each simulation, I change Rclump to see if Mach stem size is limited by Lcool !and/or Rclump.
All runs have M = 7 and Lcool = 12.7 AU.
1 computational unit = 1 Lcool.
separation distance d = 4 Rclump
Rclump/Lcool | At 50 years | Mach stem size |
---|---|---|
100 | ![]() | regular reflection, unstable? |
20 | ![]() | regular reflection |
10 | ![]() | regular reflection |
2 | ![]() | approximately 0.5 Lcool |
1 | ![]() | approximately 1 Lcool |
0.5 | ![]() | forming single bow |
So it does appear that when cooling is strong (Rclump/Lcool ≥ 1), the Mach stem size is limited by Lcool or it's even smaller. However, I'm not convinced that when cooling is weak, the Mach stem size will be limited by Rclump.
Update 03/28/16 - Eddie
- 3-D pulsed jets are still running, only hydro run is finished. Beta = 1 run will be done in about a week on bluehive, beta = 5 is 75% done and beta = 0.4 is only 50% done on davinci. The queues on davinci take way too long, so I'm going to move the beta = 0.4 run to bluehive as well. It'll take a while since I can only get 120 cores, but at least it'll be making progress.
- Ready to finish response to referee report on clumpy paper. Need to have meeting with Pat.
- Finished first draft of cooling paper, ready to revise and add figures. Working on some cooling curve plots now.
Update 03/11/2016 - Eddie
- Worked on cooling paper: abstract, intro, methods sections finished. Might've finished the whole thing had I not focused so much on recent referee report. I can release these sections for edits, and continue working on the rest.
- Worked on referee report for clumpy paper. I think most of it is resolved, and I don't think we have to change the paper too much. However, we still need to have a co-author discussion about gammas, cooling, and Mach stems. Below is a test simulation I ran to show that simulating a more "realistic" clump does not change our conclusions.
- 3-D jet simulations are still going. Here's a status update:
Beta | Frame (out of 120) | Est. Wall Time Remaining |
---|---|---|
Inf | 103 | 10.7 hours |
5 | 76 | 2.1 days |
1 | 88 | 17.4 days |
0.4 | 55 | 4.2 days |
The beta=1 run will take significantly longer because it's on bluehive, but the other runs wait in the queue for so long in between restarts on davinci that all of these runs may finish around the same time in 3 weeks or so. Honestly don't see this project finishing in time to be part of my thesis, especially since we wanted to do additional runs with random pulses, less frequent pulses, no pulses, etc.
Update 02/18/16 - Eddie
- Monday, I was busy gathering/generating some old data for Pat for the Mach stem paper. Revisions according to referee report should have been submitted yesterday.
- Been putting my thesis together and writing the cooling paper which will also be a chapter in my thesis.
- Running the 3-D pulsed jets on davinci. Might move to bluehive since the queue on davinci is taking forever. Below is the beta = 1 run. I changed the scaling on the emission to make it brighter. It doesn't look like there is enough resolution to separate the H-alpha and [S II] lines within a clump, but we do see H-alpha (green) leading the clumps (yellow).
Update 02/08/16 - Eddie
- Clumpy paper is done. Time to submit?
- Writing cooling paper, and thesis.
- Running 3-D jet sims. The only thing I don't like so far is that the magnetized runs seem to have pressure (probably magnetic pressure) pushing radially outwards at the injection region. See images and movies below…
Update 01/25/16 - Eddie
Nothing new to post in terms of images or movies, so here's an update on what I'm currently working on:
- Clumpy paper: Received another round of comments from Pat. Time to discuss these and figure out how best to edit the paper.
- Cooling paper: Writing up the Methods section which includes a bunch of new equations that I haven't included in any of my other papers. I'm giving myself an arbitrary deadline to have a completed first draft by the end of February.
- 3-D jet simulations: These are running on Rice's Davinci. Looks like it takes about 1 week per run (not including wait times).
- Thesis: Started putting this together. Not much writing done yet, but I have a formatted template now.
- Job search: I'm going to give a talk at the LLE on Wednesday Feb. 9th. Still need to decide which of my projects I want to present. Also, need to schedule an interview with Paul Drake's group (U. Mich.)
I think that's everything. I'm pretty busy these days…
Update 01/11/16 - Eddie
- Finished another round of edits for 3-D clumpy paper
- Got the 3-D pulsed jets running faster now.
- On 576 cores, I got the MHD run w/ beta = 1 to run 94 out of 120 frames in just 8 hours.
- This was at a lower resolution of 32 cells per jet radius, and the jet propagated to a distance of about 23 jet radii.
- I suspect this run would take about 6 more hours to get the full 120 frames. Meaning that a run at the desired 64 cells per jet radius should take a little over one week to finish. This is doable.
Still having issues with Visit on Grass. I hope Rich can fix the issue soon, so that I can make movies more easily.
Update 11/30/2015 - Eddie
- I've applied to several jobs now, but need to do more.
- Working on rewrites for paper. Plan on finishing it by the end of this week.
- MHD version of 3D pulsed jets was having some restart issues, but I fixed it. Below is an image from a run with an initially purely toroidal field with beta = 1, with 64 cells per jet radius. This run (6 out of 120 frames) took 3 hours on 576 cores. The current standard output estimates that the run will take about 3 weeks to finish. Perhaps it's a better idea to do a longer run at a lower resolution (32 cells per jet radius) before continuing with this production-level run. I just don't want to waste time and resources if there turns out to be some difficult-to-explain-weirdness in the jet at later times.
Update 11/16/2015 - Eddie
- working on edits for 3D bow shock interactions paper
- applying to jobs
- still trying to figure out waviness in 3D precessing jets
- Things I have tried so far:
- primitive limiters instead of characteristic
- lower CFL
- outflow object inside lower z boundary
- velocity fade in radial direction
- no pulses
- H viscosity
- apply diffusion
- smaller precession angle
- Noticed some potentially odd things:
- jet boundary looks tilted in first frame, likely due to precession in magnitude of the velocity in z-directio
- the magnitudes of the velocity in every direction x, y, z seem to be precessing. I don't think this is what we want. I think vz should be constant in time (except for the pulsations of course), and vx and vy should change but in a way such that at any given time they are constant throughout the outflow object. In the current implementation, there appears to be a phi dependence on magnitude.
- Things I have tried so far:
Update 11/09/2015 - Eddie
- Made image from HH 2 for my 3D clumpy paper. This is to compare with the multi-clump simulation since both show a shock "sheet".
- We wanted to do another comparison figure with HH 34 to show lateral motion, but this is looking a bit trickier than I thought it would be because the .fits images are not aligned with each other. Waiting to hear back from Pat to see if there is an easy way to deal with this.
- Having issues with Grass and with my account on the Rice machine, so I don't have all of the data from some of my new 3D pulsed jet runs. Below are images from a 0 deg precession and a 1 deg precession run. I have simulation data for a run with increased refinement around the outflow object, and a run where the outflow object is inside the bottom domain boundary; I just can't get to the data right now. Below is an image comparing runs with 0, 1, and 3 deg precession from left to right respectively.
- Lastly, I worked on my CV and I've been busy for the past few days applying for some national lab jobs.
Update 11/02/2015 - Eddie
Using ds9 and gimp to generate a couple of figures for my 3D clumpy paper that will show an HST image and a simulation image side-by-side. Below is an example:
3D Jets
I've tried several things, and I still see the waves in all cases:
- For some reason, the exact Riemann solver does not work well. Code stops early due to nan in flux.
- Added velocity fade from 0 to Rjet. This is not a precession fade, just a velocity magnitude fade.
- Turned off characteristic limiters, use primitive limiters instead.
- Turned on ApplyDiffusion.
- Lowered target CFL from 0.3 to 0.1.
- Turned off pulsations.
Below is an image from the run with no pulses:
Update 10/19/2015 - Eddie
Below is image/movie showing the hydro 3-D pulsed jet: slice of x-z plane, slice of y-z plane, and column density.
It might be helpful to review how the outflow object is set up: OutflowObjects
My jet uses a radius = Rjet/2, thickness = Rjet/4, base = Rjet, and open_angle = 0. The outflow object steps on some cells, but the rest of the lower boundary is extrapolated; should probably change this to reflecting.
Things To Try
- use different Riemann solver
- adjust CFL
- adjust boundary conditions of outflow object (radius, thickness, fade, etc.)
- turn off pulsations
- have precession fade from 0 to Rjet
UPDATE
Things i have tried thus far:
- For some reason, the exact Riemann solver does not work well. Code stops early due to nan in flux.
- Added velocity fade from 0 to Rjet. This is not a precession fade, just a velocity magnitude fade. Still see waves.
- Turned off characteristic limiters, use primitive limiters instead. Still see waves.
- Turned on ApplyDiffusion.
- Lowered target CFL from 0.3 to 0.1.
- Turned off pulsations.
Update 10/13/2015 - Eddie
2D Mach stem results with low M
low M = 1.55
TMS, SB = Transient Mach Stem, develops into a Single Bow shock
SMS = Stable Mach Stem
RR = Regular Reflection
gamma | d | Result |
---|---|---|
5/3 | 24 | RR |
5/3 | 23 | SMS |
5/3 | 22 | SMS |
5/3 | 21 | SMS |
5/3 | 18 | SMS |
5/3 | 12 | SMS |
5/3 | 11 | SMS |
5/3 | 10 | SMS |
5/3 | 9 | TMS, SB |
1.4 | 21 | RR |
1.4 | 20 | SMS |
1.4 | 19 | SMS |
1.4 | 18 | SMS |
1.4 | 17 | SMS |
1.4 | 15 | SMS |
1.4 | 12 | SMS |
1.4 | 11 | SMS |
1.4 | 10 | SMS |
1.4 | 9 | SMS |
1.4 | 8 | TMS, SB |
1.4 | 6.5 | TMS, SB |
1.4 | 6 | TMS, SB |
1.4 | 5 | TMS, SB |
1.4 | 4.5 | TMS, SB |
1,2 | 18 | RR |
1,2 | 17 | SMS |
1,2 | 16 | SMS |
1,2 | 15 | SMS |
1,2 | 13 | SMS |
1.2 | 12 | SMS |
1.2 | 9 | SMS |
1.2 | 8 | SMS |
1.2 | 7 | TMS, SB |
1.2 | 6 | TMS, SB |
1.2 | 5 | TMS, SB |
1.2 | 4 | TMS, SB |
3D Pulsed Jets
Ran the MHD, beta = 1 model farther at 32 cells/Rjet.
Started a production-quality hydro run at 64 cells/Rjet.
3D Clump/Bow Interaction Paper
- Having some trouble with ds9 and the .fits images. I just want images that look like the ones from Pat's 2011 paper, but I can't seem to get the scaling quite right.
3D Pulsed Jet on DAVinCI
I tested the new pulsed jet module on Rice's machine DAVinCI, and it seems to work well. This pulsed jet is 3D, MHD, with non-equilibrium cooling.
The code ran pretty well, producing chombos at a steady pace with no restarts. Here are some numbers to show how well this run did on DAVinCI:
Resolution | 32 cells/Rjet |
---|---|
Cores | 576 |
Run Time | 2.5 hrs |
Cost | $0.49 |
% Run Completed | 30% |
The run does slow down as it moves along, so it's hard to say how long the full run would take at this resolution. The code estimates 11.2 more hrs after frame 30, but the actual time is probably closer to 24 hrs. If I assume that the time between chombos grows linearly, then I can estimate the total runtime to be about 24 hours.
Thus, the total run at 32 cells/Rjet would take 24 hours on 576 cores, costing 4.70. This is the maximum number of cores that I can use on DAVinCI, and I can only submit for a maximum of 8 hours at a time. So even at this low resolution a run is going to take more than a day if you include wait time and restarting. If we double the resolution to 64 cells/Rjet, this run could take 2 weeks or more instead of 1 day, and the cost would be more than 50.
UPDATE
Two new models: hydro (no MHD), and an MHD with beta = 1 (same as above) but with precession. Precession amplitude is 3 degrees and the precession frequency is 2.5 times the pulsation frequency.
Update 09/21/2015 - Eddie
New 2-D Mach stem runs with lower M
I tested gamma = 1.4 with 4 different separations: 4.5, 5, 6, 6.5 rclump. I chose these separations because they are right on the edges of the different regimes (see figure below).
Velocity was lowered from 50 km/s to 15 km/s resulting in Mach number being decreased from 5.18 to 1.55. I had to lower the maximum of the density scale from 35 to 16 since these are weaker shocks (number density in computational units of 10,000 cm-3).
The table below summarizes the Mach stem behavior for the relevant models after 75 years. To characterize this time in terms of clump-crossing times, it would be about 40 for the original M = 5.18 runs and about 12 for the new M = 1.55 runs. (TMS = Transient Mach Stem, SB = Single Bow, SMS = Stable Mach Stem, RR = Regular Reflection)
d (rclump) | M = 5.18 result | M = 1.55 result |
---|---|---|
4.5 | TMS, SB | TMS, SB |
5 | SMS | TMS, SB |
6 | SMS | TMS, SB |
6.5 | RR | TMS, SB |
As you can see, with a lower M, we now get a single bow in every model. I think this is because there is less pre-shock ram pressure to keep the enhanced post-shock thermal pressure at bay. In regards to Mach stems, I see two possibilities: 1) Achieving a SMS is impossible, and the behavior will transition to RR at larger separations, or 2) the SMS region has just shifted to larger separations. Either way, I want to do some more runs to see if I can find the SMS regime for this lower M.
Other Stuff
- Debugging the error in the projections. The error I'm getting now is not where I expected it to be which is confusing. It is in data_declarations.f90 where the Info object is initialized. Maybe I need to run with more cores to get by this part of the code, and trigger the real error?
- Still haven't received any confirmation about my guest account on the Rice machines.
- I'm in the process of moving some data off of BlueHive to make room for some 3D pulsed jet runs.
- PAPER EDITS!!! There is a lot to be done, so I'm mainly busy with this now.
Update 09/14/2015 - Eddie
- re-applied for a guest account on Rice machines, waiting to hear back from Pat
- I'm working on some low M, 2-D Mach stem simulations for Pat. This is for the referee's report; to see how low M affects Mach stem formation in simulations.
- Working on re-generating some of the paper figures for my 3-D clump/bow shock paper. Mainly want to increase resolution of images, but I'm having some issues. See images below…
3D Clumpy Paper Figures
I had some issues regenerating the projections at higher resolutions, so I just used the old projections. I think this is fine since you can't see the high resolution in a paper anyways. I did get it to work for the one zoom-in image which is arguably the only one that this is important for. Below are the 11 figures going in the results section; the other 4 figures which I did not post here are in the theory and methods sections.
fixed resolution issue for figures
I fixed the resolution issue. Turns out that when you use a camera object, the code automatically sets the resolution to be equal to the maximum number of cells in the x direction. Thus, if your y or z directions have a larger extent, your data will be coarsened.
Some of my 3D clump simulations have a y domain that is almost 3 times the size of the x domain. This is why I said that the resolution should be 80 cells per clump radius, but the images show approximately 30 cells per clump radius.
We should probably change this in the code so that the camera resolution matches the actual maximum resolution of your data. The old image is on the left, new one on the right.
Update 08/18/2015 - Eddie
- applied for guest account on Rice machines
- we will be using DAVinCI: 2304 cores on 192 nodes (12 cores per node)
- scratch directory for running jobs, work directory quota is 2TB per group, can use Globus for file transfer
- can use 48 nodes (576 cores) at a time, but only for 8 hours at a time
- there might also be some fees, but I'm not sure yet. From what I found online, it might be a flat annual fee of 100.00 + 340.00 per million SUs
- made some figures for the 3-D clump paper
- Adam, I'm taking all of your suggestions from your latest email, and adding one more figure. In addition to "orientation comparison", I want a "viewing angle" comparison figure.
- of the figures below, the enlarged one requires the most discussion…
i = reddened, post-shock [S II], cooling regions
ii = bow shock intersection point, possibly Mach stem
iii = post-intersection point "cone"
iv = unstable, fragmenting bow shock
v = the "froth"
Update 08/13/2015 - Eddie
- waiting to hear back from Andy about getting on Rice machines
- need to start generating figures for the results section of the clump/bow shock paper
- need to start some higher resolution 3D pulsed jet runs
- worked on some 1-D radiative shock models for Pat since he supplied us with a new cooling table
- The new z cooling is very very close to the old z cooling in terms of total energy losses. It appears to be slightly stronger or slightly weaker depending on the exact parameters.
- DM cooling is clearly much weaker in the range of the z cooling table, but stronger at high temperatures.
- I did a convergence study to see how resolution affects the initial post-shock temperature of one of my models. The model appears to begin converging for resolutions > 30 cells/AU.
Update 07/28/2015 - Eddie
- new movies: MachStems
- best ones to look at are 2-clump comparisons
- the column density movies are useful for watching the transmitted shocks move through the clumps, and the movies also support the idea that the secondary "messy" emission is due to material shedding off of the clumps
- working on the paper; mainly the theory section to include discussion of various instabilities in these sims.
- Adam, I'll send you an up-to-date version of everything I've written so far on Friday?
- I could continue making emission map and column density movies, but I don't know how much more information we'll get. I have generated 55 movies already, out of a potential total of 92. I'll do some more, but I might not do all 92. I think my focus should shift to reading and writing now. I could also start some new 3-D pulsed jet runs.
- I'm also starting to pull out some specific details from Pat's 2011 third epoch paper. One thing that stuck out was this image:
In the paper, Pat writes "One can also imagine the structure in the froth as arising from the irregularities in the shape of the wing of the bow shock as it sweeps up material". The bow shocks in our simulations show irregularities due to the strong cooling, and I believe this is the Vishniac instability (a thin layer bounded by shock and contact).
Update 07/20/2015 - Eddie
Emission Maps
- Run C is the only one left that needs emission maps.
- More movies posted: MachStems
Clump Instabilities
- 4 processes: transmitted shock compresses clump, RT, KH, radiative cooling
- cloud-crushing time: where is the clump/ambient density ratio, is the clump radius, and is the incident shock velocity
- RT time: , and the dominant mode is of order the clump radius, so the RT instability will grow on time scale equal to
- KH time: , and again the dominant mode is of order the clump radius so the KH instability will grow on same time scale
- Radiative cooling is not occurring within the clump because the clump is already very cold (2.5 K) to satisfy pressure equilibrium with the ambient. Even with the transmitted shock, the clump will not reach a temperature large enough to trigger radiative cooling (T > 1000 K).
M | (yrs) |
---|---|
7 | 40.66 |
10 | 28.46 |
15 | 18.98 |
The simulation time is 75 years, and references state that a clump is typically compressed in 1-2
Bow Shock Instabilities
- Vishniac instability: typically cited for blast waves. Occurs when a slab bounded by a shock and a contact discontinuity is thin enough.
- NTSI: usually referred to in colliding flows. A cold slab bounded by two shocks will undergo NTSI if the slab is thin enough.
Update 07/13/2015 - Eddie
Not much to add since late last week, just a few things:
- progress on generating the rest of the emission maps is steady and on-going (~1 week until they are all finished): MachStems
- updated the page with a couple of helpful diagrams
- I'll put up some more movies tomorrow.
- will also start generating some column density movies this week.
- re-reading Pat's 2011 3rd epoch paper in closer detail
Update 06/29/2015 - Eddie
- Astronum and my talk went very well. I got a lot of good compliments, and talked with a lot of people.
- Amira vs. Disperse paper with Federrath?
- Finished more 2-D Mach stem runs for Pat and the paper: ehansen06222015
The complete set of models result in the following figure:
- Work continues on generating all the emission maps for my 3-D clump runs. Here is the updated table: MachStems
- I have a reservation on BlueStreak right now to finish up run E. It should be done by next Monday. Then, I just need to finish up the emission maps for runs D and E.
- For next week, I will start generating movies of all the emission maps that I do have.
- My goal for July is to finish the 3-D clump paper and start some new 3-D pulsed jet simulations.
new 2D Mach stem runs
With gamma = 1.4, clump separation = 5.5 rclump, higher clump/inflow density contrast of 2e5. 200 frames for a total simulation time of 150 years.
With cooling, clump separation = 3 rclump:
New gamma = 1.2 runs, with d = 3, 4, 5:
New Rotated Emission Maps
I implemented rotation for my emission maps, so now they can be generated at any angle rotated about the y-axis (phi-direction).
A couple of the emission maps were initially created without using a camera object such as the 0 deg inclination and 0 deg rotation map. However, I found that it was difficult to plot this next to the maps that were created with a camera object. The camera objects change the scaling in a non-trivial way and it's not easy to plot them next to non-camera-object-generated maps in visit. So I'm regenerating a couple of emission maps so that they all use a camera object.
Below are the images/movies for one of the 3-clump models (model K).
Meeting Update 05/18/2015 - Eddie
- Update on status of 3-D clump runs: MachStems
- Need to add some lines of code to rotate emission maps about y-axis (phi direction). Then, I will generate some new images/movies for one of the models (3-clump model K), and we'll see what we get before doing this for all models.
- inclination angles = 0, 30, 60 deg
- rotation angles = 0, 45, 90 deg
- Emission movie from multiclump run:
- Running a 2-D Mach stem simulation (the one with separation distance d = 5.5 rclump, gamma = 1.4) with higher clump/ambient density contrast (106 instead of 104) and for a longer time (150 years as compared to previous 75) to see if Mach stem grows and then smooths into one continuous bow shock.
- AstroNUM abstract?
Meeting Update 05/11/2015 - Eddie
Mach Stems
- Fixed an error in my Mach stem problem module: ehansen05072015.
- Now, I'm rerunning all the simulations, and so far so good. Here is an update on the status of these runs:
- I'm hesitant to keep using Stampede since we're down to 300,000 SUs and Marissa needs to use it. So my plan is to finish some of the Stampede runs on Comet, but Comet is down at the moment.
- My reservations on BlueStreak and BlueHive just started, so now I can have more things running at the same time.
Other Stuff
- Rich and I are in the process of moving some of Shule's old data off of our machines to make more space.
- AstroNUM is just around the corner. We need to finalize this abstract, so I can start putting a presentation together.
- All indications show that the new pulsed jet module works, so as soon as the Mach stem runs are done I can start some higher resolution pulsed jet runs.
new emission map after fixing error
There was an error in the code I was previously using such that the ambient wasn't being initialized with the correct ionization fraction. I have since fixed the error and am in the process of rerunning all the sims. I made a movie of comparing the old emission with the new emission.
This is run J from the table of runs: ehansen05042015. Baowei is running A, B, and C on Gordon. Runs D, E, and F are approximately 50% done on Stampede. Run K just started, runs L and M will start once my BlueStreak reservation begins. Runs G, H, and I are sitting in the queue on BlueHive. If they don't start soon, then they will start once my reservation starts on that machine.
Update 05/04/2015 - Eddie
Mach Stem Runs
Below is an update on the status of all the runs. I added an "Inclination Angles" column to denote if the emission data at different inclination angles (0, 30, 60, and 90 deg) have been produced for that run.
Run | # Clumps | Defining parameters | Machine | Frames (out of 50) | Inclination Angles |
---|---|---|---|---|---|
A | 2 | d = 2.5, M = 7, 10 | Stampede | 50 | |
B | 2 | d = 4.5, M = 7, 10 | Gordon | 49 | |
C | 2 | d = 6.5, M = 7, 10 | Gordon | 50 | |
D | 2 | d = 2.5, M = 7, 15 | Gordon | 23 | |
E | 2 | d = 4.5, M = 7, 15 | Gordon | 26 | |
F | 2 | d = 6.5, M = 7, 15 | Gordon | 24 | |
G | 3 | d = 2.5, a = 30 | BlueHive | 44 | |
H | 3 | d = 2.5, a = 60 | BlueHive | 50 | |
I | 3 | d = 2.5, a = 90 | BlueHive | 30 | |
J | 3 | d = 4.5, a = 30 | BlueStreak | 50 | Yes |
K | 3 | d = 4.5, a = 60 | BlueStreak | 50 | Yes |
L | 3 | d = 4.5, a = 90 | BlueStreak | 50 | Yes |
M | 10 | BlueStreak | 45 |
I expect the multiclump run (run M) to finish tonight. Once the remaining 3 clump runs finish on BlueHive (runs G and I), that will free up the nodes to do the different inclination angles. The Gordon data has been transferred to BlueHive where it will have to stay because we don't have any more space on our local machines.
The remaining Gordon runs are taking too long, and the machine appears to be having file system issues. I'm going to move these to Stampede and finish them there. However, we may not need many more frames for the 2-clump, M=15 runs (runs D - F). Below is an image of run D at frame 22/50.
The bow shock has already moved past the slow clump, so perhaps this is enough. Runs E and F are a little different. It takes longer for the fast bow to move past because the separation is greater. Run E might be fine with 30 frames, and run F may need the full 50. Doing "partial" runs like this will speed things up.
Below is the emission image/movie from run A.
Jet Module
I fixed the EMF profile in the jet module, so the magnetic field looks correct now. Below is 60% of a run that I ran on BlueHive on the afrank nodes. It won't run any further for some reason, but I think it's an issue with the afrank nodes, and not the code.
The run below has the outflow object resolved to 32 cells per jet radius, and the rest of the grid is resolved up to 16 cells per jet radius.
Updated Outflow Module
The outflow module has been updated with a few new features, so that I can use it for my 3D, pulsed, magnetized, cooling jets. I basically added 3 new profile objects to the outflow object:
- SpeedProfile: Before, you could only give the outflow object one outflow speed, and you could make it oscillate sinusoidally if needed. Now, the outflow object will first check to see if there is a SpeedProfile. This is a profile of the outflow speed as a function of time. I put this in here so that I can eventually set up a profile of random velocities to use.
- TempProfile: The original outflow object just takes one temperature for the entire outflow. Now, with the TempProfile, you can give it a profile of temperature with respect to radius. This was necessary to use the hydromagnetic pressure equilibrium profile that I have used in my jets.
- MagProfile: The original outflow object took one value for the peak toroidal magnetic field. Now, you can pass it a MagProfile which a profile of toroidal magnetic field strengths as a function of radius. The functionality for other field geometries is not there yet, but I, and many other projects, would only need toroidal fields anyways. We can implement helical fields later.
I built a new PulsedJet problem module which uses this new outflow object set up. In my module, I define an appropriate temperature and magnetic field profile. I have not done the random velocity profile yet, but that should be pretty straight forward. I plan on putting all of these profiles into profile_functions.f90 so that anyone can easily use them.
Below is a very low resolution example of a 3D, magnetized, pulsed, cooling jet. It is a 2D slice of density.
The bow shock seems pretty wide to me, and I have not yet figured out why. I have checked my temperature profile several times, and it is correct. Perhaps there is some confusion between thermal pressure and internal energy within the outflow object. I'm going to investigate this further just to make sure everything is working accurately. The magnetic field part looks like it's working perfectly, and I haven't seen any restarts yet.
This simulation is pretty slow since I'm just using 16 procs on alfalfa. I'll let it run further in order to see if there are any more bugs that need to be worked out.
UPDATE
I fixed a bunch of stuff related to the new outflow module:
- Got the velocity pulses working.
- Edited the new pulsed jet problem module to handle hydro (non-MHD).
- Edited profiles.f90 to make it more user-friendly.
- Fixed some issues with the magnetic field part of the outflow object. This involved a time dependent issue, and also replacing use of MagProfile with !EMFProfile.
Number 4 above has been the trickiest to solve, and I'm not convinced that we have it right yet. One issue is that you cannot step on magnetic fields because that could create divergence. So what the code does is add a certain amount of field to the outflow object. However, you can't just keep adding your initial magnetic field profile because then the field will build up and you'll have too much magnetic energy/pressure.
The amount of field you need to add is equal to the amount of field that has advected out of the outflow object. Everything in the outflow object is done in terms of the EMF, so the EMF source that needs to be in the outflow object is as follows:
This should, in theory, take care of time-dependent issues. One other issue that remains unresolved is that the initial Bphi profile in the chombo does not match the Bphi profile that I supply.
Many of these points look ok, but I'm not sure why the profile has a discrepancy around r = 0.5. This might be solved by increasing resolution or subcycling.
Below are low resolution simulations with the current implementation. Much better than before, and maybe good enough, but I want to figure out the Bphi issue first.
Update 04/20/2015 - Eddie
The 3D runs for the Mach stem paper are going on several machines. Here is the status of all the runs (d = separation distance in rclump, a = angle in degrees):
Run | # Clumps | Defining parameters | Machine | Frames (out of 50) |
---|---|---|---|---|
A | 2 | d = 2.5, M = 7, 10 | Stampede | 43 |
B | 2 | d = 4.5, M = 7, 10 | Gordon | 16 |
C | 2 | d = 6.5, M = 7, 10 | Gordon | 24 |
D | 2 | d = 2.5, M = 7, 15 | Gordon | 15 |
E | 2 | d = 4.5, M = 7, 15 | Gordon | 11 |
F | 2 | d = 6.5, M = 7, 15 | Gordon | 11 |
G | 3 | d = 2.5, a = 30 | BlueHive | 0 |
H | 3 | d = 2.5, a = 60 | BlueHive | 32 |
I | 3 | d = 2.5, a = 90 | BlueHive | 0 |
J | 3 | d = 4.5, a = 30 | BlueStreak | 36 |
K | 3 | d = 4.5, a = 60 | BlueStreak | 34 |
L | 3 | d = 4.5, a = 90 | BlueStreak | 15 |
M | 10 | BlueStreak | 0 |
Notes:
- The Mach numbers are now all integer numbers: 7, 10, and 15.
- The "slow" clump remains at M = 7 for all runs. This means that the runs that contain an M = 15 clump now have larger domains to accommodate the faster clump with the same slow reference frame.
- The 10 clump run is coded up now. I will test it and run this one after all the others are finished.
- Data transfer is holding me up. I need to get the chombos/bovs from Stampede to Clover so I can visualize but this is taking a very long time.
- I have reservations on BlueHive and BlueStreak starting on Wednesday at 12:01 AM. I could only get 10 nodes on BlueHive and 512 nodes on BlueStreak for one week. This is not much compared to what I have access to already, but I won't have to wait in the queue or do so many restarts which is nice.
Below are column density maps and emission maps from run A.
Other Stuff
- I will do some more paper edits today and tomorrow and post them.
- I submitted an abstract for the CIRC poster session, and I have my poster so that's all set.
- I typed up a quick abstract for AstroNUM and attached it to this blog post. Haven't submitted it yet, please let me know what you think.
Meeting Update 04/13/2015 - Eddie
Not much to report today, because I'm still waiting on a bunch of simulations.
- Had issues with the Stampede file system over the weekend, so I had to restart the 2-clump, Mfast = 10 runs. They are still waiting in the queue, but they are almost 50% done.
- On bluehive, I'm running some low resolution tests before moving to production runs:
- 2-clump, Mfast = 15, Mslow ~ 11 which will replace the Mfast = 5, Mslow = 3.66 runs
- 3-clump, Mfast = 15, Mmed = 10, Mslow = 7
- Will continue writing the sections of the paper that I can without results
Mach stem papers
There are only 2 relevant papers on ADS with "Mach stem" in the title:
- My HEDLA conference paper: http://adsabs.harvard.edu/abs/2014arXiv1412.6495H
- Kris' hysteresis paper: http://adsabs.harvard.edu/abs/2013HEDP....9..251Y
More results are found by searching for "Mach reflection":
- Kemm 2014, interesting paper on double Mach stem problem: http://adsabs.harvard.edu/abs/2014arXiv1404.6510K
- Ivanov 2002, nice reference on regular vs. Mach reflection, experimental and numerical: http://adsabs.harvard.edu/abs/2002ESASP.487..341I
Other references refer to optics, water waves, detonation waves, etc. I don't think these are very applicable to our scenario. Papers 1-4 might have enough references within for me to use. Paper 4 has led me to some purely fluid mechanics papers:
- Chpoun 1995, wind tunnel experiments regular vs. Mach reflection, http://adsabs.harvard.edu/abs/1995JFM...301...19C
- Hornung 1982, more on transition from regular to Mach reflection: http://adsabs.harvard.edu/abs/1982JFM...123..155H
- Ivanov 1995, more on hysteresis effects: http://adsabs.harvard.edu/abs/1995PhFl....7..685I
There are a few books that give good descriptions of shock reflection phenomena which I have used as references in my previous paper as well:
- Ben-Dor: Shock Wave Reflection Phenomena
- Courant and Friedrichs: Supersonic Flows and Shock Waves in Applied Mathematical Sciences
- Landau and Lifshitz: Fluid Mechanics
There is also a derivation of the critical angle in this paper:
- de Rosa 1992: http://adsabs.harvard.edu/abs/1992PhRvA..45.6130D
Ben-Dor is definitely seems to be a key player in shock-wave reflection. He has a lengthy review paper just on shock wave reflection in different regimes:
- Ben-Dor 1988: http://adsabs.harvard.edu/abs/1988PrAeS..25..329B
I've skimmed through these papers, but I will be reading them more closely over the next few days as I start to write my paper. Some key things that I have learned:
Historically, there have been two minimum critical angles for Mach stem formation based on two different criterion. The one that I have been using is known as the "detachment" criterion (alpha_d), and it is commonly cited in text books. The other criterion is the "von Neumann" or "mechanical equilibrium" criterion (alpha_n), and it is a pressure argument which leads to an angle less than the detachment criterion. Many papers talk about a dual-solution region where the angle is alpha_n < alpha < alpha_d. In this region, both regular reflection and Mach reflection are theoretically possible. The figure below illustrates this point:
Ben-Dor's review paper suggests that for steady flows, the mechanical equilibrium criterion is the correct one to use, and that regular reflection only occurs above alpha_n in what he calls psuedo-steady flows. One complication that I see is that alpha_n is much more dependent on M than alpha_d. Thus, strong shock approximations might not be appropriate, and I would have to use a critical angle formula that is also dependent on M. I have yet to find a good formula for that like I did for alpha_d.
Ben-Dor also has an region showing no reflection in his figures, but there doesn't seem to be a good explanation of the transition criterion or a useful formula anywhere.
Experimental results vary. Some agree with alpha_n, and some agree with alpha_d, while others show hysteresis effects and support the dual-solution region argument.
In conclusion, I might be quoting a minimum critical angle of ~ 40 deg, when in actuality it could be more like 30 deg. The issue could be avoided for Mach numbers ~ 2.2, since alpha_n = alpha_d at this M. However, then I would have to use the exact formula for alpha_d that contains M since a strong shock approximation would no longer be valid for M=2.2.
Sensitivity of H-alpha Routine
The H-alpha routine in the code determines the amount of H-alpha emission (neutral hydrogen de-excitation line from level 3 to level 2). It is dependent on electron number density, temperature, ionization fraction, and hydrogen density. There are 2 main components to this emission: a recombination term and a collisional excitation term.
Collisional excitation is what most people usually think about when talking about H-alpha emission: electrons collide with neutral hydrogen and the excited hydrogen atom goes through the H-alpha transition. The recombination term comes into play when you have ionized hydrogen that wants to recombine. When it recombines, there is a nonzero probability that it goes through the H-alpha transition. So, if you have a lot of ionized hydrogen at a low temperature, there should be a lot of recombination and the recombination term can dominate the H-alpha emission.
Below are two plots showing the H-alpha emission as a function of temperature for two different ionization fractions. The hydrogen number density is the same in both plots, and for simplicity the gas was assumed to be all hydrogen. nH = 5560 was chosen because it is the average of the post-shock densities of two simulations which I will talk about.
Notice how in the x=0.99 plot, the recombination term is much stronger and dominates the H-alpha at low temperatures longer. This is because 99% ionization is unrealistically high for these temperatures, so the gas should rapidly recombine and thus go through many more H-alpha transitions.
3D Mach stem simulations
For reference, see ehansen02242015.
For the 3D Mach stem simulations, we noticed a significant decrease in H-alpha emission from the M = 10 runs to the M = 5 runs. The best way to explain the difference is by comparing the temperatures in a couple of these runs. Left is M = 5, right is M = 10, both at separation distances of 3 rclump.
I measured the temperatures in the bow shock intersection regions and got about 6900 K for M = 5, and 19,700 K for M = 10. This is approximately a factor of 3 which seems reasonable given that postshock temperatures are approximately proportional to vshock2
These simulations have relatively low ionization fractions throughout, and at the bow shock intersection points, the ionization remains near 0.01 for both runs. So we should look at the x = 0.01 H-alpha plot.
It would appear that these simulations have hit the sweet spot. There is a sharp increase in the H-alpha going from 6900 K to 19,700 K. This is where the collisional excitation term takes over and dominates. I estimate a jump of about 2 or 3 orders of magnitude which is consistent with the emission maps.
Also, look at the maximum temperatures in the simulations, and then look at the x = 0.01 H-alpha plot. This explains why the M = 5 model does not show a wide range of H-alpha intensities like the M = 10 model does.
In conclusion, the H-alpha routine is working just fine. The parameters I chose for the simulations just happen to illustrate the temperature sensitivity in our code.
Proposed Runs for 3D Mach Stems
2-clump runs
We vary the separation and Mach numbers.
Run | Separation (rclump) | Mfast, Mslow |
---|---|---|
A | 2.5 | 10, 7.32 |
B | 4.5 | 10, 7.32 |
C | 6.5 | 10, 7.32 |
D | 2.5 | 5, 3.66 |
E | 4.5 | 5, 3.66 |
F | 6.5 | 5, 3.66 |
On bluehive, I used about 3800 SUs per run for 40 cells/rclump resolution. Adding a level of refinement to get 80 cells/rclump will increase SUs by approximately a factor of 16 (a little less due to filling fractions). Accounting for the fact that Stampede cpus are slightly faster than bluehive (2.7 GHz vs. 2.4 GHz) leads to approximately 52,000 SUs per run. For 6 runs, this is a total of 312,000 SUs.
3-clump runs
Clumps stay at M = 2.5, 5, 10 for all runs. The M = 2.5 and M = 5 clumps remain in same plane at a fixed separation distance, and the M = 10 clump is positioned at some angle and distance from the M = 2.5 clump.
Looking down the barrel, the 3-clump setup is like this:
Run | Separation (rclump) | 3rd Clump Position (deg) |
---|---|---|
G | 2.5 | 30 |
H | 2.5 | 60 |
I | 2.5 | 90 |
J | 5 | 30 |
K | 5 | 60 |
L | 5 | 90 |
We would get longer run times because filling fractions are higher with 3 clumps. My best guess is that these runs could take about 77,000 SUs each. A total of 462,000 SUs for 6 runs.
10-clump run
Clumps at random positions with varying separations and velocities, but all moving in same general direction.
With many clumps on a large grid, this could be a very expensive run. I'm guessing it would be at least 100,000 SUs if we still want 80 cells/rclump.
In total, my crude estimates give on the order of 900,000 SUs for the entire study. Again, we're playing with a factor of ~16 here if we want the increased resolution. So if we keep the resolution at 40 cells/rclump, then we can do all these runs with less than 60,000 SUs. Maybe we only do the 2-clump runs at the increased resolution and the rest at the original, lower resolution.
Meeting Update 03/31/2015 - Eddie
2D Mach Stems
Finished measuring actual angles of intersection: https://astrobear.pas.rochester.edu/trac/blog/ehansen03202015. Notice how these angles do not follow a trend.
Below is a side-by-side image of the gamma = 1.4, d = 6.5 and d = 3 models. Just by eyeballing it, it appears that the angles are not significantly different. The predicted angles should be approximately 20 deg different, whereas the actual angles only differ by about 8 deg. I conclude that these actual initial angles do not determine the result. Perhaps an angle at a later time would be a better predictor, but I don't see how to explain that theoretically.
Development Meeting 03/24/2015
Development for AstroBEAR 3.0
Below is a table of all the development tasks that we want to complete for AstroBEAR 3.0. An "X" denotes that the task is completed.
Ticket | Summary | Development | Testing | Documentation | Tested by test suite | Owner |
---|---|---|---|---|---|---|
#150 | 2.5D and 1D Self Gravity | Erica | ||||
#154 | 2D and 2.5D Sink Particles | X | Jonathan | |||
#243 | Implicit Heat Conduction | X* | Baowei | |||
#253 | Script for Automatic Testing on Local Machines | Jonathan | ||||
#264 | Parallel IO | X | Jonathan | |||
#305 | Krumholz Merging Algorithm | Shule | ||||
#321 | Simple Line Transfer | X | X | Jonathan | ||
#356 | Flux Limited Diffusion | X | X | Jonathan | ||
#359 | ZCooling | X | X | Eddie | ||
#360 | EOS Cooling Routines | Jonathan | ||||
#362 | Refinement Module | X | X | Eddie | ||
#363 | Refinement Buffering | X | X | Jonathan | ||
#366 | Disk Module | Jonathan | ||||
#367 | Emission Routines and Cylindrical Projections | X | X | Eddie | ||
#370 | Outflow Module | Eddie | ||||
#375 | MUSCL in multi-dimensions w/ AMR | X | X | Jonathan | ||
#380 | Update Problem Modules for New Test Script | Eddie | ||||
#416 | Restarts, regridding, static levels, maxlevels, etc. | X | Jonathan |
Meeting Update 03/23/2015 - Eddie
We need an abstract for AstroNum. I thought it would be useful to review our group's previous AstroNum talks, so here is what I found:
2011: Efficient Parallelization for AMR MHD Multiphysics Calculations: Implementation in AstroBEAR
- AstroBEAR intro
- Distributed tree, level threading, dynamic load balancing, sweep method, scaling result
2012: Multi-physics with AstroBEAR 2.0
- AstroBEAR intro
- Parallelization: threaded AMR and scaling results
- Multi-physics examples:
- interaction between magnetized clumps and radiative shocks
- anisotropic heat conduction
- magnetic towers
- disk formation around binaries
- effect of inhomogeneities within colliding flows
2015: Multi-physics with AstroBEAR 3.0?
- AstroBEAR intro
- new code stuff
- new multi-physics
- emphasis on my recent jet paper?
2-D Mach Stems
https://astrobear.pas.rochester.edu/trac/blog/ehansen03202015
Possible Runs for Stampede
- 3-D pulsed jets could easily use up the rest of the computational time. However, these are not ready to go yet. If I shift my focus, I could probably have this ready to run sometime in April.
- 3-D Mach stems. Right now, the runs take approximately 3 days each on 120 cores on bluehive which is fine. They are 40 cells/rclump. If we want higher resolution, Stampede could be a good option. We need to think more about exactly which runs we want for the paper, but these would be ready to go right away.
- An 80 cells/rclump run might use up about 100,000 SUs on Stampede? (rough estimate)
Intersection angles for all recent 2D Mach stem runs
As a reminder, the purpose of this continued 2D study is to find/verify the maximum intersection angle for Mach stem formation. My previous paper focused on the minimum angle.
Below is a table of the critical angles for a few cases of gamma.
gamma | Min Angle | Max Angle |
---|---|---|
5/3 | 37.38 | 41.81 |
1.4 | 41.65 | 58.05 |
1.2 | 47.67 | 69.73 |
1.01 | 68.73 | 85.92 |
I did a total of 11 runs with either gamma = 5/3 or 1.4. The initial intersection angle is controlled by changing the clump separation distance, and thus the minimum and maximum critical angles can be tested. Below is a useful table showing which runs exhibited Mach stems. The "Predicted" column is what we expect just based on the angle, and the "Actual" column is what we see from the simulations.
run | gamma | separation | Predicted angle | Actual angle | Predicted result | Actual result |
---|---|---|---|---|---|---|
A | 5/3 | 13 | 33.81 | 45.73 | regular reflection | tiny Mach stem, practically regular reflection |
B | 5/3 | 12 | 35.08 | 41.55 | regular reflection | stable Mach stem |
C | 5/3 | 11 | 36.59 | 42.78 | regular reflection | stable Mach stem |
D | 5/3 | 10 | 38.40 | 39.80 | Mach stem | stable Mach stem |
E | 1.4 | 6.5 | 38.97 | 31.75 | regular reflection | regular reflection |
F | 1.4 | 5.5 | 42.84 | 34.79 | Mach stem | stable Mach stem |
G | 1.4 | 5 | 45.16 | 32.38 | Mach stem | stable Mach stem |
H | 1.4 | 4.5 | 47.77 | 39.05 | Mach stem | transient Mach stem |
I | 1.4 | 4.25 | 49.20 | 31.81 | Mach stem | transient Mach stem |
J | 1.4 | 4 | 50.73 | 34.88 | Mach stem | transient Mach stem |
K | 1.4 | 3 | 57.87 | 39.66 | Mach stem | single bow |
The gamma = 5/3 runs:
The gamma = 1.4 runs:
Meeting Update 03/09/2015 - Eddie
I did some investigation into the intersection angles in my previous 2D Mach stem sims. I want to start thinking about the problem more time dependently.
Below is the gamma = 1.4, d = 5 simulation which results in a fairly wide Mach stem.
I think there is a good story here, and this is my interpretation:
- Bow shocks intersect.
- Bow shocks get "puffed" up due to increase in pressure in intersection region.
- The intersection angle crosses the minimum critical angle and a Mach stem forms.
- The Mach stem grows until the intersection angle hits the maximum critical angle.
- The Mach stem shrinks/grows a bit, almost like its breathing but it stays confined by the maximum critical angle.
Using this interpretation as a guide, I measured the relevant angles at different frames. I found that the angle right before Mach stem formation was about 42.75 deg (analytic theory says minimum should be 41.65 deg). Also, the intersection angle when the Mach stem stalls was about 58.33 deg (Pat's work says this maximum critical angle should be about 58.05 deg). This seems promising!
This seems like a good interpretation of the widest Mach stem scenario, but we need this to be consistent with the other models as well. There are 3 other possibilities:
- When the clumps are too close together, the Mach stem will grow beyond the maximum critical angle and disappear.
- For intermediate separations, the Mach stem does not always grow to the maximum angle. It can stop growing at some intermediate angle and be fairly stable. The d=5.5 model had a Mach stem that grew and stalled at an intersection angle of about 49.64 deg.
- For large separations, the intersection angle never gets above the minimum critical so a Mach stem cannot form.
Scenario 3 results in regular reflection, and we've seen that before. Scenarios 1 and 2 can be seen in the sims below:
3D Mach stems
I had to change the scale for the emission map with the M = 5 runs. There are some key differences here…
The Halpha is so weak that clearly the [S II] dominates. The Halpha maximum is only about twice as strong as the minimum, and the ambient is somewhere in between.
Meeting Update 02/24/2015 - Eddie
3-D Mach Stems
All 5 runs that I wanted to do for the 3-D Mach stem study are basically done. I still need to make a few more emission maps from different angles.
One issue that I'm running into is plotting the different Mach number runs with the same scale. The peak H-alpha emission is approximately 3 orders of magnitude stronger in the M = 10 cases compared to the M = 5 cases. See image/movie below.
I think I need to plot the two models with different scales, or maybe we should consider running different models…
One idea that I had was perhaps testing something other than velocity effects. We could leave M = 10, use the two different separations (d = 3, 6), and implement an initially constant magnetic field perpendicular to the flow. This would become a 3-D study on the effects of magnetic fields on Mach stems and intersecting bow shocks.
Perhaps this is too much, since we did not do any 2-D runs with magnetic fields. We could stick to simply studying more separation distances while leaving M constant at 10.
Below are emission and column density maps with 0 deg inclination for all M = 10 models. From left to right: single, d = 6, d = 3.
Other Stuff
- Last week we started discussing how to implement magnetic fields into the outflow module. Perhaps it would be good to have a separate meeting for this this week.
- I/We need to write an abstract for the AstroNum talk.
Update 02/03/2015 - Eddie
2D Mach stems
gamma = 1.4, d = 5 rclump results in wide Mach stem
3D Mach stems
Now have column density plots for d = 3, M = 10 case. I also changed how the angled emission maps are produced; inclination angles of 0 deg and 90 deg do not use a camera object. It will be tricky to look at the different angles side-by-side, but at least we can now be sure that any blacked-out regions are actually in the data and not being caused by a camera object.
Inclination | Density | Emission |
---|---|
0 | ![]() |
movie | |
30 | ![]() |
movie | |
60 | ![]() |
movie | |
90 | ![]() |
movie |
Other Stuff
Meeting Update 01/26/15 - Eddie
3D Mach stems
Ambient temperature lowered to 1250 K leads to a Mach number of 10 for the moving clump. Resulting emission is more complex.
2D Mach stems
For gamma = 1.4, still have not found a steady state Mach stem wider than the one in the d = 5.5 case. The "dimple" is more pronounced, but the reflection shocks go away.
Oblique Shocks
With the angle at 50 deg, this run should form a Mach stem.
Meeting Update 01/20/2015 - Eddie
3D Mach stems
Haven't done that much since last Thursday. See my most recent blog post: ehansen01152015.
I'm waiting on a new run that has Tamb = Twind = 1250 K. This will make the moving clump effectively at M = 10 as opposed to 5 in my previous runs. So we'll get to see how higher Mach number affects the evolution of the clump and its emission.
2D Mach stems
I have a couple of more runs for the maximum Mach stem formation angle study. If you recall, I'm trying to find the limit at which I get the widest possible Mach stem with the reflection shocks intact. For the gamma = 1.4 case, this appears to be somewhere between clump separation distances of 4 and 5.5 rclump. The new runs have d = 4.25 and 4.5. I will post those images/movies soon.
I also completed a single clump run, so that I can redo my bow shock shape calculations.
Oblique Shocks
I built an Oblique Shock problem module with the intent of studying Mach stems. I have a pretty good set up that generates regular reflection perfectly, but the Mach reflection runs still look a bit dicey.
Other Stuff
- Hope to finish my Qual brief this week
- In regards to the new machines that we want…I got Dave's email about the financial system being all messed up. I think we can hold off for another month or so on the larger, more expensive visualization machine, but I really think we should get Zhuo a new workstation asap despite the financial system's woes.
3D Mach stem update
These are the base parameters that I have used for previous simulations. The ambient temperature could be lowered if we want higher Mach numbers for the clumps.
Ambient/Wind | Clumps | |
---|---|---|
density (1/cc) | 1e3 | 5e5 |
velocity (km/s) | 27.2644 | 0, 10 |
temperature (K) | 5e3 | 10 |
So the density contrast is 500, and the clumps are initially in pressure equilibrium with the ambient.
We are in the reference frame of the slow clump, which is why the ambient is moving and the slow clump is given velocity = 0. Thus the slow clump has a Mach number of approximately 3.66, and the fast clump is at approximately M = 5.
The total simulation time ~ 50 years.
The clump separation is 3 clump radii.
I have completed 2 runs: the base run described above, and a run with twice the clump densities. Thus the higher density run has a nclump/nambient ratio of 1000. Below are the emission maps from the 2 runs at 4 different inclination angles: 0, 30, 60, and 90 degrees.
We do not have to continue with runs that have different density contrasts. Perhaps it would be better to study how clump separation and velocity affects the emission. I have a grid of runs below that we might try:
Clump Separation (rclump) | Fast Clump M |
---|---|
3 | 5 |
3 | 10 |
6 | 5 |
6 | 10 |
So this is just 4 runs to start. I would want to leave the slow clump alone, so it would remain at M = 3.66. The exact numbers can change, but I was basically thinking of doing what I already have + investigating a faster clump + investigating a clump that is farther away.
Also, do we want to move the clump so far away that it does not travel through the slow clump's wake? Or perhaps this is desirable, and we want to see an even closer approach?
Meeting Update 01/12/2015 - Eddie
2D Mach stems
See recent blog post: ehansen01092015
New results for gamma = 5/3:
I changed the density scale because I think this scale more clearly shows the Mach stems and post-shock regions. It appears that all 4 runs have Mach stems, but I only expected d = 12 and d = 11 to have them. It may be just because these separations are too close to the critical values. The d = 13 Mach stem shrinks and is very close to regular reflection which is expected. The d = 10 Mach stem grows and may be approaching the other limit where the bows merge. I'll have to run more models with d < 10 and d > 13 to really see the upper and lower limits.
I think it will be difficult to quantitatively prove or support Pat's calculations on maximum Mach stem formation angle. However, these simulations do clearly show, qualitatively, that a maximum angle does exist which explains why Mach stems grow and then disappear once the intersection goes above a certain angle.
It might be true that Mx approaches 1 at the triple point when a Mach stem just barely forms with the reflection shock still intact. This is at the maximum angle limit where the widest possible Mach stem forms. For this reason, I will be doing some more runs with gamma = 1.4 in hopes of finding and confirming this limit. So my runs will have separation distances between 4 and 5.5.
3D Mach stems
I completed a 3D run that will serve as the base comparison model. The emission map result is below:
This does look different than the previous run which can be seen at ehansen11252014. Fixing the initial ionization fraction changed things, but I don't know how much it is responsible for the differences.
Our initial plan was to explore the secondary bows that form behind the moving clump. I will investigate this by changing the clump velocity and/or the clump density. To start, I have a run going with twice the clump density.
This study will be somewhat slow moving because these runs each take about 2.5 days on 120 cores on bluehive. The files are also quite large (~500 GB per run), and disk space is limited on our local machines. For this reason, I will probably not save all of the chombos, so I will only show emission maps which come from the much smaller BOV files. Unless, anyone has a better idea to work around these issues?
Other Stuff
- I have made progress on some development things that I will share tomorrow.
- I started working on my Qualifying Brief, but I need to focus more on that for the next couple of weeks.
finished 2D Mach stem runs
Finished runs for gamma = 1.4 cases. See ehansen01052015 for more details.
The d = 5.5 case is the only one that I feel confident in saying forms a Mach stem, and the d = 6.5 case clearly does not. This is consistent with my previous work on minimum Mach stem formation angles.
Despite this, it would appear that the reflected shock (the shock that is reflected at the triple point) does not have a lateral Mach number > 1. This is unexpected, and I'm not sure what it means.
I expected to see a Mach stem in the d = 4 case, but it's hard to tell is this is one or not. One thing that we might need to consider here is what happens to a Mach stem if the reflected shock hits the clump.
Below are the same plots with slightly larger domain for d = 4, 3 since the bow shocks were at the edge of the previous images.
Maximum Angle for Mach Stem Formation
After last week's conversation with Pat, I started testing his calculations for the maximum Mach stem formation angle. My recent paper on Mach stems verified the minimum angle. If we look at both equations graphically, we can define a region where Mach stem formation is possible:
Below is a table of the critical angles for a few cases of gamma.
gamma | Min Angle | Max Angle |
---|---|---|
5/3 | 37.38 | 41.81 |
1.4 | 41.65 | 58.05 |
1.2 | 47.67 | 69.73 |
1.01 | 68.73 | 85.92 |
For my paper, I had graphically determined the shape of the bow shock for a few cases of gamma (could not do this for gamma = 1.01). This allowed me to convert minimum angles into maximum separation distances.
I used Pat's formula to determine minimum separation distances. If the clumps are closer than this minimum distance, the bow shock becomes a smooth continuous shape with no Mach stem. Below is a table of critical separation distances in units of clump radii.
gamma | Max Separation | Min Separation |
---|---|---|
5/3 | 12.47 | 10.79 |
1.4 | 6.20 | 3.62 |
1.2 | 4.58 | 1.81 |
Note that the gamma = 1.2 case has a minimum separation distance < 2 which cannot be tested. A separation < 2 means that the clumps are overlapping.
So I plan on testing the gamma = 5/3 and 1.4 cases. I have a set of 8 runs queued up as follows:
run | gamma | separation | Mach stem? |
---|---|---|---|
A | 5/3 | 13 | No, above maximum |
B | 5/3 | 12 | Yes |
C | 5/3 | 11 | Yes |
D | 5/3 | 10 | No, below minimum |
E | 1.4 | 6.5 | No, above maximum |
F | 1.4 | 5.5 | Yes |
G | 1.4 | 4 | Yes |
H | 1.4 | 3 | No, below minimum |
Below is what I have so far for run H. It seems consistent with what we were predicting, but we need the other runs to determine if Pat's formula is correct.
Pat also wanted to look at lateral Mach number plots. I generated a couple snapshots using two different scales:
new 3-clump runs
I implemented a rotation angle, so that we could view the emission maps from another angle. I also fixed the initial ionization fraction of the clumps, so they don't start off bright in H-alpha as they did before.
The emission maps below are both rotated 90 deg about the y-axis, so we are now looking along the x-axis. And they are also tilted, or inclined, 45 deg toward the observer.
The first one is the usual set up where all 3 clumps are in the same plane. The second one has the slower moving clump offset, so from this viewing angle you see this clump on the right.
Below is a slice of the planar set up showing density. This is just to get an idea of physically what is being simulated. The offset set up is the same except the slower moving clump on the left is rotated 45 deg about the stationary clump in a direction that moves it to the "background".
For the offset set up, you can't see all clumps via a slice, so 3D is required. I have to play around with 3D rendering a bit more, but hopefully this at least gives you a visual for the set up.
Meeting Update 12/15/2014 - Eddie
As a reminder, here are the parameters I used for my 3-D, 3-clump set up:
Ambient/Wind | Clumps | |
---|---|---|
density (1/cc) | 1e3 | 5e5 |
velocity (km/s) | 27.2644 | 0, 5, 10 |
temperature (K) | 5e3 | 10 |
So the density contrast is 500, the clumps are initially in pressure equilibrium with the ambient, and the fastest clump leads to a Mach number of approximately 5. Total simulation time ~ 50 years.
- Fixing the ionized sulfur species in the wind boundary condition solved the issues I was having with my emission maps. Here are some maps at different viewing angles:
- I ran a 1-D radiative shock model with parameters from my 3-clump set up and found a cooling length of approximately 4.999e13 cm = 3.342 AU. As we have done in the past, this is defined as the distance for the gas to cool down to 1000 K. This means that my cooling length resolution for the most recent runs is 13.37 cells/Lcool.
- This cooling length calculation is for the strongest possible shock (plane-parallel shock, velocity of fast-moving clump). Most of the shock structure on the grid will have cooling lengths longer than this due to oblique angles and slower moving or stationary clumps.
- At this resolution, the run took almost 3 days (69.38 hrs) on 120 cores on bluehive. This equals approx 8326 SUs. You might recall from last week, that this same run took 77.3 hrs on 96 cores (7420 SUs), so slightly slower but less expensive which is expected.
- I have updates on MaxSpeed tracking that I will share tomorrow, and then it will be time to focus my development efforts on the new Outflow object module.
Meeting Update 12/08/2014 - Eddie
- Got more file space on bluehive which will make post-processing emission maps much easier and quicker.
- Finished implementing MaxSpeed tracking. Need to analyze the results, will share at tomorrow's code meeting.
- Mach stem HEDLA paper revisions submitted.
- Discovered and fixed (hopefully) and error in my Mach stem module:
- A 3-D 3-clump run at 40 cells/rclump took ~77.3 hrs on 96 cores on bluehive. This is about 7420 SUs. Not bad at all, so we could reach higher resolutions if we need to. I'll figure out the cooling length as soon as I get a chance. The 3-clump sim that I had is re-running with the fix to the wind BC.
Meeting Update 11/25/2014 - Eddie
- ticket #431 (tracking density and pressure protections) is ready to merge.
- I tested it on one of my pulsed jet simulations and you can see it working.
- ticket #433 (tracking S species) is complete, and it will be ready to merge later tonight.
- I thought there might be an issue because I saw some strange behavior in the emission maps from my 3-D 2-clump Mach stem runs. However, I think the cooling length was just poorly resolved, because increasing the resolution got rid of the strangeness.
- ticket #436 (tracking max speeds) still needs some work.
- Work will be slow for the next few days because of the holiday, but I think I can finish this before next Monday.
Meeting Update 11/17/2014 - Eddie
- Finished 2.5-D pulsed jet paper revisions. The referee had asked about how a convolution with the instrumental psf would change my emission maps, so I wrote a program to calculate a convolution. Here are images of the beta = 1 model:
- The 3-D Mach stem simulation with cooling finished. Below is the density movie.
- Need to work on HEDLA paper revisions.
- Also working on setting up camera to make angled emission maps of my 3-D Mach stem runs.
Meeting Update 11/11/2014 - Eddie
- Had some issues with running 3-D Mach stem problem with cooling.
- It seems that the code can only use the exact Riemann solver with the ideal equation of state (and not a multi-species EOS which is what I use for Z cooling). The code immediately gets restarts due to nan in flux.
- There are 11 tracer fields: one for each clump, one for the wind, HI, HII, HeI, HeII, HeIII, SII, SIII, SIV. It's likely that the nan is being caused by one or more of these.
- For now, I have resubmitted the run to the queue using the default (HLLC) solver.
- On the afrank nodes, my 3-D Mach stem with cooling run failed after 7 frames. I got srun errors that I didn't really understand. I'm thinking it's a memory issue, so I switched back to the standard queue. Here is what I have so far:
- Altered my problem module to use 3 clumps now, and also implemented a camera object so that I can make emission maps based on inclination angle. This will be useful for comparing with observations. This just needs to be tested now.
Meeting Update 11/03/2014 - Eddie
- 3-D co-moving clumps, no cooling, gamma = 1.01, nclump/namb = 500, 40 cells/rclump, M = 5 shock for moving clump, difference in relative clump velocities = 10 km/s. Took 16.7 hrs on 96 cores on bluehive. Used default solver (HLLC).
- For 2.5-D pulsed jets, I am just going to regenerate figures with frame 96 and finish the paper revisions.
- Will also do paper revisions for HEDLA Mach stems paper this week.
- This is a busy week for PHY 101, but I should still have time to do development if we want to start having our weekly development sessions.
Density and Pressure Protections
In physics.data, there are options for protecting density and pressure. These control the Protect routines in EOS.f90. This routine is called various times throughout the simulation to make sure the density does not go below MinDensity and the temperature does not go below MinTemp.
There is also a routine called PrimProtect in EOS.f90. This routine is also called to protect density and pressure. The difference is that it is called on the q-array when the code is calculating predictor fluxes and final fluxes. Due to the way the code is written, PrimProtect cannot make use of the different protection options in physics.data.
Below, I have described how each routine handles density and pressure protections.
Protect
If the density < MinDensity, then a variety of things can happen depending on the options set in physics.data. Because momentum and density are related, they are both protected.
Momentum
- momentum is conserved and nothing happens
- momentum is set to zero
- momentum is set using average nearby velocity (default)
Density
- density is set to MinDensity
- density is set to the minimum nearby density
- density is set to average nearby density (default)
Then, no matter what options are chosen, the code sets the energy to zero for some reason.
Pressure
This protection occurs when the temperature < MinTemp.
- temperature set to MinTemp
- pressure set to minimum nearby pressure
- pressure set to average nearby pressure
- temperature set to minimum nearby temperature
- temperature set to average nearby temperature (default)
PrimProtect
Since data about neighboring cells is not available to the routines that call PrimProtect, the protection is more limited. As before, the momentum and density protection occurs when density < MinDensity and the pressure protection occurs when temperature < MinTemp.
Momentum
- momentum set to zero
Density
- density set to MinDensity
Then, the temperature is set to MinTemp.
Pressure
- temperature set to MinTemp
Meeting Update 10/27/2014 - Eddie
Last Week
- Finished (for now) the simulations of colliding Al plasma flows: ehansen10212014
- Did some nice little simulations of the NTSI for journal club: ehansen10232014
- Implemented protection tracking for silent protections. These are the protections that occur for the calculations of the predictor fluxes and final fluxes. Have some questions on why these protections do what they do.
- Got my MHD models of the 2.5-D pulsed jets to run. Beta = 5, beta = 1 completely finished, but beta = 0.4 only made it to 96/100 frames and is stuck.
This Week
- For the beta = 0.4 model, which is stuck at frame 96…
- It produces a new error that we haven't seen before in these simulations: "Src routine produced a NAN". The NAN is in the species tracers.
- It tries to reduce the time step and restart, but eventually shows "Min Timestep Reached, Stopping, …", and the code exits. I haven't seen this before either. Is this a relatively new feature for the purpose of stopping a run that is essentially stuck so that it doesn't burn SUs?
- I'm trying to get my final 4 frames, but I'm not sure if it will work out. It's not the end of the world if I can't get this to run. Worst case scenario: I just use frame 96 from all of my runs for comparison and for the figures for the paper.
- Once I have the beta = 0.4 issue resolved, I will have everything I need to complete the paper revisions.
- Working on implementing tracking of max speeds. Not just max speed of the velocities on the grid (sound speed and vx, vy, vz), but max speed of all velocities used internally within the code (e.g. the Riemann solvers).
- Next thing to do is the Mach stem paper revisions.
More colliding Al flows
Changed the set up a little bit so that there is no gap now (no ambient). The perturbation is now pseudo-random. The bottom half of the domain is initialized with either a density or velocity perturbation which contains multiple sinusoidal modes with random phases. Below is an image of an initial 5% amplitude density perturbation.
I ran 4 models: density and velocity perturbations, 5% and 10% amplitudes. All simulations were run to 300 ns.
Density Pert. | Velocity Pert. | |
---|---|---|
5% | ![]() | ![]() |
movie | movie | |
10% | ![]() | ![]() |
movie | movie |
I also made lineouts of density and temperature across the middle of the cooling region.
Density Pert. | Velocity Pert. | |
---|---|---|
5% | ![]() | ![]() |
movie | movie | |
10% | ![]() | ![]() |
movie | movie |
No cooling
For comparison, here is a 10% velocity perturbation run without cooling.
![]() | ![]() |
---|---|
movie | movie |
Inflow perturbed
One idea that we had was to keep the inflow from the bottom boundary perturbed with the same perturbation that is initialized in the bottom half of the domain. This drives the perturbation for a longer time and should lead to more dynamical effects. The perturbation is again 10% on velocity.
![]() | ![]() |
---|---|
movie | movie |
Meeting Update 10/20/2014 - Eddie
Wrote a program to compute convolution of emission map with Airy disk. This could be included in the code, probably in projections.f90, but I'll talk with Jonathan later about the best way to implement this if we want it.
Here is the Airy disk that I used for my convolution:
To see how this changed my emission map, here is an image of the original jet alongside the new emission map after the convolution. This process basically blurs the image a bit to more accurately represent what the telescope would see. In this case, this is what HST would see if my jet was located at Orion.
Back to Sulfur tracking results
Since implementing Sulfur tracking, I have only been able to successfully complete the hydro run for my 2.5-D pulsed jets. Below is an image showing the original emission map alongside the new one. As expected, allowing for SIII and SIV species means less SII and thus "greener" (H-alpha/[S II] is higher) emission.
Still having restart issues with the MHD runs. It looked like lowering the CFL might have helped at first, but the restarts still occur at later times. The restarts also occur at different times depending on how many cores are used. This leads me to believe that the real problem might be how the grid is distributed and how the AMR is refining.
Colliding Al Jets
Using typical lab parameters, I ran simulations of colliding Al jets. To simplify the problem, I ran in 2-D with the jet width = domain width. This essentially reduces to colliding two slabs.
The slabs are initialized with opposing velocities, and the bottom slab has an initial 10% sinusoidal velocity perturbation on its interface.
I ran two cases: the first is adiabatic, and the second has a special Al cooling table.
Adiabatic
As expected, the slabs collide and create a hot central shock that propagates back towards the boundaries. The central region remains hot and thus at a low density. The interface maintains an awkward shape due to the perturbation, but this shape is fairly constant.
Cooling
Note that the legends are much different. This is because I'm trying to highlight the regions of high density and low temperature. With this scaling , you don't see much in the density until about halfway through the simulation. Thus the cooling_rho.gif movie is only of the second half of the simulation.
The simulation starts out differently because cooling occurs immediately at the interfaces of the "jet heads". This quickly morphs into a similar configuration as the adiabatic case in which there is a hot central region with shocks propagating back towards the boundaries.
However, since cooling is now on, this central region eventually cools and forms dense regions. There is a periodicity to this because of the chosen perturbation.
There are times at which these high density regions look "clump-like". The key here is in the perturbation. The high velocity regions of the interface collide first and can thus start cooling sooner. Their shock strength is also slightly higher which leads to stronger cooling. This is where high density knots can form.
So unless the jet heads collide at the same time perfectly along the entire interface, you should always get high density knots within the post-shock region.
NTSI
I found that the perturbation that I was using was either too short-lived or not strong enough. I altered it a bit, and got some nice looking instabilities (namely, the Non-linear Thin Shell Instability).
For this simulation, I put a 15% velocity perturbation on both of the colliding flows. The average velocity and the continuously injected velocity from the top and bottom domain boundaries is 60 km/s. As a reminder, here are all the relevant parameters:
nslab = 1e-5 g/cc tslab = 10 eV vslab = 60 km/s (initially +/- perturbation) namb = 1e-8 g/cc tamb = 0.1 eV runtime = 100 ns domain size = 2 cm x 2 cm resolution = 64 x 64 cells + 3 AMR levels
The ambient is sufficiently low density and low temperature to play the role of a vacuum.
Note that in the following image/movie, the legend is constantly changing. This is because the cooling keeps increasing the maximum density and I wanted to always have the densest regions in red. The densities are in computation units (so in units of 1e-5 g/cc).
Once you enter the highly non-linear regime of the NTSI, you get filaments of high density scattered throughout the cooling region. Below is temperature in eV.
And just for comparison, here is the density image/movie for a simulation without cooling. Due to the shape of the perturbation, you get an interesting pattern, but you can see that it does not have the NTSI. Thus, so there is no high density fragmentation.
2-D test run for Mach stem project
I think I found an excellent regime for my Mach stem runs.
density contrast X = 500 shock strength M = 5 run time t = 25 years
With 2 clumps: one stationary and one moving at 10 km/s. For the image/movie below, I used gamma = 1.01 and a horizontal clump separation of 4 rclump (center-to-center). The effective resolution is 80 cells/rclump.
The run only took about 12 minutes on 64 cores. It will take longer to run with cooling and in 3-D, but it should be very doable at a decent resolution.
Code Meeting 10/07/2014 - Eddie
HLLC Solver
- The solver passed all 5 tests from Ch. 6 in Toro. The test parameters are given on my wiki page: u/ehansen/buildcode
- I tried adding transverse velocity (vy) to see if it would affect the solution to the Riemann problem. I used test 4 which is a shock/contact/shock test.
- I tried running my clump simulation with the clump and ambient/wind moving in the x-direction instead of the y-direction as before. This would test if the sweep method is causing asymmetry. If we find that the sweep method is the culprit, it would still be strange for the HLLC solver to produce more asymmetry than the Exact solver. Perhaps it is the combination of the sweep method with HLLC that is bad.
Pulsed Jet Restarts
- The hydro run finished, but it took almost 35 hours on 256 procs. I know that my original runs used more procs, so they ran in less time, but this still seems a bit high to me.
- I found 18 restarts due to High CFL, which seems higher than normal.
- How much slower should the code run because of the SII, SIII, and SIV tracers?
- For the MHD runs, there are a lot of restarts that will make it impossible to finish. I tried the beta=5 model, and I only got 20 out of 100 frames in 48 hours on 256 procs. The estimated wall time remaining got up to about 1 month.
- These are 'restarts due to nan in flux', which do not occur during the hydro run.
- I attached a list of restarts that show position, time, CFL, maxspeed,etc. It seems like most of them occur at the head of the jet, but it's a bit hard to tell without overplotting things, and I didn't have enough time to figure out how to to that.
Meeting Update 10/06/2014 - Eddie
- Implemented sulfur tracking: ehansen10042014
- new 2.5-D jet runs: hydro run has finished, working on MHD runs now
- testing 2-D clumps with different parameters for 3-D Mach Stem runs
- testing HLLC solver in 1-D
- using tests straight from Toro which I used 3 years ago on my own hydro code: u/ehansen/buildcode
Tracking Ionized Sulfur Species
See ticket #433 for a brief description.
I found recombination rates from Nahar 1996. There is a nice table with many more species that could potentially be used in the future.
Ionization rates were a bit more complicated. I used Arnaud & Rothenflug 1985. They give some complicated formulas for the rates of several species. Again, the other species could be implemented in the future. I wrote a separate fortran program which uses these formulas to generate a table of ionization rates. The generated table is in a format that is compatible with astrobear's table and chemistry routines.
I tested my program with H ionization, and I was able to recover the table that we currently use in astrobear. This program is so useful that I'm going to check it into the main repo so that we can use it in the future for other species.
Here are the ionization rates for SII and SIII. The x-axis is log(T) and the y-axis is the rate in units of 10-9cm3s-1.
These follow a trend that is common to ionization rates, and it makes sense that SIII peaks at a higher T than SII. The peaks don't seem that far apart because the ionization potentials (23.4 and 35 eV) are pretty close in log(T) space. Below are the recombination rates for SIV and SIII. Same units as above except the rates are in log space now.
These were straight from a table, no calculations necessary which is quite convenient. Makes sense that SIV recombination rates are a bit higher since it should be harder to ionize SIIII / easier to recombine SIV. This is consistent with the ionization rates as well.
Jets
I tested the implementation with a low resolution, non-MHD, pulsed jet. Below are plots of SII, SIII, and SIV density (from left to right in computational units).
The plot makes sense qualitatively. SII peaks inside the clumps, aka cooling regions. SIII peaks at the edges of clumps and on the wings, aka shocks. There seems to only be SIV on the spur shocks where it is hot enough. I did find a typo in the SIV recombination table which was extremely high at one data point, so there's a chance that the SIV density should be slightly higher but I think it still looks good. We didn't expect much SIV anyways.
Despite that faulty data point, I think it is fair to say that sulfur tracking is working properly. The emission maps should not have changed much, since all I had to do was multiply the output by SII/Stot where Stot = SII + SIII + SIV. The emission maps really need higher resolution to be accurate, so there's no sense in looking at the ones from this run. I am going to set up production runs now.
Meeting Update 09/29/2014 - Eddie
- Did some interesting testing on velocity/momentum smoothing for clumps: ehansen09262014. What do we want to be the "correct" way? Should I remove velocity smoothing in our main development branch?
- Tested the 3-D clump set up (see below)
- Plan on doing 2.5-D pulsed jet paper revisions tonight.
3-D Clump
- single clump (density contrast = 2000)
- no cooling, gamma = 5/3
- ambient/wind and clump in motion (M ~ 10)
- velocity smoothing is still on
- uses exact Riemann solver, linear interpolation, CFL = 0.5, 0.2, 0.5
- effective resolution = 40 cells/rclump
Below is a 2-D slice of the density:
Still playing around with 3-D rendering. It takes a long time on my machine, so I may try using interactive jobs on bluehive from now on.
The above simulation completed 28.5/100 frames in 8 hours on 64 procs on Stampede. The remaining wall time estimates lead to a total run time of about 2 days for the entire 100 frames. On 64 procs, this is about 3000 SUs. Using Z cooling will at least double this to 6000. Adding a 2nd clump will generate more gradients and significantly increase run time (I'm guessing by a factor of at least 10). Best case scenario is that these runs take about 60,000 SUs at this low resolution of 40 cells/rclump.
I should also note that the above numbers are for a simulation time of 32.87 years, which seems reasonable when compared to observational time scales. Pat's HST observations show changes in emission features on even shorter time scales (~10 years). We may be able to improve resolution and run times if we limit our simulation parameters in such a way to make it easier on the code.
I spoke with Shule, and he said that his shocked clump simulations were at 64 cells/rclump, MHD, and took up to 2 days on 1024 procs on Kraken. So about 48,000 SUs. This seems much better than what I'm getting. It's hard to compare hydro + Z cooling vs. MHD + no cooling, but maybe if I tweak my parameters I can get up to this resolution.
Clump Behavior: Velocity/Momentum Smoothing
After comparing the last set of simulations with moving vs. stationary clumps, I decided to investigate the role of velocity/momentum smoothing at the clump boundary. I found a small possible error in the code. Momentum gets smoothed by two factors of f where f is a tanh function. This happens because both density and velocity are smoothed by a factor of f. I think that it might be better to only smooth density, since smoothing velocity would make the code more dependent on reference frame. I took out the extra velocity smoothing, but I am not sure what affect this will have if any.
I looked at the initial frames for a few cases in velocity and momentum. There is the original smoothing, the corrected smoothing, and no smoothing. Below are the images for momentum.
Stationary (orig) | Moving (orig) |
---|---|
![]() | ![]() |
Moving (correct) | Moving (no smoothing) |
![]() | ![]() |
Below are the movies showing density.
Stationary (orig) | Moving (orig) |
---|---|
![]() | ![]() |
movie | movie |
Moving (correct) | Moving (no smoothing) |
![]() | ![]() |
movie | movie |
Some strange stuff happens without any momentum smoothing. This is most likely due to the fact that density is still smoothed, so there is a mismatch between density and momentum.
The corrected smoothing case looks a bit more like the stationary case which is good. However, there is velocity smoothing even in the original stationary model. So I should really do a new stationary run without the extra factor in velocity smoothing, and compare that to the moving clump also without that extra smoothing factor. We want these two simulations to look the same since the code should be Galilean invariant. Below are the final results with corrected smoothing:
Stationary (correct) | Moving (correct) |
---|---|
![]() | ![]() |
movie | movie |
So taking out the extra velocity smoothing didn't affect the stationary model as much as it affected the moving model. However, I would still argue that the correction is an improvement. I don't think we should expect the stationary and moving models to be exactly the same, since you will get a different amount of numerical diffusion when you have things moving on the grid differently. Stay tuned for 3D runs.
Meeting Update 09/22/2014 - Eddie
Clump Behavior
Before moving on to 3-D simulations, I wanted to see how setting the clump in motion affected the clump behavior. I set the clump velocity at 10 km/s, and reduced the wind by the same amount to keep the relative velocity the same. I tested both the low density contrast (X = 500) and high density contrast (X = 2000) set ups.
There was also a small error in my previous simulations. I forgot to change the ambient temperature to 5000 K as I had stated. The temperature was actually at 1000 K, which means that the shock was not Mach 10, it was more like Mach 22. This doesn't change what we found out about the different solvers and what not. The new set of runs below are with the correct ambient temperature of 5000 K =⇒ M = 10 shock.
Low Contrast | High Contrast | |
---|---|---|
Stationary | ![]() | ![]() |
movie | movie | |
Moving | ![]() | ![]() |
movie | movie |
These were all done with 2nd order interpolation and the exact Riemann solver, which still seems like the way to go.
Other Stuff
- The 3-D version of the high contrast, moving clump simulation above is currently running. If that looks good, I can think about using a 2nd clump again, and Z cooling as well.
- Need to redo the 2.5-D pulsed jets emission maps at higher resolution for the film thingy on Friday
- Working on 2.5-D pulsed jets paper revisions
Clump Behavior: Alexei's Setup + High Density Contrast
The same setup as before (ehansen09162014), except now the clump density has been increased by a factor of 4 to make the density contrast 2000. This time, I only tested the 3 different Riemann solvers with 2nd order accurate interpolation. Also, the CFL vars have been lowered to 0.5, 0.2, 0.5 for all runs.
Exact | HLL | HLLC |
---|---|---|
![]() | ![]() | ![]() |
movie | movie | movie |
Some interesting facts about run times…These all ran on 64 procs on Stampede.
Solver | Run Time (mins) |
---|---|
Exact | 66.88 |
HLL | 180 (to frame 92.5) |
HLLC | 68.34 |
Some Observations
- HLLC is still very asymmetric
- Exact has become more asymmetric
- HLL is now most symmetric but is also most diffusive and for some reason the code is much slower with HLL
So Exact still looks like the best option even though a high density contrast has led to more asymmetry. Is it possible that we are pushing the accuracy limits of the code with such a high density contrast?
I think the next logical step is to put the contrast back down to 500 and put the clump in motion. Then we will know how clump motion affects the simulation independent of the high density contrast.
Clump Behavior w/ Alexei's Setup
Going back to basics, with a simulation setup that we know should work. Here are the parameters:
nwind = 1000 cm-3 Twind = 5000 K vwind = 74.53 km/s ==> M = 10 nclump = 500,000 cm-3 ==> X = 500 Tclump = 10 K (pressure equilibrium) vclump = 0
I tested different Riemann solvers and different interpolation orders. I was only able to test the Sweep method, because MUSCL is only implemented for fixed grid and 1 proc.
Also, the HLL linear run had the most trouble. It only got through 61/100 frames in 2 hours on 64 procs. This may have to do with whatever happened around frame 45…looks like a density protection maybe. All the other runs got through the full 100, and usually in much less than 2 hours with 64 procs.
Exact | HLL | HLLC | |
---|---|---|---|
1st Order (Godunov) | ![]() | ![]() | ![]() |
movie | movie | movie | |
2nd Order (linear) | ![]() | ![]() | ![]() |
movie | movie | movie |
I also did a run with the 'default' solver settings (sweep, HLLC, linear), but with a lower CFL. The original CFL variables were 1, .5, .5, and the run below has .5, .2, .5.
Some Observations
- All the Godunov runs are more diffusive which is expected.
- Out of all 3 Godunov runs, only the HLLC solver produces asymmetries.
- All of the 2nd order runs become asymmetric at some point with perhaps the Exact solver being most symmetric and HLLC most asymmetric.
- Lowering the CFL did not seem to make a significant difference.
Meeting Update 09/15/2014 - Eddie
- A decent chunk of my time is being spent on diagnosing the strange behavior in my clump simulations.
- Checked AMR vs. fixed, and cooling: ehansen09092014
- Tried 2.5-D: ehansen09112014
- Looked back at Kris' shocked clump paper, and I am working on calculating different time scales. I'm curious to see if there are any instability time scales, or perhaps the cooling time scale, that are much shorter than the clump destruction time scale.
- Ran some simulations to check for convergence in my 2.5D pulsed jet simulations: ehansen09122014. Just need to come up with a better metric to show convergence.
- Implemented tracking for pressure and density protections. See image below where the cells that have triggered protections are marked in black. I will post another blog post later to further explain this new feature.
Convergence of 2.5D Pulsed Jets
This is related to one of the comments by the reviewer of my 2.5D pulsed jets paper. Basically, he wanted to know if our simulations had converged to a solution. So I did a quick resolution study. From left to right is increasing levels of AMR.
The resolutions are so high for all of these simulations that you can't really see a difference here. In other words, there are more computational cells than there are pixels to display. I'll post an enlarged image later of one of the clumps.
UPDATE
Below is the clump located near z = 19 Rj from the above image.
2.5D Single Clump Tests
Since the instabilities at the head of my 2-D clumps are not going away, we decided to try 2.5-D. I ran two cases, one where the clump is still in the center of the grid (like simulating the cross section of a donut), and one where the clump is centered on the z-axis to take advantage of the imposed symmetry (like simulating the cross section of a ball).
I'm beginning to think that these instabilities may be physical. If I'm right then perhaps the reason I see them in my simulations when others have not is due to the large density contrast between clump and wind. Stay tuned as I figure out if there's a reasonable explanation for the instabilities.
Strangeness in Clump Simulations
In order to figure out what is going on with these clumps, I am simplifying the problem to 1 moving clump. We will also be looking at a clump that is initialized in a moving ambient.
Here is the run with all of the original parameters. This uses Z cooling, 80 cells per rclump (4 levels of AMR), Hvisc, and lapplydiffusion are on. Additionally, interporder = 2, lcharlimiters is on, lrestartonpressureprotections and lrestartondensityprotections are both off, and cfl_vars = 1, 0.5, 0.5.
Test Runs
Here are all the models we talked about yesterday. They all have an effective/actual resolution of 80 cells/rclump.
- Zcool/nocool, AMR/fixed models are self-explanatory
- lower cooling = all densities decreased by factor of 10
- lower contrast = clump density decreased by factor of 50
- increase diffusion = increased DIFF_ALPHA in solver.data from 0.1 to 0.2
Model | Zcool, AMR (orig) | Zcool, fixed | nocool, AMR | nocool, fixed |
---|---|---|---|---|
Image | ![]() | ![]() | ![]() | ![]() |
Movie | movie | movie | movie | movie |
Model | Zcool, AMR (orig) | lower cooling | lower contrast | increase diffusion |
---|---|---|---|---|
Image | ![]() | ![]() | ![]() | ![]() |
Movie | movie | movie | movie | movie |
My thoughts:
- Zcool, fixed still shows the weird stuff at the head of the clump. It starts off symmetric, but then develops asymmetries which is odd. Thus, it would appear that the AMR is not seeding any instabilities.
- nocool, AMR/fixed both look better but not perfect. So although much of the problem may be caused by the cooling, I may still need to add some diffusion to get rid of the problem entirely.
- lowering the cooling by lowering the density everywhere helped. I suspect that not resolving the cooling length may directly lead to the instabilities, or it causes pressure protections that lead to the instabilities. One thing to try is to turn on restarts for pressure protections.
- lowering the contrast by lowering the clump density did not help. The clump was destroyed much more easily, so it's hard to draw any useful conclusions from this run.
- increasing the diffusion helped a lot. This run looks nice and smooth. We just have to ask ourselves at what point does adding more diffusion become nonphysical.
UPDATE
Tried turning on restarts for pressure protections. Restarts are triggered almost immediately and often. This will prevent the run from ever finishing, so I had to stop it. I don't see any nans, so it is the protections are presumably caused by very strong cooling that brings the temperature below the mintemp of 1d-10. Perhaps increased resolution could help, or lowering the mintemp and/or raising the floortemp.
They are being triggered at the point (-0.41, 9.12) which is near the lower left edge of the clump. Maybe there is something with the softening or refinement buffering that might help.
Meeting Update 09/08/2014 - Eddie
Not too much to report since last Thursday. I've been really busy setting up the PHY 101 course. The organization for this course has been a complete mess, but all the confusion should clear up this week which will free up more time for me to get back to research.
Mach Stems
I have a few runs queued up on Stampede to test a moving clump through an ambient. Not sure why they haven't started yet. Perhaps there is something going on with Stampede that I'm unaware of. ehansen09042014
This Z cooling run took about 90 minutes on 256 processors on Stampede (80 cells per clump radius). A no-cooling run took about 50 minutes.
The other test runs will be testing parameters related to diffusion.
Intersecting Shocks
I am going to work on the theory behind shock interactions today. I need to use shock jump conditions to figure out parameters for a post shock region behind two intersecting oblique shocks.
New Outflow Object
I am also testing the new outflow object. The diverging wind option quit with an error that I am working on debugging. I can expand on this topic tomorrow.
Mach stems w/ moving clumps
I'm in the process of setting this up. Here is an example image/movie of what the simulation is:
So the functionality is there, but I need to get the spacing and timing right so that we get a satisfactory simulation to look at and analyze.
Update
The above simulation was using Z cooling. The one below is without cooling…
Meeting Update 08/25/2014 - Eddie
New Outflow Object
The new outflow object module is starting to come together. I successfully ran Martin's problem module with an outflow object of a simple 2.5-D hydro jet.
Next, I will test some of the velocity parameters and the diverging wind. There are some things that were already built into the outflow object module that I have excluded for now. I will have to re-implement those later.
Other Stuff
- Need to revise the 2.5-D pulsed jet paper
- Will start 2-D Mach stem runs of clumps with relative velocities
Meeting Update 08/18/14 - Eddie
Basically just working on two things right now:
- the jet/outflow object module (ticket #370)
- this is more of a code issue that I would like to discuss briefly tomorrow
- revisions to my 2.5-D pulsed jets paper
Science Meeting Update 07/20/2014 - Eddie
For reference: ClumpClump/Morphology
Mainly just working on the Mach stem paper for the HEDLA proceedings. I hope to have it done within the next few days. I did some analysis on the shapes of the bow shocks from my simulations to confirm that they agree with the theory about Mach stem formation…
From analyzing the density maps, I was able to get a functional form for the bow shocks which end up looking like this:
Then, I can figure out what angle they would make if they intersected another bow shock. The theory of Mach stem formation says that they form at and above a critical angle which depends on gamma. For the values of gamma that I have used, these angles are:
gamma | phi crit |
---|---|
5/3 | 74.75 |
1.4 | 83.31 |
1.2 | 95.34 |
1.01 | 137.46 |
Now, going back to the bow shock shape, I was able to pinpoint where these critical angles would occur. In other words, I figured out the critical separation distance for the 2 clumps at which Mach stems should begin to form. These are approximately:
gamma | d crit |
---|---|
5/3 | 12.5 |
1.4 | 6.2 |
1.2 | 4.6 |
1.01 | ? |
I can't figure out d crit for the gamma = 1.01 case because the bow shock is unstable. The simulations are consistent with the theory because of which simulations we see Mach stems:
gamma | d | Mach stem? |
---|---|---|
5/3 | 15 | no |
5/3 | 10 | yes |
5/3 | 5 | yes |
1.4 | 10 | no |
1.4 | 5 | yes |
1.2 | 5 | no |
All simulated models not listed in the table above did not show evidence of a Mach stem.
Intersecting Bow Shocks: Extended Boundaries
This is the gamma = 5/3, d = 15 model. I extended the left and right boundaries to be 20 clump radii away from each clump, and this was enough to ensure that the bow shocks do not hit those boundaries.
I am in the process of making side by side comparisons of all the other models. I will post them to this page: ClumpClump/Morphology
Makefile.inc for darter
I have these lines in my .bashrc:
module purge
module load PrgEnv-intel
module load fftw mercurial cray-hdf5 cray-mpich
ulimit -s unlimited
Clump and Bow Shock
Ambient
- density = 103 1/cc
- temperature = 104 K
Clump
- density = 107 1/cc
- temperature = 103 K
Wind
- density = 5 x 103 1/cc
- temperature = 104 K
- velocity = 50 km/s
Also, gamma = 1.01, and the effective resolution is 160 cells/rclump.
The next image/movie is the same run, but with a different color bar.
Fixed Grid
No AMR, still 160 cells/rclump. I made 2 different movies of the same run, both are of density. The only difference is the color bar.
Science Meeting Update 06/09/2014 - Eddie
- Last week I mainly worked on fixing my 2.5D pulsed jet simulations, and regenerating the figures for the paper.
- I have also been working on the Al colliding flows problem with Francisco. The most recent simulation looks better than this one from last week: http://astrobear.pas.rochester.edu/trac/blog/ehansen06072014.
- The small scale structures and explosion went away when I imposed a static mesh.
- However, there is still an odd asymmetry.
- I am also working on the Mach stems project. I need to get some good runs completed so I can write the HEDLA proceedings. Right now, I am simply running a wind over a clump to see if I can get a nice bow shock. My previous simulations looked somewhat strange. This is queued up on bluestreak.
Al Colliding Flows: Adiabatic vs. Cooling
The adiabatic case is not very interesting, so I only posted the density movie. The flows collide and form reflection shocks which propagate all the way to the boundaries.
I ran the cooling case for longer time. Now each frame is equal to 10 ns. There is a lot going on in this simulation. The small structures that form due to cooling look nice and are believable. However, there is a strange explosion near the end of the simulation. My intuition tells me that this is not physical, but we should discuss this further and figure out why the code produced this result.
There is some asymmetry between the top and bottom halves of the grid which is also curious. The inflow at both of the boundaries are the same. It looks like the bottom half of the grid cools more than the top half.
I also included an enlarged image of the first frame of the explosion. The image shows pressure, temperature, velocity, and density. This might have been caused by AMR effects, so I will try a run that has a static fine mesh in the collision region.
density maps for new 2.5D pulsed jet sims
Resolution is twice as high as the old runs, and the cooling is stronger because of the error that I fixed in the temperature calculation.
It is interesting that the hydro and beta = 5 cases are so different. You might expect the runs to approach the hydro case as you increase beta, but this is not what we see here.
The runs are (from left to right): beta = inf (hydro), 5, 1, 0.4.
New Colliding Al Plasma Run
This is the same setup as before: http://astrobear.pas.rochester.edu/trac/blog/ehansen05302014.
The only difference is this is with an updated cooling table, so I had to change the scaling on the legends. I also turned on h_viscosity to help get rid of any nonphysical carbuncles.
Science Meeting Update 06/02/2014 - Eddie
- High resolution 2.5-D pulsed jet sims are finished. I am now just working on regenerating the figures for the paper. Below is an image and movie comparing the old hydro run to the new hydro run.
- I'm having trouble running the 2.5-D random velocity model on stampede. It always quits before it starts time-stepping, and I'm not sure why. It should only take a day or two to run once I figure out the problem.
- Colliding flows simulations are working: http://astrobear.pas.rochester.edu/trac/blog/ehansen05302014
- Francisco sent me an updated cooling table this morning which I will implement, and then I will rerun the same setup.
Colliding Alumin(i)um Plasmas
The point of these simulations is to study cooling instabilities that may be present in Francisco's experiments of colliding jets. To simplify the problem, the simulations are just 2-D colliding flows with the appropriate Al cooling. The initial parameters are as follows:
T = 10 eV
rho = 10 ug/cc
v = 60 km/s
Z = 6 (ne/ni the electron to ion number density ratio)
The grid is 2 x 2 cm and 64 x 64 cells with 6 levels of AMR. The total runtime is 100 ns.
In the experiments, the "ambient" is a vacuum. To get around this in astrobear, I initialized the flows such that they already fill the domain. So the flows immediately collide at the beginning of the simulation, and they continue to flow in from the top and bottom boundaries. The left and right boundary conditions are periodic.
Below are images and movies of the temperature in eV, and also the density in ug/cc. They are enlarged to look at the region of interest.
Meeting Update 05/27/2014 - Eddie
Science
- The new high resolution 2.5-D pulsed jets (5 models) are running on bluehive2 and stampede.
- Wrote a problem module to study cooling instabilities in Francisco's colliding jets experiments. Now working on implementing an aluminum cooling table.
- Need to analyze Mach stem simulation data, and possibly do some more runs to write HEDLA proceedings.
Code
- Made a quick fix to the code (ticket #401)
- Need to look at Martin's jet/outflow module and combine my work with his.
- Julio should have the code now. I will give him some time to get set up and play around with it, and then I will contact him at the end of the week to see where he is at.
new resolution numbers for 2.5D pulsed jet paper
Now that I fixed that typo in the cooling routines, the non-equilibrium cooling will be a bit stronger. This is good news for the 2.5-D sims as long as I am still able to resolve the cooling length, and it will make my resolution more impressive when compared to other papers. I redid the cooling length calculations and updated the following table in my paper:
This will, however, be bad news for the 3-D sims. We've already seen that it is very difficult to get high resolution on these runs. The shorter cooling length means we will be that much less resolved. We may have to reconsider the problem set up, and do something that the code can handle more easily.
Science Meeting Update 04/28/2014 - Eddie
- Found an error in the code that has been affecting cooling. Here is a plot that shows the temperature profile for a shock model before and after the error was fixed.
- This fix will affect any simulations that I have done recently or in the past that use cooling. This includes the pulsed jets and mach stems projects that I've been working on.
- I reran the 2.5-D hydro pulsed jet model that I used in my paper and compared the new results to the old. The old model is on the left, and the new model is on the right…
The overall morphology doesn't look too different. As expected, the big differences are in the small scale structures and the head of the jet.
- The high resolution 2-D mach stem runs are done, but again the Z cooling models are with the old version of the code which has the aforementioned error. I suspect that rerunning these would make the Z cooling models look even closer to the gamma = 1.01 models. These high resolution Z cooling sims took a very long time on bluestreak, but I'm not really in any hurry so I should be ok to rerun.
Science Meeting Update 03/31/2014 - Eddie
- Working on high resolution 2-D mach stem runs. Ruka has 16 of them, and I have 4. The 4 that I have are zcooling, so they take a bit longer (3 submissions per run on bluestreak). They will probably be completed sometime next week, but it depends on the queue times.
- The 3-D pulsed jet run is coming along, although it seems to be slowing down quite a bit. I have noticed that the slow down is not primarily due to restarts. The last 24 hr run that went through generated only 3 frames, and there was only 1 restart within that run time. If I can only get a few frames every 24 hrs, then it might be time to switch tactics as I would need about 500,000 SUs to complete the simulation. Here's the density image/movie of what I have so far:
- Need to look back on some old stuff about zcooling and my 1-D radiative shocks, so that I can write an informative email to John Raymond.
- Will try the emission analysis for my 2.5-D pulsed jets on the OI line like Pat suggested. I built that line into the code a long time ago, so it's just a matter of turning a switch on and rerunning astrobear in post-processing mode which should not take very long.
- I am also working on updating all of the problem modules for the new test suite.
Code Meeting Update 03/25/2014 - Eddie
Pulsed Jet Restarts
My latest run is just finishing on Kraken. It made it to 30 frames in 24 hours on 960 cores. There were 10 restarts triggered, 9 of which were high CFLs. Only one of them was the restart due to nan in flux. It looks like one of the toughest points was going from frame 24 to 25 which had 2 restarts and took almost 4 hours, but it eventually made it through.
I expect the frame output rate to slow down over time since the jet propagates into the grid and triggers more refinement. But the restarts make the run go even slower. It is possible that I could get this run to finish as is, given enough cpu hours, but I'm not sure how much we want to invest.
Right now, a good frame step that has no restarts, takes about 1.5 hrs. In the beginning, they only took about 10 minutes. If I assume that this trend will continue linearly, then at best the rest of the simulation would take 220 hours which is a little more than 9 days. With this many cores, this is 211,200 SUs. This is a best case scenario estimate. These numbers will increase due to triggered restarts, the need to physically restart the code every 24 hours will lead to inefficiency, and also the frame output rate probably decreases faster than linearly.
I have not yet had a chance to analyze the chombos from this run, but I attached the standard output to this post.
Czarships
Testing
We are still in the process of completing the new test suite which no longer uses bear2fix. I believe all the routines are in place, and now it is just a matter of setting up each of the problem modules to be compatible with these new routines.
Local Resources
Visit has been working very slowly on my machine (Grass). If I work on clover, ssh to grass, and then run visit, it is much faster. So I'm guessing that the bottleneck is on Grass' graphics card, and thus it has difficulty displaying things locally.
I don't know if this is something Rich can fix, but I think it might also be time for a new machine. Can we get the ball rolling on figuring out what we want and purchasing it?
Science Meeting Update 03/24/2014 - Eddie
- I have Ruka helping me do some mach stem runs. Specifically, we are working on the 2-D, high resolution runs (160 cells/rclump). When those are finished, we will move on to 3-D.
- I again worked on the emission map analysis for my 2.5-D pulsed jets, but the results did not look very promising. It seems that the FWHM calculation does not support our interpretation of the morphology of the jet and its internal working surfaces. We might be able to do something with total emission.
- I have a new pulsed jet run in the queue on kraken. This is 3-D, beta = 1, with random velocity pulses.
- There is still some work that needs to be done on the test suite which I didn't get to last week. I will focus more attention here and get this done soon.
Science Meeting Update 03/18/2014 - Eddie
I've been working on analyzing the emission maps from my 2.5D pulsed jets, but I'm not having much luck. The emission does not seem to be as well-behaved as we had hoped. Below is an image of the Halpha for one of the internal working surfaces:
The Halpha seems to peak outside the jet radius instead of at the jet center. This doesn't directly contradict what we were saying in the paper, but it makes the analysis that we want to do impossible. I'm referring to the full width half max idea.
So I started thinking that maybe there could be an error in the emission map projection that Jonathan and I implemented. Unfortunately, the chombos do not contain any emission data, so I do not have the raw data from before it gets projected. The best I could do for now is to look at densities and temperature.
Neutral H:
Temperature:
Based on the H density, you might think that Halpha emission peaks at the center, but the temperature says the Halpha emission should peak outside the jet radius. So it's very hard to tell what the emission is supposed to look like.
UPDATE
I think I figured out how to get the raw data for the emission. Here is what I got for Halpha:
Science Meeting Update 03/03/2014 - Eddie
- I got a low resolution 3-D mach stem simulation to run on bluestreak, and now I'm running a higher resolution set up with periodic boundary conditions. If this works fine, then I'll move to kraken and do my production runs.
- I ran a 2.5-D pulsed jet with random velocity pulses. It looks more like real observations because you get more shocks and clumps close to the source. The density images are being made right now, and I'll do the emission images later.
- Following Pat's comments, I've gone through the abstract and intro of my paper. I still have a lot of work to do on the figures and emission line analysis.
- I'm reading through a couple of the papers that Pat pointed me to last week. Namely the 2007 paper that has the random velocity pulses, and the 1993 paper that he did with Raymond which was on emission from 1-D models of pulsed flows.
Testing Report 02/252014
Jonathan, Baowei, and I have been working on updating the test suite. Baowei has been focusing on the buildproblem script, repository, and wiki implementation while Jonathan and I have been working on creating new testing routines. The goal is to no longer use bear2fix and give astrobear its own testing capabilities. Here's how it will work:
- We will now have a test object that handles error calculations and image production. I have created a wiki page that details how to use this object: TestObjects
- Each problem module in the test suite needs to have a few things
- a file named 'layout0000.dat' which contains verified data that new data gets compared to (this is instead of a reference chombo)
- a parameter 'ltest' in globabl.data (default is F) so that the buildproblem script can turn on/off testing
- in problem.f90, the test object needs to be created and certain parameters may need to be initialized.
- there should also be a reference image 'refimage.png' that we can use for visual verification
Science Meeting Update 02/24/2014 - Eddie
- working on polishing the 2.5D pulsed jets paper
- I was testing a 3D mach stem set up, and got it to run on bluestreak, but bluestreak is down for a couple days for maintenance so I can't access the data yet.
- I want to start reading up on jet launching. Christian Fendt, who I met in Heidelberg, has done some work on this so I'm looking at some of his papers. Here's one that focuses on the effects of magnetic diffusivity: http://adsabs.harvard.edu/abs/2002A%26A...395.1045F
Testing Report 02/18/2014
- We only do weekly testing on bluehive for the official revision. Our latest official revision of AstroBEAR2.0 is 1245:6ea1a968486d, which was checked in on Aug 19th. So it's basically testing the old code over and over. We need to start testing our development version of AstroBEAR3.0.
- The current testing suite revision is 161 — for AstroBEAR2.0, which was checked in On July 23rd. This was copied to botwin when we moved the wiki from clover to botwin
- The testing script buildproblem needs to be updated: a) The testing results are saved to clover; b) but we lost access to /cloverdata/ from other local machines possible due to failure of one disk on clover.
- The personal current-tests page which you can check the testing result on different machines by [[CurrentTests( username_machine)]] need to be updated, as again the results are saved in clover and we lost direct access to the machine
- Baowei working on updating the testing script and wiki page as he needs to test Jonathan's fix to Marvin's ticket.
- Want to go over the modules that are currently tested, and make sure they are up to date. Also need to know which new modules we will have in 3.0, and which of those should be tested regularly.
Current Problem Modules Tested in 2.0
- BasicDisk
- BonnorEbertSphere
- Bondi
- BrioWuShockTube
- CorotatingBinary
- EinfeldtRarefaction
- FieldLoopAdvection
- FieldLoopRestart
- FillingFraction
- GravitationalCascade
- HydroWaves
- IsoHydroWaves
- IsoMHDWaves
- IsotropicTurbulence
- Marquee
- MHDWaves
- MolecularCloudFormation
- MomentumConservation
- MultiClumps
- OrbitingParticles
- RadiativeInstability0
- RadiativeInstability05
- RadiativeInstability1
- RadShock
- RTInstability
- SingleClump
- SlowMolecularCloudFormation
- SodShockTube
- ThermalInstability
- UniformCollapse
2.0 Problem Modules Not Tested
Module | Should it be tested? | |
---|---|---|
1DWaves | ||
2DWaves | ||
Binary | ||
Christina_Original | ||
CurrentSheet | ||
GaussDiffusion | ||
GravoTurbulence | ||
HydroStaticStar | ||
jets | no | (this is basically Martin's jet module which was renamed CoolingJet in 3.0) |
MTI | ||
MultiJets | ||
StreamDisk | ||
TrueLoveProblems |
New 3.0 Problem Modules
Module | Should it be tested? |
---|---|
PulsedJet | |
MachStems | |
CoolingJet | |
EmissionMaps | no |
Science Meeting Update 02/17/2014 - Eddie
Mach Stems
- all of the 2-D 80 cells/rclump runs completed on BlueStreak
- I did some estimates for the SUs required for 3-D runs: ehansen02122014
1-D Rad Shocks
- Anna made the comparison plot, and it showed that implementing the Raymond ionization table did not make much of a difference. There must be a different reason for why our models are different.
- I thought to check the time resolution, but when I looked at one of my models I found that the time step was less than a nanosecond. If I used the right units, I think typical ionization rates lead to a very rough estimate of about 1023 ionization events per cell per time step near the shock front.
Pulsed Jets
- I couldn't get visit to open the chombos on kraken, but I did get it to open the .bov files. Below are the emission maps for the hydro case (left) and the beta = 1 case (right).
hydro | ![]() |
---|---|
beta = 1 | ![]() |
Other
- thought about memory and run time scaling: ehansen02132014
Memory and Run Time Scaling
It seems intuitive and obvious that the amount of memory used (M) should scale with the total number of cells . For a 3-D problem,
So if the total memory M is proportional to N, then it follows from the above equations that
You might also think that the run time
scales in the same way, but that is not quite correct. Memory usage is a function of how much information you have and run time is a function of how many calculations you have to do on that information. The single most important factor determining run time is how many time steps ( ) are required to update each cell. So it now makes more sense to write an expression like this:
The CFL condition is needed to calculate how the number of time steps scales:
where C is the CFL number and
is the maximum wave speed. C and can be ignored since they are constants, and then we can write
So if we have a run where we are trying to get to some final time
, we have to take a certain number of time steps.
Combining this equation with our previous equation for
, we get the final result:
Below is a very simple example that illustrates how these equations work. For simplicity, the 1 x 1 x 1 simulation has all of its values explicitly set to 1.
Grid | N | dx | dt | |||
---|---|---|---|---|---|---|
1 x 1 x 1 | 1 | 1 | 1 | 1 | 1 | 1 |
2 x 2 x 2 | 8 | ½ | ½ | 1 | 2 | 16 |
4 x 4 x 4 | 64 | ¼ | ¼ | 1 | 4 | 256 |
The scaling equations can also be generalized to problems of any dimension ndim.
Estimated SUs for 3-D Mach Stems
Here are some stats for the 2-D runs. These were all run on bluestreak at a resolution of 80 cells/rclump. I meant to run them on 128 cores, but I accidentally ran on 128 nodes which is equivalent to 2048 cores. This means that I should have better efficiency for my 3-D runs, so these estimates should be upper limits.
gamma | d (rclump) | run time (s) | SUs (thousands of CPU hrs) |
---|---|---|---|
5/3 | 20 | 33916 | 19.3 |
5/3 | 15 | 31027 | 17.7 |
5/3 | 10 | 23635 | 13.4 |
5/3 | 05 | 22282 | 12.7 |
1.40 | 20 | 34698 | 19.7 |
1.40 | 15 | 31184 | 17.7 |
1.40 | 10 | 29120 | 16.6 |
1.40 | 05 | 22086 | 12.6 |
1.20 | 20 | 41328 | 23.5 |
1.20 | 15 | 37236 | 21.2 |
1.20 | 10 | 32170 | 18.3 |
1.20 | 05 | 23586 | 13.4 |
1.01 | 20 | 52246 | 29.7 |
1.01 | 15 | 52459 | 29.8 |
1.01 | 10 | 47344 | 26.9 |
1.01 | 05 | 28754 | 16.4 |
Zcool | 20 | 63.9 | |
Zcool | 15 | 61.2 | |
Zcool | 10 | 54.9 | |
Zcool | 05 | 35.6 | |
TOTAL | 524.55 |
If I used the same filling fractions and extended these numbers to 3-D, then that would be like calculating the estimated run time for simulating rods instead of spheres. I went ahead and did it this way anyways, because the grid in the z-direction for these problems is relatively small compared to the other dimensions. So these estimates should be an upper limit.
To get the estimate for the 3-D runs, we need to compare the expected number of cells updates per root step.
2Dupdates = N0*(1 + 8*ff(0) + 64*ff(0)*ff() + 512*ff(0)*ff(1)*ff(2) + …)
3Dupdates = N0*(1 + 16*ff(0) + 256*ff(0)*ff() + 4096*ff(0)*ff(1)*ff(2) + …)
These can be written more generally as sums:
where n is the level, N_0 is the total number of root level cells and ff(i-1) is the filling fraction of level i-1. Note that ff(-1) = 1. It is also important to keep in mind that 3-D cell updates are approximately twice as expensive as 2-D cell updates, so…
3DSUs/2DSUs = 2*(3Dupdates/2Dupdates)
After using the above equations, this is what I get for the 3-D runs:
gamma | d (rclump) | SUs (millions) |
---|---|---|
5/3 | 20 | 22.3 |
5/3 | 15 | 20.3 |
5/3 | 10 | 15.5 |
5/3 | 05 | 14.8 |
1.40 | 20 | 22.5 |
1.40 | 15 | 20.3 |
1.40 | 10 | 19.0 |
1.40 | 05 | 14.6 |
1.20 | 20 | 26.8 |
1.20 | 15 | 24.1 |
1.20 | 10 | 20.8 |
1.20 | 05 | 15.5 |
1.01 | 20 | 33.7 |
1.01 | 15 | 34.2 |
1.01 | 10 | 30.8 |
1.01 | 05 | 18.7 |
Zcool | 20 | 74.2 |
Zcool | 15 | 71.5 |
Zcool | 10 | 63.9 |
Zcool | 05 | 41.6 |
TOTAL | 605 |
Clearly this is not feasible. Even if I could account for better efficiency in 3-D and perhaps lower filling fractions, I don't think I will be able to do all of these runs at this resolution. If I decrease the resolution to 40 cells/rclump then my total drops to 54.78 million SUs. I can improve this further by decreasing my domain in the x-direction and imposing periodic boundary conditions. I have to alter the resolution slightly to do this, so this would be for 38.4 cells/rclump. For this set up, my total drops to 35.56 million SUs. To summarize:
Resolution (cells/rclump) | Periodic BCs? | Total SUs (millions) |
---|---|---|
80 | no | 605.0 |
76.8 | yes | 390.9 |
40 | no | 54.78 |
38.4 | yes | 35.56 |
It is important to remember that these are upper limits, and there are several factors that will bring these numbers down:
- The 3-D filling fractions are probably a bit lower than what I used.
- The efficiency should be better for the 3-D runs. The 2-D runs typically had efficiencies of 15%-60%. In a perfect scenario, my efficiency would be 100% which is a significant increase.
- These estimates are for BlueStreak. Kraken, for example, has processors that are approximately 1.625 times faster (2.6 GHz vs. 1.6 GHz).
If I do some very rough estimates to try to account for these things, I can get the total SUs down to 9.85 million. Perhaps it would be possible and reasonable to do a subset of the runs at this resolution on Kraken. Below, I computed the average SUs and run time required for a run with 2048 cores on BlueStreak and Kraken.
Machine | SUs (thousands) | Run Time (days) |
---|---|---|
BlueStreak | 800.2 | 16.3 |
Kraken | 492.4 | 10.02 |
With run times this long, we run into the issue of having to do restarts due to the limited wall time on these machines. This will bring the efficiency down and increase the SUs and run time required. There is also the extra time spent waiting in the queue which probably doubles the time it takes to get a run completed.
Science Meeting Update 02/10/2014 - Eddie
- ran a 1D shock model with a Raymond H ionization table, sent the data to Pat, and am waiting on comparison results
- the 3D MHD pulsed jet run on kraken slowed down quite a bit, so I took a different approach
- I'm trying a very low resolution run, and then I will do a restart at a higher resolution to get some high res late frames
Mach Stems
- below I am showing the d = 5 runs from left to right: Z cooling, and gamma = 5/3, 1.4, 1.2, 1.01 respectively
- on bluestreak I have completed 6 out of 20 of the runs at the 80 cells/rclump resolution. The remaining 14 are still in the queue, but should run and finish this week.
Meeting Update 02/03/2014 - Eddie
Mach Stems
- finished the low res (20 cells/rclump) runs for stationary bow shocks with different gammas and separation distances
- gamma = 1.01, 1.20, 1.40, 5/3 with separations of d = 5, 10, 15, 20 for a total of 16 runs
- currently running new set of sims with Zcooling (set gamma back to 5/3)
- will visualize and post images/movies of all the completed runs this week
- also this week, I will move to bluestreak and run the same simulations with 2 more levels of refinement to give me production runs at 80 cells/rclump
- finally, I will also start testing low res, 2-D runs with a moving bow shock
3-D Pulsed Jets
- hydro simulation on Kraken finished, and beta = 1 is chugging along
- beta = 1 has a few restarts, but it is still producing more frames than previously
- it is currently at frame 40, and I think getting it to frame 50 will be good enough
- I hope to create the emission maps for these runs as soon as possible so Pat can use them for the HST proposal
Cooling Issues
- looked at the Raymond code and found some discrepancies in the ionization and recombination rates
- not sure if this discrepancy is large enough to account for what we've seen in our 1-D radiative shock models
Other
- wrote an abstract for upcoming HEDLA conference
- need to proofread my 2.5-D pulsed jets paper again
- moved office around so now we can make use of Adam's old futon
Raymond vs. AstroBEAR - Hydrogen Ionization and Recombination
Ionization (H —> HII)
I went through the Raymond shock code and figured out how it calculates the collisional ionization rate for H. I then plotted this (CION, black) as a function of temperature alongside the table that we use in AstroBEAR (AB Rate, red). In AstroBEAR, the units we assume on these values are cm3/s.
The two curves match pretty well which is great. Unfortunately, this also means that the discrepancy between the codes is elsewhere and still unknown. If we have the units wrong in AstroBEAR, then we are doing the scaling wrong which occurs after this ionization rate calculation. Also, the Raymond shock code has photoionization, but AstroBEAR does not.
I also made a plot of the ratio of the two rates. Even though they are on the same order of magnitude, you can see that the Raymond rate is often 10% - 60% higher than the AstroBEAR rate within the temperature range that is relevant to my simulations. I'm not sure if this difference is significant or not.
Recombination (HII —> H)
I did the same thing for HII recombination. This time, I put the y-axis on a log scale to make the difference more visible.
The two rates only differ at very high temperatures. My simulations do not use such high temperatures so this is fine.
Meeting Update 01/27/2014 - Eddie
- Worked mostly on paper edits. I should be finished with them tomorrow if not later today.
- Completed some different mach stem runs: ehansen01232014. I am going to put together a complete set of runs for different gammas and separation distances. I will have to rerun some of my old set ups, but other runs are already completed. I will also start exploring simulations in which I give one of the clumps a velocity.
- Another thing on my to do list this week is to dive into the Raymond shock code. More specifically, I will be looking for subroutines dealing with ionization and recombination and trying to compare these to what we have in astrobear.
- Lastly, I am going to see how far I can get with my 3D pulsed jet simulations. I want to make emission maps for a hydro run and for a beta = 1 run. Then, these could perhaps be used for Pat's HST proposal.
some new 2D mach stem runs
I have started looking at some different cases for my 2D mach stems with 2 clumps. These all have the lower boundary extended to avoid the boundary effects we saw before.
So far, I have run the
case with clump separations of 20, 15, 10, and 5.I will do runs with gamma = 1.4 as well to have a total of 16 runs that fill a parameter space for different gammas and separations.
I also tried a run with the clumps offset. This is with gamma = 5/3, a separation of 15, and an offset of 4.
I won't do any more runs with offsets like this. Rather, I am going to start studying simulations in which one of the clumps has a velocity.
Meeting Update 01/21/2014 - Eddie
- 3D pulsed jets seem to be working better with the density and velocity smoothing ehansen01162014
- the cell position is now in the output with the 'restart due to nan in flux' error
- ran the 2 clump, mach stem simulations for an intermediate gamma = 1.20. I did clump separations of 5, 10, 15, and 20 rclump. Images and movies will be posted soon.
pulsed jet working better with smoothed profile
I implemented a tanh function to smooth the density and velocity radial profiles at the injection boundary. The code runs much better now since there have not been any nans in the flux calculations. Also, there have only been a couple of high CFL restarts which is not a big deal.
I started with smoothing the profiles out to the edge of the domain. Clearly, this makes a very wide jet that is not what we want. But it was a good test to make sure what I had written was correct, and to see how the code handled it.
I then ran with smoothing out to twice the jet radius. This looks like a more realistic jet, but the head of the jet looks weird. You can see what looks like secondary bow shocks that protrude beyond the head of the jet. These are located at the original jet radius (aka where the smoothing begins).
I am waiting on a new run with smoothing again out to twice the jet radius. I fixed an error with the species densities. If it looks the same, I may try adjusting the pressure a little bit to see if the protrusions go away.
UPDATE
Fixing the species densities did not change much. Not sure if changing the pressure is the right way to go, but I am trying it anyways. I did another run with the smoothing out to 1.1 Rjet. This is the closest to what I actually want, and it ran without any issues.
None of the movies are uploaded yet because I have an issue with visit. For some reason when I try to use the 'save' options in visit, the images get yellow-ed out. This seems to only be an issue when running visit remotely from clover with 3D data. Running locally on clover with 3D data is fine, and running remotely from clover with 2D data is fine. Not sure what the fix is; could be a driver issue.
Meeting Update 01/09/2014 - Eddie
Mach Stems
Ran the d = 20 simulation with an extended lower y boundary. The size of the y-domain was doubled from 40 to 80. The distance from the clump center to the bottom boundary was increased from 30 to 70.
Here is a comparison with the previous run. The run with the extended boundary is on the right, but I am not plotting the extended region.
So the bottom boundary does appear to be playing a role in regards to this nonphysical upwards motion. Now we just have to decide what to do about it. The easiest thing would be to extend the boundary (perhaps even farther than I already have), and then ignore the bottom portion of the domain. This would be ok for these 2D runs, but the 3D runs might be different. This also might not work for the runs where I put one or both clumps in motion.
3D Pulsed Jet
Jonathan implemented a new CFL and maxspeed check for the astrobear standard output. When there is a "high CFL restart", it now reports the maxspeed across all processors instead of just processor 0, and it gives the position. Here is an example:
High hydro CFL ( 0.21E+01 > 0.80E+00) found at position ( 0.1875E+00, 0.8929E-02, 0.9911E+00) on processor 3175, due to a maxspeed of 0.27E+03 - Restarting Step
In regards to my pulsed jet problem, this is great for figuring out where in the domain the maxspeed is getting too high. This is most likely due to low densities and thus high alfven velocities.
However, my simulations are also showing restarts "due to nan in flux". We thought these might be caused by the high CFLs and maxspeeds, but these restarts are not preceded by high CFLs. Attached to this post is a copy of the standard out from my most recent run on kraken.
The issue has probably evolved since we first started this problem, but the original ticket created for the 3D pulsed jet runs is ticket #289. The problem could also be related to tickets #301 and #302. For simplicity, I closed #289 and #302, and started a new ticket #330.
Meeting Update 12/16/2013 - Eddie
- made all of the images and movies from the 2 clumps simulations ehansen12052013
- also did some analytic work to show that the equation that I was using for bow shock shape does not do a very good job of matching the bow shocks that I am getting in my simulations ehansen12112013
- checked the output from my 3D pulsed jet run. The first things I checked were pressure and temperature, and the pressure does something strange at the injection boundary.
I am doing some runs with different parameters to see if I can find or fix the issue.
bow shock shape
Now that I did some simulations just looking at the shape and interactions of bow shocks from clumps, I wanted to compare with the formula that I had been using for my mach stem simulations.
It looks like the equation for the shape of the bow shock given by Ostriker (paper) does not work well for different gammas. There is not a big difference in shape between = 5/3 and = 1.01.
From my 2 clumps simulation with separation 40 rclump, you can see that the bow shocks are very different.
Furthermore, the equation does not even come close for my gamma = 5/3 run. I think the Ostriker formula assumes a radiative bow shock, which makes more sense for its shape.
The equation comes closer for gamma = 1.01, but there is still a big discrepancy. Maybe it is accurate for a gamma = 1, but I doubt it comes close enough for the purposes that I need.
This leads me to think that perhaps I have to derive my own formula for bow shock shape from my simulations in order to do a proper mach stem study with astrobear.
Meeting Update 12/09/2013 - Eddie
- Finished all 12 of the 2 clumps runs. Just working on visualizing right now. The rest of the images and movies will be up later today. ehansen12052013
- The 3D pulsed jet run finally made it through the queue on bluestreak. There was something wrong with the system, so they gave me a reservation to make it run. It still produced a lot of restarts though, so I will have to debug that this week.
- Also on my agenda for this week is to make some progress with the non-equilibrium cooling.
new simulations of bow shocks from 2 clumps
Parameters
Now that I have a clear plan for these simulations, I thought I would give a little more detail about the problem set up. These are all 2D, so the clumps are actually slices of cylinders. For simplicity, I am still referring to them as clumps.
For all simulations, these parameters are the same:
Tamb | 1000 K |
---|---|
namb | 1000 cm-3 |
Twind | 1000 K |
nwind | 1000 cm-3 |
vwind | 50 km/s |
nclump | 107 cm-3 |
rclump | 10 AU |
Ly | 40 rclump |
my | 50 cells |
AMR levels | 4 |
tfinal | 329 yrs |
The wind is supersonic with a mach number of approximately 13.5. Based on the wind speed and length of the domain, this run time is about 8.67 wind crossing times. The clump temperature is set to be in pressure equilibrium with the ambient.
Simulations
I have a set of 12 runs that I am working on. They differ by clump separation distance (d) and
. The length of the y domain is fixed, but the length of the x domain changes. The resolution stays fixed at 20 cells/rclump. Here is a list of separation distances and grid size for each type of run:d [rclump] | Lx [rclump] | mx [cells] |
---|---|---|
40 | 48 | 60 |
20 | 28 | 35 |
15 | 24 | 30 |
10 | 20 | 25 |
5 | 12 | 15 |
The run with d = 40 rclump was added to see a set up where the bow shocks did not interact much if at all. Each set up will be run with a
= 5/3 and a = 1.01 for a total of 10 runs. These runs all have open boundary conditions. The remaining 2 runs will be with d = 20 and 15, = 5/3, and they will have a reflecting boundary in the middle instead of simulating both clumps. We want to make sure that the reflecting boundary simulations behave just like the 2 clump simulations.I will post images and movies here over the next few days…
Images and Movies
bow shocks from 2-clumps simulations
I ran four simulations with 2 clumps and a wind. Each run has a different clump separation: 20, 15, 10, and 5 clump radii.
Much of the setup is the same as what I used for the mach stem runs with some slight changes to make the simulations run faster: rclump = 10 AU, nclump = 106 cm-3, namb = nwind = 1000 cm-3, Tamb = Twind = 1000 K, and vwind = 50 km/s. This means the wind is supersonic with a mach number of approximately 13.5.
Once you get down to separations of about 5 rclump, the bow shocks combine into one large bow shock. This is why my mach stem runs have not behaved as expected. The separation distances I have been testing are typically on the order of 2 rclump.
I have a couple of ideas on what to do next:
- The formula that I am using for the bow shock shape could be incorrect or not applicable in this context. I could develop my own formula from simulations.
- It might just be impossible to form mach stems by forming bow shocks right next to each other in this way. A possible way around this is to form the bow shocks away from each other and then bring them closer together.
Meeting Update 12/03/2013 - Eddie
- Marvin is starting to use NEQ cooling in his disk module. I will help with any issues that may come up.
- I need to remake the synthetic emission maps for Andrea. The previous ones show the "bubble", but what he really wants to be able to see is the jet. If there is still emission in the jet then this should just be a matter of scaling differently in visit.
Mach Stems
I have been trying different 2D runs for the mach stems. The first attempt looked the same as the slice from the 3D runs. I thought that the behavior we were seeing might be caused by the wind, so I moved the wind boundary farther away. Here is the result,
The situation here is not as clear cut as a wind of uniform density sweeping over a clump and making a bow shock. I think the fact that the wind creates a forward and reverse shock is changing the behavior from what we might expect.
I then tried setting the ambient and the wind to be the same density. I also ran the simulation for twice as long.
One mighty bow shock is formed, and the right reflective boundary (the one we are supposedly not interested in) seems to be playing a role as well.
I want to turn off the reflective boundaries and just look at the shape of the bow shock that is formed. I could also just keep running the simulation to see if I ever get to a steady state, but I still don't really like all the lateral movement going on. Any other ideas?
Meeting Update 11/25/2013 - Eddie
- Made some cool looking emission maps for Andrea ehansen11212013
- Fixed a memory leak in Marvin's code, but it did not resolve the error. Couldn't fit any other leaks, so maybe it is time for Baowei and/or Jonathan to pick this up?
- I have a test run of my 3D pulsed jet waiting in the queue on bluestreak. I want to see if it can run at least 12 frames without too many restarts before I ask for another reservation.
- I found an error in my mach stem calculations because I was using the wrong angle. I fixed that and changed the simulation set ups that I plan on running. ehansen11192013
- Extended the lower boundary on my mach stem simulation. The result was the same. I will try a 2D run next.
cool emission maps of lab jets
Just wanted to share these images with everyone, because I think they look cool.
These are emission maps of H-alpha (green) and [S II] (red). They are generated from data that Andrea provided from laboratory jet experiments. They are consistent with the concept that H-alpha marks shock fronts, and [S II] follows behind shock fronts in cooling regions.
edit to mach stem runs
I looked at Kris' paper on mach stems and hysteresis: http://adsabs.harvard.edu/abs/2013HEDP....9..251Y
Here they labeled the critical angle as the angle that the bow shock makes with a vertical line, not the included angle between the bow shocks. So my previous calculations to find the critical separation distance were actually using the wrong angle.
With this correction, the critical separations decreased. Unfortunately, some of them decreased to the point where it would be difficult to simulate. In other words, when the critical separation gets really close to 2 rclump it becomes difficult to put the clumps close together below critical because the clumps would overlap. So I had to change my runs to be 5% above and 5% below critical instead of 10% like I had before. Below is the corrected table with the grid domain and resolution for each of my runs:
Run | Mach Stem? | Domain Size (x, y, z) | Base Grid (mx, my, mz) | |
---|---|---|---|---|
A | no | 1.243, 2.486, 0.6215 | 6, 12, 3 | ←— lowest physical resolution at 77 cells/rclump |
B | yes | 1.125, 2.250, 0.5625 | 6, 12, 3 | |
C | no | 1.205, 2.410, 0.6025 | 6, 12, 3 | |
D | yes | 1.090, 2.180, 0.545 | 6, 12, 3 | |
E | no | 1.168, 2.336, 0.584 | 6, 12, 3 | |
F | yes | 1.057, 2.114, 0.5285 | 6, 12, 3 | |
G | no | 1.143, 2.286, 0.5715 | 6, 12, 3 | |
H | yes | 1.034, 2.068, 0.517 | 6, 12, 3 | ←— highest physical resolution at 92 cells/rclump |
According to these numbers, the run that I showed yesterday should not have formed a mach stem. It was at x = 1.242 which is pretty close to what I am now calling run A. Perhaps what we saw was something non-physical due to some boundary effects. I am rerunning that set up with the lower y boundary extended.
Meeting Update 11/18/2013 - Eddie
Pulsed Jets
I fixed a couple of lines in my pulsed jet problem module, and the code ran better without restarting. I am doing one more test, and if everything looks good, I will be ready to try this again with another reservation.
Mach Stems
I came up with a plan for the mach stem runs: ehansen11122013.
I changed the grid and mesh to make these runs go a bit faster…
- The mach stem should occur in the upper z plane, so I cut off part of the clump in the z-direction.
- I implemented some derefinement inside the clump and below the clump.
With these changes, I was able to complete run B in less than 1.5 days on 8 cores on Grass. When I use more cores on bluestreak, I should be able to complete all the runs in a day. I could also increase the resolution and still get the runs done in a reasonable amount of time.
Below are some density and velocity images and movies for run B:
Mach Stem Runs
Following work done by Kris, (see Mach Stems), I redid some of my old calculations.
The formulas Kris and I used in our calculations were for bow shocks from jets, but they should be applicable to bow shocks formed by a wind and clump interaction.
The theory gives a critical angle based entirely on
. Mach stems should form at angles equal to and greater than the critical angle. This does not depend on how the bow shock is formed…critical angle (deg) | |
---|---|
5/3 | 37.4 |
1.4 | 41.7 |
1.2 | 47.7 |
1.1 | 53.5 |
The set up I am doing right now is for two stationary clumps. Their symmetric bow shocks will interact exactly in between the clumps and either form or not form a mach stem. The angle formed between the bow shocks decreases with increasing separation distance. So if the critical angle is a minimum for mach stem formation, then critical separation distance is a maximum.
Based on the formula for the shape of a bow shock, you can get a critical separation distance for each critical angle. Note that the separation distance is measured from clump center to clump center. These are given in units of clump radii…
critical angle (deg) | critical separation (rclump) | |
---|---|---|
5/3 | 37.4 | 2.76 |
1.4 | 41.7 | 2.64 |
1.2 | 47.7 | 2.52 |
1.1 | 53.5 | 2.45 |
Since I am simulating one of the bow shocks with a reflecting wall, my domain has to be half the size of this critical separation. I think it will suffice to do two runs for each gamma: one above critical and one below. Here are my proposed runs:
Run | critical separation (rclump) | simulation separation (rclump) | Mach Stem? | |
---|---|---|---|---|
A | 5/3 | 2.76 | 3.036 | no |
B | 5/3 | 2.76 | 2.484 | yes |
C | 1.4 | 2.64 | 2.904 | no |
D | 1.4 | 2.64 | 2.376 | yes |
E | 1.2 | 2.52 | 2.772 | no |
F | 1.2 | 2.52 | 2.268 | yes |
G | 1.1 | 2.45 | 2.695 | no |
H | 1.1 | 2.45 | 2.205 | yes |
If I use 4 levels of AMR and aim for at least 64 cells per clump radii, I can use these resolutions:
Run | Mach Stem? | Domain Size (x, y, z) | Base Grid (mx, my, mz) | |
---|---|---|---|---|
A | no | 1.518, 3.036, 0.759 | 8, 16, 4 | |
B | yes | 1.242, 2.484, 0.621 | 6, 12, 3 | |
C | no | 1.452, 2.904, 0.726 | 6, 12, 3 | ←— lowest physical resolution at 66 cells/rclump |
D | yes | 1.188, 2.376, 0.594 | 6, 12, 3 | |
E | no | 1.386, 2.772, 0.693 | 6, 12, 3 | |
F | yes | 1.134, 2.268, 0.567 | 6, 12, 3 | |
G | no | 1.3475, 2.695, 0.67375 | 6, 12, 3 | |
H | yes | 1.1025, 2.205, 0.55125 | 6, 12, 3 | ←— highest physical resolution at 87 cells/rclump |
Meeting Update 11/11/2013 - Eddie
I did not get much work done in the past 4 days or so because I was busy packing, moving, unpacking, etc.
Pulsed Jets
The restart that I ran on bluestreak did not give me any more frames. I need to check the output to see if anything looks odd.
Mach Stems
Redid the previous run (
), but this time with density gradient refinement, and no position based refinement. Only got 80 frames out of 100 before our systems got rebooted last week.I still want to go through some of my old notes, redo some calculations, and come up with a grid of problems to run (different gammas and different separations).
Emission Maps
I am finished with the coding aspect of this side project for Andrea. Wrote a module that takes a table of x, y, z, n, T and initializes a grid. The gas is assumed to be purely H, and it uses the equilibrium ionization fraction at the given temperatures. It outputs the initial chombo and also the corresponding emission projections as a bov file and then stops.
All that is left to do is run it and give the data and some images to Andrea.
Other
An astronomy senior at RIT (Dave Anderson) is looking at UR for grad school. He is interested in astrophysics and is curious about out group. If there are no objections, I would like to invite him to next week's group meeting?
I am giving a colloquium at Geneseo on the 21st, so I need to start planning that. I suppose it should just be a general talk about what our research group does. Adam, you have done these kinds of talks before, any suggestions?
Meeting Update 11/04/2013 - Eddie
- 3D pulsed jet with beta = 1 is running on my reservation on bluestreak. We will see how far it runs by Wednesday morning.
- I am running mach stem simulations with density gradient refinement. This will take longer than previous runs, but it may be fine.
- Working on creating synthetic emission maps for Andrea. I need to figure out how to use the emission routines with ASCII data outside of astrobear. Another issue is that Andrea only has data for x, y, z, n, and T. The emission routines require ne, x, and T. Perhaps I can use some assumptions and approximations to at least produce something that is qualitatively relevant.
Meeting Update 10/28/2013 - Eddie
- Getting a reservation on bluestreak to do the beta = 1 3D pulsed jet.
- Looked back at the mach stem project I had started last year, and reran an updated version of the code…
This is one stationary clump, and I have reflected the data to make the image look like the physical situation. Now I want to simulate two clumps at the same time and give on of them a velocity.
Meeting Update 10/21/2013 - Eddie
Pulsed Jets
Beginning of last week, I made a correction to the cooling length in my pulsed jet simulations…ehansen10162013a
I also made a 3-slice movie of density for the partial run. This is 12 frames (out of 100) of the
case.I am not sure why visit made the images like this. I know it is very distracting, so I will try to fix this.
Radiative Shocks / Cooling
We still do not fully understand the discrepancy between astrobear and the Raymond shock code. It is clear that cooling is implemented differently.
After talking with Pat about cooling in astrobear vs. cooling in the Raymond code, I tried figuring out the differences analytically…ehansen10172013
Pat did not comment on my calculations, but I got the impression that this was not the correct way to compare the two different methods?
Mach Stems
Looked back at what I had started last year and some papers to refresh my memory. Will start doing some runs this week.
Cooling: internal energy vs. enthalpy
Cooling in AstroBEAR is currently handled by removing energy from the internal energy. In the Raymond code, cooling is different, but it appears to be equivalent to removing energy from the enthalpy. I have derived the resulting pressure for both cases below.
Internal Energy Formulation (current implementation)
The internal energy is defined as
.
Let
be the energy loss due to cooling. The new internal energy is,
.
Therefore, the new pressure
is.
For an ideal gas with
.
Enthalpy Formulation (proposed change)
The enthalpy is defined as
.
After cooling, we have
,
,
.
Therefore, the new pressure is
.
For an ideal gas with
.
It appears that the enthalpy formulation would result in removing less thermal pressure, and it would therefore have a higher temperature in comparison to the internal energy formulation.
correction on cooling length resolution in pulsed jet sims
So I actually made a mistake in my cooling length calculations… I was modeling a 100 km/s shock in my 1D code to get a feel for what the cooling length in my pulsed jet sims should be. I thought this was correct since the injection velocity varies from 150 - 250 km/s. However, the maximum shock velocity is actually half of that, so 50 km/s.
By comparison, a 50 km/s shock will have a smaller cooling length, so this is actually good news. It means that my cooling length resolution is better than I thought.
Below is an updated table of the resolutions used from other papers.
Paper | Base Grid | AMR | Physical Size [AU] | Jet Radius [AU] | Lcool [AU] | Cells/AU | Cells/Rj | Cells/Lcool |
---|---|---|---|---|---|---|---|---|
Mine (2013) | 42 x 420 | 7 | 2005 x 20054 | 334 | 7.29 | 2.68 | 896 | 19.53 |
de Colle (2006) | 180 x 1800 | 0 | 2005 x 20054 | 334 | 7.29 | 0.09 | 30 | 0.66 |
Raga (2007) | 16 x 64 | 10 | 1671 x 6685 | 134 | 10 | 9.80 | 1314 | 98.05 |
Tesileanu (2012) | 128 x 384 | 6 | 400 x 1200 | 20 | 0.46 | 40.96 | 819.2 | 18.81 |
This also has implications for the 3D runs that we want to do. The extremely high cpu hour estimates are for higher cooling length resolutions which we may not need. Below is an updated table of the final results I had previously posted.
Base Grid | Levels | Cells/Lcool | SUs (millions) |
---|---|---|---|
84 x 420 x 84 | 7 | 19.53 | 645.0 |
84 x 420 x 84 | 6 | 9.77 | 116.6 |
84 x 420 x 84 | 5 | 4.88 | 19.79 |
21 x 105 x 21 | 7 | 4.88 | 19.79 |
21 x 105 x 21 | 6 | 2.44 | 3.028 |
21 x 105 x 21 | 5 | 1.22 | 0.424 |
I tried a run on bluestreak with the parameters of the lowest resolution listed in the table above. I ran the
model and got 12 frames in about 18.55 hours on 4096 cores. That is approximately 76,000 SUs.Meeting Update 10/07/2013 - Eddie
Not much to report this week. Just working on the paper: added a lot of text and references, and making figures (ehansen09292013)
Resolutions Used in Various Jet Papers
Paper | Base Grid | AMR | Physical Size [AU] | Jet Radius [AU] | Lcool [AU] | Cells/AU | Cells/Rj | Cells/Lcool |
---|---|---|---|---|---|---|---|---|
Mine (2013) | 42 x 420 | 7 | 2005 x 20054 | 334 | 3 | 2.68 | 896 | 8.04 |
De Colle (2006) | 180 x 1800 | 0 | 2005 x 20054 | 334 | 3 | 0.09 | 30 | 0.27 |
Raga (2007) | 16 x 64 | 10 | 1671 x 6685 | 134 | 10 | 9.80 | 1314 | 98.05 |
Tesileanu (2012) | 128 x 384 | 6 | 400 x 1200 | 20 | 0.2 | 40.96 | 819.2 | 8.19 |
For Raga (2007), the numbers in the table make sense for a 10 level run. These values are the same as what is actually stated in the paper However, they said that they used 11 levels. So either their simulation was actually 10 levels, and not 11, or someone did the resolution calculations wrong, or there is some detail about their refinement that I am missing.
For Tesileanu (2012), I did not see what their time-dependent velocity was, so I am not sure how strong the internal shocks are. This would determine the cooling length. I will try looking at this more closely to see if I can figure it out from their data.
UPDATE
I changed the above table a bit. Before I was using an Lcool = 4 AU for my calculations. That is still a pretty good benchmark for most of the internal shocks. However, the strongest shock would be at 100 km/s. That is the peak jet velocity of 250 km/s running into the minimum velocity material which is at 150 km/s. I ran these parameters through my 1D code and got Lcool = 3 AU. This is the minimum cooling length that I would expect to see for the internal shocks.
This calculation was done for the hydro case. However, this should be a minimum for all cases (ie MHD with different betas), because magnetic fields would only increase Lcool.
I did the same calculation for the parameters given in Tesileanu (2012), and got Lcool = 0.2 AU. This is mainly due to the increased density. This is for their models which have njet = 104. They also have models at njet = 5 x 104, so these would have an even smaller minimum cooling length.
I also did the calculation using the Raga (2007) parameters. I got Lcool = 1.9 AU. This is because they also use a higher density than in my sims. However, the cooling length they reported and used was 10 AU. They did the calculation for a 40 km/s shock, while I did it for an 80 km/s shock…For a 40 km/s shock, astrobear gets 8.30 AU which is pretty close to their 10 AU.
So my question is this: When you have the sinusoidally dependent injection velocity, is it more appropriate to consider your strongest shocks to be equal to the amplitude of your sinusoidal function, or 2x the amplitude?
Characterizing the Forces on the Clumps in Pulsed Jets II - Isothermal Approx.
As discusses in the meeting today, it would be more accurate to do my calculations in the isothermal limit. So I set gamma = 1.0, and got some interesting results.
Turns out that in this limit, the force is always negative (compression regime). As a reminder, the forces are:
The magnetic field and pressure profiles are as follows:
Where
is the beta that I define in my initial conditions. is derived from and . is some intermediate radius, set to 0.6 rjet in my set up. depends on and . These are the profiles for . Outside of the profiles are such that the forces go to zero.Solving for the net force as is results in zero, but if I apply shock jump conditions first, then I can write
where I have defined
as the pressure jump and as the density (or magnetic field) jump. It is important to remember that the jump conditions are radially dependent due to the fact that the pressure is radially dependent. In the adiabatic case, this gets pretty complicated. If I enforce the conditions for an isothermal shock, then it simplifies a bit. To maintain the same temperature across the shock, the pressure jump has to match the density jump, so .
So making the problem isothermal guarantees that the force is always negative; the clump is always in the compression regime. This is consistent with what De Colle and Raga found in their 2008 paper: the non-radiative cases were dominated by thermal pressure forces, and the radiative cases were dominated by magnetic forces.
Even though I made the problem simpler, there are still some interesting things to learn from the plots…
It is easy to make sense of these plots. The force is greater for greater velocities and mach numbers because the magnetic field gets compressed more. More magnetic field means stronger pinching. The other cases of beta follow the same trend, so I will not post them.
If you only looked at the mach number plot, you would naively think that compression is stronger at larger radii. However, the plots follow this trend because mach number increases with radius. At any given moment, the velocity is constant across the shock, so you might think that stronger compression occurs at smaller radii. You would again be wrong. I made a different plot to explore this point.
This plot shows that for a given beta, you will get stronger compression with smaller radii but only to a certain point. Then, the curve turns over and the force starts decreasing. I also noticed in this plot that the different cases for beta were not in order. Strange right? So I made yet another type of plot to check this out.
Indeed, you do not keep getting stronger compression with stronger fields. Could it be that at some point you are "saturated" with magnetic field and the force really stops increasing, or is this a flaw of the isothermal approximation? In reality, the shock always increases the temperature. I played around with gamma in my spreadsheet, and I found that as I increased it to say 1.10, the above plot began to turn the other way.
So yes, for increasing gamma you get weaker compression, but you get compression that always increases with decreasing beta.
Meeting Update 09/30/2013 - Eddie
!D Rad Shocks
Redid the original 18 shock models for both x = 0.01 and x = 0.99. The low ionization models did not change that much, but improved a bit. The high ionization models improved a lot, but the temperature comparison is still a little off.
Here is an example of one of the models…
2.5D Pulsed Jets
Working on getting more references together and reading more. Fiddled around with latex this morning to get a bibliography file working, so now I can easily add all my references to the paper as I gather more. I am currently going through the PPVI review paper, and some papers by Tesileanu.
Trying to get the paper figures going, but visit has been really slow on our local machines. I thought it was the AMR in my data, but I turned the top 3 levels off and it is still slow. This might just take some more time and patience. The density image and movie look pretty cool though: ehansen09292013.
Did some analytic work on the magnetic and pressure forces on the clumps inside the jet beam: ehansen09262013. Still trying to figure out if this numerical analysis can be made more analytic. The problem is that there does not seem to be any place where I can make an approximation, so the equations are pretty ugly. Furthermore, knowing the force is one thing, but how do I use that to derive something that can be calculated from my simulations?
3D Pulsed Jets
I still have a run sitting in the queue on Bluestreak. It is a fairly large job, so that is why the wait has been so long. There is not much urgency here, but hopefully it will run this week.
The scaling tests from Stampede looked good. Thanks Baowei!
Characterizing the Forces on the Clumps in Pulsed Jets
In this problem, I am considering a magnetic pinching force caused by the toroidal field and the gas pressure force. We have:
The jet starts out in magnetohydrodynamic equilibrium so that the magnetic force and gas pressure force cancel out (F = 0). When the jet material gets shocked at an internal working surface, the magnetic field increases and the gas pressure increases. They increase such that the previous equilibrium is no longer preserved. The clump will expand if F > 0, and if F < 0 the clump will compress.
I used the MHD shock jump conditions to find the resulting F inside a clump. The jump conditions actually depend on radius in this problem. For this reason, it is extremely difficult to write an analytic expression for F. I decided to study F numerically.
Based on the above equations, you would expect to get clump compression as the magnetic field is increased. The first time I ran the numbers, I could not get F < 0. I knew this had to be wrong, so I started playing around with my initial parameters. I found that the sign of F was highly sensitive to the strength of the shock (typically characterized by the mach number M).
In general, it seems that low M shocks are easier to compress. Once the mach number goes below 1.0, the material is no longer shocked, and equilibrium can be reestablished (F = 0).
Below are a few plots for the three cases of beta that I'm running in my simulations (5, 1, 0.4).
![]() | ![]() |
![]() | ![]() |
![]() | ![]() |
The way I like to think of these plots is that the relative velocity plots will tell you about clump evolution. The mach number plots can tell you more general things about field strengths and shock strengths. Several things that we can tell from these plots:
- The clumps are only in the compression regime at certain relative velocities. Typically, smaller relative velocities (aka weaker shocks). This means that the clumps will be more susceptible to compression at later times.
- The inner regions of the clump are more susceptible to compression first. This means that as the relative velocity decreases, the inner regions of the clump will start compressing first. We could call this "inside-out" compression.
- At some point during the compression phase, the outer regions begin to be compressed with greater force than the inner regions. This seems counter to point 2. If the entire clump enters the compression regime quickly, then the "inside-out" compression might look more like "outside-in" compression, or possibly uniform compression.
- As the magnetic field strength increases (beta decreases), the compression region shifts to the right on the relative velocity plots. However, this does not mean that stronger fields are capable of compressing stronger shocks. The mach numbers would also be shifted to the right on this relative velocity plot. Despite this fact, the mach number plots show a similar trend. So yes, stronger fields can compress stronger shocks.
- At some point, the relative velocity becomes so small that a shock is no longer formed. This is the same as the mach number going to 1.0. This means that the pressure and magnetic field do not follow shock jump conditions, and thus the net force returns to zero.
- The peak magnitude of the force increases as beta decreases. This confirms that stronger fields compress with greater force.
Meeting Update 09/23/2013 - Eddie
- made some better estimates for how long 3D runs will take ehansen09172013
- working on paper figures
- need to redo some 1D radiative shock models for Anna and Pat
SUs required for 3D pulsed jet runs
I took a look at the standard out from the 2.5D simulations that were completed on Stampede. Here are estimates for the total number of cell updates required at the end of the simulation and the number of SUs for the entire run:
Beta | Cell Updates | SUs |
---|---|---|
Hydro | 181139252 | 3145 |
5 | 144062238 | 5270 |
1 | 170121573 | 7638 |
0.4 | 164839549 | 11472 |
The total number of cell updates was easily obtained from the filling fractions in the standard out. The following equations were used:
N0 = nx*ny N1 = 4*N0*ff(0) N2 = 16*N1*ff(1) = 16*N0*ff(0)*ff(1) ... 2.5Dupdates = N0 + 2*N1 + 4*N2 + 8*N3 + ... = N0*(1 + 2*ff(0) + 4*ff(0)*ff() + 8*ff(0)*ff(1)*ff(2) + ...)
Now, going from 2.5D to 3D leads to many more computational cells. There is a factor of 2 from the fact that we can no longer impose symmetry and only simulate half of the jet. So to keep the same resolution the base grid goes from 42 x 420 x 1 to 84 x 420 x 84. That is 168x more cells on the base grid. Figuring out how many cells there are at higher levels is a bit trickier.
I analyzed the 2.5D data in Visit to get an estimate for what the filling fractions would be in 3D. I used the expression
r = coord(Mesh)[0]
Then, I made a pseudocolor plot, and queried the weighted variable sum for each AMR level separately. The ratio of these sums gave me the filling fractions for each level. Then, I used these equations to get the number of cell updates:
N0 = nx*ny*nz N1 = 8*N0*ff(0) N2 = 64*N1*ff(1) = 64*N0*ff(0)*ff(1) ... 3Dupdates = N0 + 2*N1 + 4*N2 + 8*N3 + ... = N0*(1 + 2*ff(0) + 4*ff(0)*ff() + 8*ff(0)*ff(1)*ff(2) + ...)
If I assume that the number of root steps remains constant then the ratio of (3Dupdates) / (2.5Dupdates) will tell me how much longer the 3D runs will take. However, there is an additional factor of 2 because Kraken processors are about half as fast as Stampede processors. Also, there is yet another factor of approximately 2 because 3D cell updates are more expensive than 2D cell updates. So to convert the SUs columns from table 1 to table 2 you would use:
3DSUs = 4*(3Dupdates)/(2.5Dupdates)*(2.5DSUs)
Below are my estimates for the number of cell updates, and how long the 3D runs should take on Kraken (with 10 cells per cooling length).
Beta | Cell Updates | SUs (millions) |
---|---|---|
Hydro | 1178305614616 | 81.8 |
5 | 956767405342 | 140 |
1 | 1132203286417 | 203 |
0.4 | 1361329586503 | 379 |
In total, this is about 804.1 million SUs on Kraken. Clearly, this is not going to work, so we have to lower the resolution. I redid the calculation for different resolutions. Below is a table of the total SUs for all 4 runs.
Base Grid | Levels | Cells/Lcool | SUs (millions) |
---|---|---|---|
84 x 420 x 84 | 7 | 10 | 804.1 |
84 x 420 x 84 | 6 | 5 | 140.4 |
84 x 420 x 84 | 5 | 2.5 | 24.5 |
21 x 105 x 21 | 7 | 2.5 | 24.5 |
21 x 105 x 21 | 6 | 1.25 | 4.03 |
21 x 105 x 21 | 5 | 0.625 | 0.576 |
So the only one that seems feasible is the last one. Now, one thing to keep in mind is that this Lcool is based on a cooling length of 4 AU which is what we think the strongest cooling regions within the jet should be. This corresponds to the high density and high velocity shocks (1000 cm-3 and 60 km/s). There are shocks that are stronger than this, especially at the head of the jet and near the injection region when pulses enter the grid. We may never be able to resolve these fully.
So even though a resolution might show 0.625 cells/Lcool, we are still fully resolving cooling lengths of 64 AU with 10 cells. This is adequate for most of the low density (100 cm-3) shocks, and even some of the low velocity shocks regardless of density (< 40 km/s). So this should be fine for the clumps that have already traveled some distance. At the last frame of the simulation, this looks like it'd be good enough for at least 3 of the 5 clumps, maybe 4.
In other words, at this resolution, we do not resolve the cooling length at the beginning of a clump's "lifetime", but we do end up resolving it later.
UPDATE
As was discussed in today's meeting (9/23), I redid all the calculations for frame 65 rather than frame 100. The ratio of (3Dupdates/2.5Dupdates) is different at different frames due to filling fractions being different. Furthermore, during a simulation, the filling fractions in 2.5D won't necessarily change in the same way that they do in 3D. The idea here is that the filling fractions around 2/3 of the simulation time are most characteristic of the time-averaged filling fractions.
I only redid the last table, since it is the most useful. The numbers improved, but it appears that these runs still require many SUs.
Base Grid | Levels | Cells/Lcool | SUs (millions) |
---|---|---|---|
84 x 420 x 84 | 7 | 10 | 645.0 |
84 x 420 x 84 | 6 | 5 | 116.6 |
84 x 420 x 84 | 5 | 2.5 | 19.79 |
21 x 105 x 21 | 7 | 2.5 | 19.79 |
21 x 105 x 21 | 6 | 1.25 | 3.028 |
21 x 105 x 21 | 5 | 0.625 | 0.424 |
Meeting Update 09/16/2013 - Eddie
Pulsed Jets
2.5D
Runs are done and the paper is coming along nicely.
- Baowei completed the 7-level runs on Stampede. I looked at one of them last week, and it looked fine. Next step is to analyze this data and make the figures for my paper.
- The main text of the paper is more or less finished. I just need to put it in latex, and then I will send Adam a pdf later today. I will be working on the figures this week. It will also need some more references, and of course some edits.
3D
Conducting test runs on Kraken; using the rest of the old allocation which will expire at the end of the month.
- I tried a hydro run first. The resolution is 4x less than the 2.5D runs, which means in the strongest cooling regions I will only have about 2.5 cells per cooling length. The simulation was only able to advance 5/100 frames in 12,096 SUs (12 hrs on 1008 cores). This means that if filling fractions remained constant, it would take about 240,000 SUs to complete the run. However, filling fractions increase as the jet propagates, so realistically this run will take closer to 1 million SUs. We may have no choice but to back off on the resolution, at least for most of the run.
- There won't be that much time between the 2.5D and the 3D papers. Should these remain isolated publications, or should I give them titles that go together like a "series"?
1D Rad Shocks
Met with Pat last week, and we discussed the next steps that need to be taken. * He is working on figuring out some scaling factors in the Raymond code so that he can compare cooling rates to the ones that I gave him from AstroBEAR.
- I will be rerunning a set of 1D models and give the data to Anna for comparison.
- We discussed the paper that will come of this. Pat sees it as more of an AstroBEAR verification paper and a discussion of the cooling processes, rather than something that gives an in depth quantitative analysis of the magnetic field and emission lines. Hence, the bulk of the writing will be done by Pat and I.
Meeting Update 09/09/2013 - Eddie
The issues I have been having with the pulsed jets simulations seem to have been resolved. I had set up a bunch of test runs with different resolutions to try and diagnose the problem, but none of the runs triggered restarts (PulsedJets). There must be something in the data files or maybe a slight modification that I made to my code which fixed the problem.
The output of my Kraken runs looks fine (I'll be posting results on the wiki soon). Baowei has been running MHD versions on Stampede (bliu09092013). These also look good, so Baowei is going to let these runs finish.
Here is my plan for the rest of this week (in order of priority):
- Finish writing the 2.5D pulsed jets paper
- Make figures for the paper
- Use rest of the time on Martin's kraken allocation to test 3D pulsed jets
Also this week is the LLE meeting on Thursday. Hopefully, I will also get a chance to talk with Pat about our 1D radiative shock models.
Meeting Update 09/03/2013 - Eddie
Pulsed Jets
Spent most of my time last week trying to get the 2.5D pulsed jet simulations working. Met with both Jonathan and Baowei a few times. We had decided on a set-up to try on kraken that we were sure would work. However, restarts are still being triggered, and it is still unclear why.
From previous runs, the refinement looked strange. There were small patches of high refinement scattered throughout the jet beam. After changing a few parameters and adding some error flag buffering, I was able to clear this up. The latest kraken run is at the resolution that I want. It does density gradient based refinement up to level 3, and cooling time threshold refinement to level 7. You can see from the first few frames that the mesh looks reasonable…
Frame | Density and Mesh |
---|---|
0 | ![]() |
1 | ![]() |
2 | ![]() |
3 | ![]() |
4 | ![]() |
Baowei was making better progress on bluehive and on our local machines, so I decided that maybe the problem with the restarts is machine dependent. I have a run going on Grass which seems to be progressing just fine. It would probably take on the order of weeks to finish completely though since it's only running on 8 cores.
I have a run on bluestreak using 256 cores which is progressing faster than the run on Grass, but is also triggering restarts. I am currently setting up a run on bluehive, so we'll see how that goes.
Other Stuff
- Haven't heard back from Pat yet regarding the numbers I had sent him.
- I plan on finishing the paper within the next week or two. I will write it as though i have all the results and figures, which can be added in later. Maybe this will just be a 2.5D paper though?
- Went to a Geneseo physics alumni event last Saturday. We were celebrating the physics department's 50th anniversary, and also the fact that Geneseo now leads the country in number of physics bachelor degrees awarded yearly. There were some other interesting things as well…
- I briefly met their new astronomy professor, Annie Pellerin.
- My work that I did as an undergrad there is referenced in a new textbook by Koberlein and Meisel (a former professor of mine), and it is titled Astrophysics Through Computation: With Mathematica Support.
- They're looking for potential colloquium speakers and asked me if I would be interested. So I might be giving a colloquium there either this semester or in the spring.
- RIT has a student named Dave Anderson who is a senior this year. He's interested in going to graduate school for astronomy, and he's considering UR. More specifically, he is interested in our research group and would like to meet all of us to see if this might be a good fit for him.
Meeting Update 08/19/2013 - Eddie
Pulsed Jets
I ran into a couple of errors with the 2.5D pulsed jet runs last week. After some troubleshooting, I finally figured out that it was memory issues due to the projections (emission maps). Jonathan has some ideas for improving the code, but I think the best thing for me to do right now is to lower the number of levels for the projections.
So my simulations will be refined on cooling time at 42 x 420 + 7 levels, but the projections will be at 42 x 420 + 4 levels (or 672 x 6720 since projections are actually fixed grid images).
I did a quick test to make sure this would work. Previously, my runs were stalling at the first frame dump, so I couldn't get any output. By decreasing the plot level for the emission maps, I got 4 frames in less than 5 minutes on 1080 cores. So it looks like these runs will work now.
1D Radiative Shocks
I exchanged a bunch of emails with Pat today. There are some discrepancies between our results, but they follow the same trends. We are working on sorting out some of the discrepancies. I sent him the fully ionized model that he requested, and I am also going to be rerunning all of the original 18 models and sending him the data tomorrow.
Lab Jets
Francisco got back in touch with me and referred me to a paper on plasma simulations from Spain. These simulations use Kr and Xe plasmas, while Francisco's experiments are of an Al plasma jet entering an Ar medium. There are definitely some parallels here, so I will read this paper more closely to learn more and see if any of it might be useful for astrobear.
resolution needs for my 2.5D jet sims
The run information for my original 2.5D runs was purged from my Kraken directories, so I don't have any run-time estimates from those.
I did a 32 x 320 + 3 run on Grass today. It took 42 minutes on 8 procs (5.6 cpu hrs). This is with the new cooling time refinement criteria.
Kraken processors are similar to Grass processors in terms of speed (2.6 GHz vs 2.5 GHz). The resolution that I think I need is 42 x 420 + 7. This should take 441 times as long, thus requiring about 2374.62 cpu hrs.
Meeting Update 08/12/2013 - Eddie
- I got updates to my pulsed jet code (thank you Jonathan and Baowei), and now I will test these changes on Kraken. I'm also still working on the paper.
- In regards to the 2.5D runs…the most recent resolution that I ran these at was 32 x 320 + 4. This gave an effective resolution of 3.92 AU/cell, and 85.33 cells/jet radius. According to a plot that I made a while ago (https://astrobear.pas.rochester.edu/trac/astrobear/blog/ehansen04232013), the minimum cooling length for typical internal shocks in this set up should be approximately 4 AU. Therefore, if we use these numbers, to ensure 10 cells per cooling length we would need 10 times the aforementioned resolution. This would be 40 x 400 + 7.
- Anna compared some of out 1D radiative shock models and this is a plot that she sent me.
The behavior of our results look similar which is good. One thing to note is that my code for whatever reason always gives slightly longer cooling distances. I think this makes sense if the code that Anna is using has more emission lines and therefore stronger cooling.
Also, she said that our emission line fluxes were off by an order of magnitude, so she's looking into why.
Meeting Update 07/23/2013 - Eddie
- PPVI went really well: met a lot of new people, shared my poster, got my emission map movie in Adam's review talk, had a really good time in Germany!
- Haven't had too much opportunity to review the 3D pulsed jet stuff. I know Baowei and Jonathan have been working on it and updating the ticket, so I will check that out soon. I'll be in email and phone contact while I'm in California for the next two weeks, so I'll try to help out as much as I can.
- I re-ran some 1D radiative shock models, reformatted the data, and sent it to Pat. Hopefully, that is all he needs, but if he needs more I can probably do any additional work for this fairly quickly when I get back.
Meeting Update 07/08/2013 - Eddie
Just got back into the swing of things today, so not much to report in terms of results. Here's where I'm at now with my various projects:
Pulsed Jets Paper
The paper is coming along fine. I can only work on the 2.5D parts, since the 3D runs are still blowing up.
3D Jets (precessing)
Implemented a precessing injection velocity and tried running a few simulations. I have low resolution runs for theta = 3 and 1 deg. The precessing period is 1/3 of the pulsation period. The resolution is 24 x 120 x 24 + 1.
I tried a higher resolution (3 levels) on kraken for theta = 2, 1, and 0.1 deg, but they failed to complete. I backed down to 2 levels and am trying the theta = 1 deg run again.
3D Jets (testing minimum temperature)
At the meeting two weeks ago, Jonathan mentioned that some of my simulations had these "flashes" in the temperature plot. We suspected it might be pressure protections, so I tried a run where I increased the minimum temperature. That run completed and the data is on Kraken waiting to be visualized.
3D Jets (fixing the explosions)
We still don't know for sure what caused the explosions, but we've seen some interesting grid effects at the injection boundary. Today and tomorrow, I am going to implement a sub-sampling routine into my problem module. This should be a much more accurate way of initializing radial profiles on a Cartesian grid.
3D Jets (colliding)
We heard back from Sergei and Andrea, and they are going to give me some initial parameters soon that I can use to run some comparable simulations.
1D Radiative Shocks
Just got an email from Pat today. He needs the data in a specific format to make plots, so I will work on reproducing that stuff for him soon.
Injection Boundary for 3D Pulsed Jets
This is related to ticket #289.
It looks like the point at which things really blow up was at a restart (restarted from frame 33, explosion at frame 34), so maybe there is something about my code that makes restarts a bad idea. Regardless, the simulation was at a hault before I tried doing the restart.
Meeting Update 06/24/2013 - Eddie
- made movies of the 3D pulsed jet runs
beta | image | movie(s) |
---|---|---|
hydro | ![]() | movie |
5 | ![]() | movie |
1 | ![]() | full movie; expansion (slowmo) |
0.4 | ![]() | movie |
- I tried increasing the size of the domain for the beta = 1 run, but the expansion started happening sooner…
- I also altered the refinement at the injection region in two ways: I used the cylindrical shape option instead of a rectangular prism, and I increased the shape size by 50% in all directions. The shape of the mesh didn't change (it's still rectangular), so I don't know if the cylindrical option did anything differently. These runs haven't finished yet, but I can already tell that the beta=0.4 run seemingly worked better. The image below looks like an improvement, but the simulation on kraken has come to a hault just like the old one did.
- Made image/movie of emission maps for the 2.5D runs
- working on results section of paper and adding more references
Meeting Update 06/17/2013 - Eddie
- Worked on paper: practically done typesetting up to the results section for which I am just starting the text now.
- 2.5D runs ran fine on Kraken, so those are done.
- 3D runs: hydro and beta = 5 runs finished, but the beta = 1 and beta = 0.4 runs have not been able to finish. I posted a ticket for this: #289. Kraken is down at the moment, but I will attach the standard out and other files to this ticket once Kraken is back up.
- Finished a colliding jets run at a higher resolution. I had emailed Sergei a few weeks ago to let him know that I started on this project, but I never received a response.
- Back to the 1D radiative shock models: I don't remember hearing back from Pat when I finished the runs and posted all those images: ehansen04222013. Is there anything I have to do for this project right now, or am I just waiting for something to happen on Pat's end?
Meeting Update 06/03/2013 - Eddie
- Dealing with a variety of problems for my pulsed jet runs: ehansen05312013. The errors that I have gotten on bluestreak don't make sense to me. They reference a line in the code that looks perfectly fine.
- Have most of the text written for abstract, introduction, and methods sections of my paper. Just need to typeset it with equations and references.
Issues with Pulsed Jet runs on Different Machines
2.5D Runs
We are doing these runs at a higher resolution now. Approximately 3.5x the resolution used by Raga. I have tried running these on alfalfa (16 procs) and bluestreak (256 procs). Below is a summary of what is going on…
Beta | Alfalfa | BlueStreak |
---|---|---|
Inf (hydro) | completed in ~11 hrs | restarted from frame 34, then from 35, then from 88 and completed (can't tell exactly what total runtime was, but looks to be ~7 hrs) |
5 | quits at frame 6.6 with error message (see attachment) | completed in ~11.6 hrs |
1 | ? | restarted from frame 79, now cannot get past frame 82 with restarts |
0.4 | ? | restarted from frame 34, now cannot get past frame 61 with restarts |
All the runs that I have had to restart on bluestreak quit with the same error. See the attached files/images for the standard output from the beta = 0.4 restart, and also the corresponding error messages from core 64.
3D Runs
I have been trying to do these runs on kraken on 480 procs…
Beta | Kraken |
---|---|
Inf (hydro) | completed in ~1.6 hrs |
5 | restarting from frame 89, but last estimated wall time before restarting was > 2 days |
1 | restarting from frame 51, but last estimated wall time before restarting was > 2 weeks |
0.4 | restarting from frame 40, but last estimated wall time before restarting was > 2 days |
There are no error messages from these runs. They just start taking forever because they keep requesting restarts. I attached the standard output from the restarted beta = 5 run.
UPDATE
Turning off refinement for B-field gradients helped the run on alfalfa get past frame 6.6, so I'm trying this on the other machines as well for the runs that have not yet completed.
Meeting Update 05/20/2013 - Eddie
I was really busy last week, and out of town for a few days, so I didn't get as much done as I would have liked. Here's where I'm at right now…
- My 3D pulsed jet runs are having memory issues on both bluestreak and kraken. I ran valgrind but found no memory leaks, so I'm not sure why this is happening. I could probably get some lower resolution runs through, but I definitely can't resolve the cooling length if this doesn't get fixed. I made a couple of tickets so that Jonathan, Baowei, and I can work on this issue.
- I altered my jets module to handle colliding jets. I also implemented a switch to do co-rotating magnetic fields and counter-rotating fields. Here are some low resolution images and movies which show that my modifications worked:
- I didn't get any paper writing done last week, but I should have enough free time this week to make some progress.
Some cool colliding jets
So it wasn't too difficult to modify my pulsed jet module to handle colliding jets. I did a quick low-resolution, hydro run to test out my modifications, and it looks pretty cool.
Density | Emission Map |
---|---|
![]() | ![]() |
movie | movie |
Next, I have to edit the magnetic field part, but this shouldn't be too difficult either.
Meeting Update 05/13/2013 - Eddie
- 3D pulsed jets are working…ehansen05102013. I have a set of runs waiting to go on blue streak.
Paper Outline
- Abstract
- general overview of paper
- Introduction
- mention observations
- motivation for this research
- what questions does this paper address:
- simulations such as these have been done before, but only for H-alpha emission, and only in 2.5D…what new information do we obtain by also looking at [S II] emission?
- also, what new information do we obtain by extending the simulations to 3D? are "nose-cones" still present?
- how does magnetic field strength affect the clumps and emission within a pulsed jet?
- can we use H-alpha and [S II] emission to estimate magnetic field strength in YSO jets?
- Numberical Simulations
- Methods
- equations being solved with astrobear and type of cooling
- 2.5D
- define the problem set-up and parameters
- 3D
- more brief description of problem set-up
- Methods
- Results
- 2.5D
- density plot, a few lineouts (rho, v, T), emission map
- 3D
- density slice, a few oineouts of the slice (rho, v, T), emission map
- 3D image
- 2.5D
- Discussion
- more qualitative than quantitative
- differences between 2.5D and 3D
- Conclusions
- summarize the paper
- relate simulations back to observations
- lingering questions and future work
- are we likely to learn anything new from simulating other magnetic field geometries?
- also, is there anything else to learn from looking at other emission lines such as [N II] or [O I]?
- References
3D pulsed jets appear to be working!
So it was a long and annoying task, but Jonathan and I finally got it working. Most of the changes I had made to my module to go from 2.5D to 3D were correct. The two biggest changes that were made today which appeared to make it work were:
- Doubling the problem domain. Now simulating half the jet instead of a quarter.
- Switched interpolation order from 3 to 2.
Jonathan suspects there might be a bug in the 3rd order interpolation scheme which caused assymetries in the z-direction in my simulations. I'm going to run a couple more tests and post a ticket on this.
Here is a low resolution movie of the H-alpha and [S II] emission map from a beta = 5 run.
Meeting Update 05/09/2013 - Eddie
- Star Formation Jamboree was fun even though there were no talks, and hardly even any mentioning, of jets
. Nevertheless, I had a good time talking with Fred Adams, Ralph Pudritz, and many other people.
- Made a poster yesterday for the CIRC Poster Session, so come by Goergen tomorrow anytime from 10-12 and VOTE FOR MY POSTER! It is by no means perfect, but it will suffice for now. I will probably use a different, but similar, one for PPVI, especially since I will have 3D runs before then.
- The 3D runs are still producing infinities. I've tried different things to fix it, but nothing has worked yet. I think with Jonathan and/or Martin's help, this should be working by tomorrow if not by the end of today.
- For the 1D runs, we are now just waiting to hear back from Pat. I'm gonna go ahead and do some analysis on the emission lines to see if I can better quantify the effect of the B-field.
Meeting Update 04/30/2013 - Eddie
- Finished the drift analysis for my 1D radiative shock models. You can see that the drift is indeed quite small. ehansen04262013
- Submitted introduction slide for the Star Formation Jamboree.
- Completed hydro run for 3D pulsed jet. The MHD runs quit due to divergence protections. I'm going to try initializing the field in a different way to see if I can fix this issue.
The run here is low resolution, a fixed grid of 45 x 450 x 45 which is 4 times less than the effective resolution of the original 2.5D runs.
Drift of the 1D radiative shocks
I looked at the drift of the shock front for one computational time for all 18 models. One computational time is different for each model since the length scale is different for each model, so in order to compare apples to apples I decided that vdrift / vshock was the appropriate diagnostic.
The drift was not measurable for vs = 60 km/s. Even in the worst case, the drift velocity was only 0.065% of the shock velocity. Therefore, I think this effect is negligible. Obviously, if we ran the simulations for many computational times, the shock would eventually drift off the grid but that kind of simulation will take too long and will never be necessary. Still, there must be a very small error somewhere. I suspect it is somewhere in the chemistry or scaling, but again, the error is so small that we don't need to worry about it.
Meeting Update 04/23/2013 - Eddie
- Did all 18 1D radiative shock models ehansen04222013
- Trying to get 3D pulsed jet runs working. Had login, compile, and submission issues on Kraken and BlueStreak. Those are resolved now, but I am still having some runtime issues. Hopefully, they get resolved today.
1D Radiative Shock Models
Here is a table of the 18 models that I ran with plots of the hydrodynamic variables and the emission lines.
vs (km/s) | log(n) (cm-3) | beta | hydro | emiss |
---|---|---|---|---|
40 | 2 | 5 | ![]() | ![]() |
40 | 2 | 1 | ![]() | ![]() |
40 | 2 | 0.5 | ![]() | ![]() |
40 | 3 | 5 | ![]() | ![]() |
40 | 3 | 1 | ![]() | ![]() |
40 | 3 | 0.5 | ![]() | ![]() |
50 | 2 | 5 | ![]() | ![]() |
50 | 2 | 1 | ![]() | ![]() |
50 | 2 | 0.5 | ![]() | ![]() |
50 | 3 | 5 | ![]() | ![]() |
50 | 3 | 1 | ![]() | ![]() |
50 | 3 | 0.5 | ![]() | ![]() |
60 | 2 | 5 | ![]() | ![]() |
60 | 2 | 1 | ![]() | ![]() |
60 | 2 | 0.5 | ![]() | ![]() |
60 | 3 | 5 | ![]() | ![]() |
60 | 3 | 1 | ![]() | ![]() |
60 | 3 | 0.5 | ![]() | ![]() |
vs (km/s) | log(n) (cm-3) | beta | hydro | emiss |
---|---|---|---|---|
40 | 2 | 5 | hydro_v40_n2_b5.png | emiss_v40_n2_b5.png |
40 | 2 | 1 | hydro_v40_n2_b1.png | emiss_v40_n2_b1.png |
40 | 2 | 0.5 | hydro_v40_n2_b0.5.png | emiss_v40_n2_b0.5.png |
40 | 3 | 5 | hydro_v40_n3_b5.png | emiss_v40_n3_b5.png |
40 | 3 | 1 | hydro_v40_n3_b1.png | emiss_v40_n3_b1.png |
40 | 3 | 0.5 | hydro_v40_n3_b0.5.png | emiss_v40_n3_b0.5.png |
50 | 2 | 5 | hydro_v50_n2_b5.png | emiss_v50_n2_b5.png |
50 | 2 | 1 | hydro_v50_n2_b1.png | emiss_v50_n2_b1.png |
50 | 2 | 0.5 | hydro_v50_n2_b0.5.png | emiss_v50_n2_b0.5.png |
50 | 3 | 5 | hydro_v50_n3_b5.png | emiss_v50_n3_b5.png |
50 | 3 | 1 | hydro_v50_n3_b1.png | emiss_v50_n3_b1.png |
50 | 3 | 0.5 | hydro_v50_n3_b0.5.png | emiss_v50_n3_b0.5.png |
60 | 2 | 5 | hydro_v60_n2_b5.png | emiss_v60_n2_b5.png |
60 | 2 | 1 | hydro_v60_n2_b1.png | emiss_v60_n2_b1.png |
60 | 2 | 0.5 | hydro_v60_n2_b0.5.png | emiss_v60_n2_b0.5.png |
60 | 3 | 5 | hydro_v60_n3_b5.png | emiss_v60_n3_b5.png |
60 | 3 | 1 | hydro_v60_n3_b1.png | emiss_v60_n3_b1.png |
60 | 3 | 0.5 | hydro_v60_n3_b0.5.png | emiss_v60_n3_b0.5.png |
Meeting Update 04/16/2013 - Eddie
- My pulsed jet module is working in 3D. Now, I need to move to a larger machine and do some runs. I will try both BlueStreak and Kraken.
- I am working on generating the data and plots for the 18 models of 1D radiative shocks. There will be a blog post before the end of the week containing my version of the plots. If everything looks good, the raw data will be sent to Pat so that he can create his own plots for comparison.
- There is still a small amount of drift in the 1D shock simulations. This is not a problem for generating the data and plots since those will come from the initialization frame. However, it would be nice to quantify this drift a bit more, so I plan on making some sort of plot to characterize this. Plotting vdrift/vshock vs. vshock for a few different shock velocities makes the most sense with multiple lines for different pre-shock densities, or betas.
Meeting Update 04/08/2013 - Eddie
- Did cooling length analysis on one of my 2.5D jet simulations: ehansen04042013
- Made movie of 2.5D lab jet. Shown is the total density, the aluminum (jet), the argon (ambient, and the mesh.
- This week, I want to focus on getting some 3D runs done.
Pulsed Jet Simulations and Cooling Length
I'm setting up a diagnostic variable, so that I can see what the cooling length is everywhere on the grid in my pulsed jet simulations.
I have… (cooling length) = (velocity)*(internal energy) / (cooling strength)
Where cooling strength is in units of energy/time.
This should give the value to me in computational units. My simulation has 30 cells per cu at the highest level. Say I want at least 10 cells per cooling length, so if my cooling length drops below 0.33 then I'm under-resolved.
The only thing I"m not sure about is if it's correct to use the velocity, because it will be 0 in the ambient. So the cooling length would be zero everywhere, but maybe that's ok and we just ignore the ambient?
Below is an early snapshot of the hydro run with the Raga set up. The hydro case should have the strongest cooling since magnetic fields tend to puff up cooling regions. You can see that the cooling length starts at 0 because it's in the ambient. Using this criteria the region directly behind the bow shock, and behind that first internal shock appear to be under-resolved.
Meeting Update 04/01/2013 - Eddie
- Got the new 1D radiative shock module running now, but the shock is not steady. It moves to the left, which is in the direction of the preshock gas. It's not a very fast movement, so it looks like something is just a little off. I'm trying different things to see if I can get it to stay put, but the initial profile looks good. I think it's safe to run some 1D models for Pat now.
- For the 1D models, I will basically do a small subset of the original data set that Pat had sent earlier this year. The idea is to avoid doing 100+ runs and just give Pat a handful of runs to compare astrobear to the raymond shock code. So my grid will be vs = 40, 50, 60 km/s, log(n) = 2, 3, beta = 5, 1, 0.5, x = 0.01. That's a total of 18 runs, but if Pat wants to look at preionization as well, then I will add 18 more runs with x = 0.1.
- I started implementing 3D capability into my pulsed jet module. I'll be running and testing that out this week. I also need to put some more thought into the specific runs that I want to do. I see the main focus of my paper being on the effect of magnetic fields on emission in a pulsed jet. Emission maps from a 3D simulation like this have not been presented in the literature to the best of my knowledge.
- I did a 2.5D run of the lab jet for Francisco. I implemented tracers for the aluminum jet and the argon ambient, and the code uses a separate cooling table for each. The tables come from an old paper, but this is what Francisco was using to estimate the cooling strength of the experiments. I will generate some images and movies, and send Francisco an update. I also think it might be worthwhile to see if I can find better tables that span the temperature range of the simulation.
Meeting Update 03/18/2013 - Eddie
- made movie of emission maps for all cases of beta for the 2.5D pulsed jets (de Colle & Raga set up)
- merged Jonathan's new cooling routines with my revision. Sorted out a bunch of merge conflicts and got it to compile.
- updated my radiative shock module to track all hydrogen and helium species. Not getting a steady shock, but I think I've isolated the issue and I should have that fixed by tomorrow.
- ran a 2.5D jet with lab parameters from Francisco. jet radius = 1 mm, jet velocity = 100 km/s, runtime = 500 ns. The jet is an aluminum plasma, and the ambient is an argon gas. The cooling functions for these elements come from an old paper by Post et al., 1977. This simulation is simply using a cooling rate that adds the Al and Ar functions. I think I can improve this by adding tracers and using a density weighted function.
Meeting Update 03/04/2013 - Eddie
Progress Last Week
- finished first round of 2.5D pulsed jet runs with the Raga set up ehansen02272013
- set up 2.5D projections to create emission maps
- started the strong cooling jet project again (Francisco & Imperial College)
- implemented new cooling table and set up the problem with the experimental parameters
Stuff for This Week
- finish coding NEQ cooling
- adding / adjusting a few more tables for some of the He processes
- run some 1D shock models with the finished cooling routines to compare with Pat
- possibly rerun the Raga set-up for 2.5D pulsed jets, or move on to new simulations (maybe random pulses or Kris' set up)
- get the strongly cooling jet with lab parameters running and analyzed
Pulsed Jet Emission
I now have 4 different 2.5D pulsed jet runs, each at 2 different resolutions. The high resolution is at the same resolution that de Colle & Raga used. The low resolution is half of that. The low resolution runs seem less reliable in terms of cooling. The higher resolution runs show more structure, and the jet head evolution is different. This must be due to resolving the cooling length. In the image below (hydro runs), the low resolution run is on the left, and the high resolution run is on the right:
Since this resolution dependency also affects the emission, I have chosen to only do emission maps from the high resolution runs.
Emission Maps
Using the new 2.5D projections object, I was able to create emission maps. In theory, I could make maps for any emission line that I want, but I am only going to show Halpha and SII (6716 + 6731) for now. There are a couple of different ways to visualize the emission data, but I find the truecolor plots to be pretty neat. The scaling for the truecolor plots is done in a way that Jonathan described in a previous post (johannjc02202013).
Halpha is red, SII is green, and where they mix is yellow.
Beta | Inf (hydro) | 5 | 1 | 0.4 |
---|---|---|---|---|
Emission | ![]() | ![]() | ![]() | ![]() |
Clump Emission Analysis
Figure 4 of de Colle & Raga shows a more quantitative analysis of how the magnetic field affects the Halpha emission for each of the 5 clumps within the jet. They plotted the ratio (MHD / Hydro) of the total Halpha emission for each clump. The result is pretty conclusive; as B increases, Halpha increases. I have tried constructing a similar plot with the low resolution runs, but the results don't follow a conclusive trend, so again I am only going to show the plot for the high resolution runs.
The total Halpha emission of a clump is found by summing the Halpha emission inside a circle of radius 1.4 centered at each clump. The clump positions are taken to be z = 6.2, 12.5, 19.37, 25.7, 32.53.
The numbers for the high resolution runs didn't make much sense either. Some clumps showed increase in Halpha with increase in B, while others showed decrease or little change in Halpha. I need to rethink how the Halpha and other emission lines are being calculated, or maybe come up with a better diagnostic.
2.5D Projections
Jonathan helped me alter the projections object to handle 2.5D data. The data is revolved around the symmetry axis (cylindrical z-axis = cartesian y-axis), and then summed along the line of sight (cartesian z-axis). This can be done for any field defined in fields.f90, but it is especially useful for creating emission maps.
How it Works in Astrobear
The projections object was also altered to handle an array of fields, instead of just one. Now you can have several projections in one .bov file, whereas before we would have a separate .bov file for each projection. Here's how you would initialize it:
CALL CreateProjection(Projection) Projection%dim = CYLINDRICAL_PROJECTION Projection%field(1)%id=Halpha_Field Projection%field(1)%name="Halpha" Projection%field(2)%id=SII_6716_Field Projection%field(2)%name="SII_6716" ...
The CYLINDRICAL_PROJECTION parameter tells astrobear that this is 2.5D data, and it needs to revolve the data before integrating. Projections are always output to .bov files, but this is handled a little differently depending on how many fields your projection object has. If there is only one field, then the .bov file will contain the field name, and the variable in visit will appear as the field name. If there is more than one field, then the .bov file will contain the name 'projections', and the variable in visit will appear as a vector or array named 'projections'.
Visit
Visit can be a little quirky when handling a multi-component variable from a .bov file. I have seen that in visit 2.2.2, you cannot see each component separately. Instead, you have to write expressions to access each component individually. If the projection object has 2 or 3 fields, visit will read the variable projections as a vector. However, visit does not know how to use a 3-component vector within a 2D data set. So let's say for example that the 3 fields within your projection object are density, pressure, and temperature. In visit 2.2.2, you would have to write these expressions:
density = projections[0] OR density = array_decompose(projections,0) pressure = projections[1] OR pressure = array_decompose(projections,1) temperature = array_decompose(projections,2) ! notice no choice for the 3rd component
If your projections object contains more than 3 fields, then visit will read the projections variable as an array, and then you have to use the array_decompose function for all elements of the projections array.
Jonathan found that visit 2.5 is a little better. This version of visit gives you the components of the projections vector or array individually. It names them consecutively as comp00, comp01, etc. There is still no way of naming the components by their appropriate field names before visit reads the .bov files. However, at least in this version of visit (and I would assume newer versions as well), you can rename the components with expressions without having to worry about that vector and array_decompose nonsense.
Meeting Update 02/25/2013 - Eddie
The past couple of weeks have been completely devoted to fixing non-equilibrium cooling, 1D radiative shocks, and 2.5D pulsed jets.
Jonathan has done an awesome job of breaking down the DM cooling curve into its components. Now we can more accurately use equilibrium values to reproduce the DM curve, or we can use non-equilibrium values to get a more dynamic cooling curve due to ionization/recombination effects. In comparing with models done by Pat and his student Anna, it has become evident that our cooling routines probably have different physics. Or at the very least we are not putting the same importance or weight on various cooling processes. A closer look at the code that Pat is using will tell us more in the near future.
I was able to reproduce the pulsed jet simulations done by de Colle & Raga. The similarities in the plots look very good ehansen02212013. The emission maps need a little more work.
Pulsed Jets - de Colle & Raga set up
I'm rerunning my pulsed jet simulations, because the previous runs did not appear to be affected very much by the magnetic field. I"m using the exact set up used by de Colle and Raga, 2006. In the interest of time, my runs are at ½ the resolution of de Colle and Raga.
Quantity | Value |
---|---|
Ambient Density | 100 cm-3 |
Jet Density | 500 cm-3 |
Ambient Temperature | 10,000 K |
Average Jet Velocity | 200 km/s |
Jet Velocity Amplitude | 50 km/s |
Jet Velocity Period | 50 years |
Total Run-time | 500 years |
Jet Radius | 5 x 1015 cm |
Domain Size (r, z) | (3 x 1016 , 3 x 1017) cm |
Resolution | 45 x 450 + 1 |
Beta | inf (hydro) | 1 (weak field) | 0.4 (strong field) |
---|---|---|---|
log(density) | ![]() | ![]() | ![]() |
emission | ![]() | ![]() | ![]() |
For the lineouts…solid is hydro, dashed is beta = 1 (weak field), dotted is beta = 0.4 (strong field). The head of the jet, which is to the left, is only included in the 1/beta plots.
Quantity [cu] | Image |
---|---|
density | ![]() |
ionization fraction | ![]() |
velocity (z-dir) | ![]() |
temperature | ![]() |
1/beta (r = 0) | ![]() |
1/beta (r = 0.6 Rjet) | ![]() |
Pulsed Jet Images
I finished the first round of pulsed jet simulations and generated some images. The pulsations are generated by a sinusoidal function with the average jet velocity being 60 km/s. The ambient density is at 100 cm-3 and the jet is injected with a density of 1000 cm-3. The three types of runs that I am showing differ in magnetic field strength. For the pseudocolor plots they are ordered hydro, beta=1, and beta=0.4 from left to right. The lineouts are straight down the axis of the jet, and for those we have the solid lines (hydro), dashed (beta=1), and dotted (beta=0.4). I should also clarify that the head of the jet is to the left for the lineouts.
Quantity | Image |
---|---|
log(density) | ![]() |
log(density) | ![]() |
Ionization Fraction | ![]() |
Velocity (z-direction) [km/s] | ![]() |
Temperature [K] | ![]() |
Emission (multiplied by cell area) | ![]() |
Emission (scaled to maximum) | ![]() |
What is this nonequilibrium cooling nonsense?
Jonathan and I have met a few times now to discuss nonequilibrium cooling since he is implementing the molecular hydrogen routines and variable gamma. While trying to understand how NEQCooling modifies DMCooling, we realized that it was not being modified correctly. First, let me define some of the main parts in the different types of cooling.
Hex_eq = equilibrium hydrogen excitation Hex_neq = nonequilibrium hydrogen excitation DM_em = metal excitation via electron/metal collisions (from DM) DM_Hm = metal excitation via hydrogen/metal collisions (from DM) Z_em = metal excitation via electron/metal cololisions (from Pat)
Now, here is how the different types of cooling are related to each other:
DM = Hex_eq + DM_em + DM_Hm NEQ = DM - Hex_eq + Hex_neq Z = NEQ - DM_em + Z_em
NEQ
For the NEQ case, I think there was some confusion in the past and this was not handled correctly. Hydrogen excitation was being double counted in a weird way. There was a routine for an equilibrium ionization fraction table, which should have been used to get Hex_eq, but it was not used and I could not find a table. To fix this, I implemented a function to calculate the equilibrium ionization fraction:
! solves for equilibrium ionization fraction of hydrogen (derived from Saha equation) FUNCTION EqIonFrac(Temp, nH) REAL(KIND=qPREC) :: EqIonFrac, Temp, nH, b b = (2d0*Pi*emass*Boltzmann*Temp/planck**2d0)**(1.5d0)*EXP(-13.6d0*ev/(Boltzmann*Temp))/nH IF (b>10d10) THEN EqIonFrac = 1d0 ELSE EqIonFrac = (-b + SQRT(b**2d0 + 4d0*b))/2d0 END IF END FUNCTION
The equation for EqIonFrac is quadratic, and b is used with the quadratic formula to get EqIonFrac. The IF statement is in there because I learned a valuable lesson about limits and large numbers. You can see as T becomes large, the exponential in b goes to 1 and b scales as T1.5. So b also becomes large as T becomes large. Due to rounding, the computer takes the square root part to be approximately b as b becomes large, and thus EqIonFrac to be zero. However, if you take the limit of EqIonFrac as b goes to infinity you find that it is 1 (which makes more sense physically). I found that 1010 is a good limit to use for b in that IF statement.
Z
The implementation of Z is slightly more complicated. Pat's cooling tables have a relatively limited range. Outside of this range, Z should revert back to NEQ. This transition needs to be smooth, so DM_em is not completely subtracted from NEQ and Z_em is not completely added to NEQ. These two components are weighted depending on the temperature.
Another issue is that in the code, we do not have a good way of separating DM_em and DM_Hm from the DM curve. So for Z cooling, Z_em is actually replacing DM_em + DM_Hm. This is okay as long as DM_Hm is small in the range of Pat's table. According to DM, DM_Hm is only important when the ionization fraction is low (approximately < 0.001). In Pat's table, the ionization fraction never goes below 0.01, so DM_Hm is not important and can safely be ignored.
Below is an example of how the fixes affected the temperature profile behind a particular shock. DM is black, NEQ is blue, Z is red…the old version is dashed, and the new version is solid. vs = 30 km/s, namb = 100 cm-3, B = 0 (hydro run)
Meeting Update 02/11/2013 - Eddie
Pulsed Jets
I altered my 2.5D jet module a bit and added forced refinement of the injection region (it's at 16 cells per jet radius). I completed 4 runs: beta = inf (hydro), 5, 1, 0.4…need to visualize the data now.
Rad Shocks
Spent almost all weekend practically rewriting my 1D rad shock module because I found that at longer run times the shocks were not steady. I rewrote the jump conditions in a different way, and was much more explicit about tracking the electron number density. I also started accounting for electron mass. For a long time I couldn't get it to work, but I just figured out this morning that it was due to the velocity scale. My density was being scaled much more accurately now, but in physics_control.f90 VelScale ~ SQRT(1/mu). In all other scale factors (except rScale), mu cancels out so I actually don't need to multiply the temperature by mu like I originally thought in cooling.f90. I still need to run some tests, but the shock looks rock steady now.
Updates on Radiative Shocks and Pulsed Jets
1D Radiative Shocks
I have been trying to come up with a better way to figure out the cooling length and thus the appropriate length scale for these runs. While I was reading through the literature, I found a better way to calculate the immediate post-shock conditions, and I made some neat plots.
The points represent some of the actual runs that I am doing based on the grid of parameters that Pat sent me. Each line represents a different shock velocity. You can see that stronger shocks get closer to that maximum compression ratio of 4. I should note that this is the initial compression ratio across the shock jump; with cooling, the density compresses much more behind the shock. Also, a decrease in beta (stronger field) results in a smaller compression ratio. Furthermore, the field has no real noticeable affect until around beta = 10, and then a much greater affect as it gets close to and less than 1.
Similarly, the immediate post-shock temperature decreases with increasing field strength.
Figuring out the cooling length is much less straightforward. This is because it depends on the cooling function lambda. Even if I had a simple function for lambda (which I don't), the cooling length that you might calculate at the shock front is not what the actual cooling length will be. This is because lambda, as well as all other variables, are changing as you move farther behind the shock. The only solution I see is to figure out the cooling length numerically on a case by case basis with trial and error. However, that's essentially what my simulation does and it defeats the purpose of trying to predetermine the cooling length.
2.5D Pulsed Jets
I found a very useful paper that covers a lot of what I am trying to do (de Colle & Raga 2006). There are several main differences that my research will explore:
- I will be using Pat's cooling tables which, as a reminder, include forbidden line emission.
- The differences between helical and toroidal field geometries.
- I have implemented the capability for pseudo-random jet velocities instead of sinusoidal pulses.
- Results will show other emission lines besides H-alpha.
- Eventually, moving to full 3D simulations might change the evolution.
I have a few runs of these jets with different beta. These runs are not pulsed yet, they have a constant jet velocity of 60 km/s, and they are characterized by their minimum beta.
minimum beta | density | beta | emission |
---|---|---|---|
INF (hydro run) | |||
346.6 | |||
38.2 | |||
3.11 | |||
0.026 |
A note on pressure profiles and beta
A lot of times, you'll see these jets defined by a plasma beta as:
Meeting Update 02/04/2013 - Eddie
- Rerunning 1D radiative shocks with a fixed length scale and adjusted plot ranges to send to Pat. Also, working on an analytic form for cooling length dependent on magnetic field strength. I'll use that to add the length scale as another parameter that changes for each run.
- For the 2.5D pulsed jets, I fixed the magnetic field and pressure profile following the form used by Lind in his 1989 paper. So now I'm doing some runs with varying maximum magnetic field strength. Below are images and movies of density, H-alpha, SII, and beta. This particular jet is not pulsed (v = 60 km/s). I changed some of the other parameters as well: namb = 100, njet = 1000, Tamb = 1000 K. Tjet is no longer set to a constant, it now follows the aforementioned pressure profile.
2.5D Pulsed Jet Runs
My initial thought to inject a jet was to have a velocity profile that decreases with increasing radius, and have pressure and magnetic field be constant inside the jet. However, I decided to try it a different way following Lind 1998. They keep the jet velocity constant, and give the pressure and magnetic field a radial profile. Here is a first try at the run with very strong magnetic field (707.11 uG). I also changed the size of the grid to see more of the evolution of the jet.
Meeting Update 01/28/2013 - Eddie
1D Rad Shocks
Done with formatting plots for 1D radiative shock runs. Some of the models don't look quite right, so I have to investigate a little further as to why. Hopefully, it's just a plotting issue where I will have to adjust the yrange in gnuplot and perhaps the length scale of the simulation. Here are some examples (v=30 km/s, n=100 cm-3):
B (uG) | Hydro Variables | Emission Lines |
---|---|---|
1 | ![]() | ![]() |
30 | ![]() | ![]() |
2.5D Pulsed Jets
Coding and initial runs for the 2.5D pulsed jets are finished. Here are the parameters and a few movies of the H-alpha emission. Since AMR is on, the H-alpha emission calculated for each cell is multiplied by the cell area. The magnetic field is helical.
Density (cm-3) | Velocity (km/s) | Temperature (K) | |
---|---|---|---|
Ambient | 500 | 0 | 8000 |
Jet | 5000 | 60 (initial), 20-100 (random pulses) | 10000 |
B (uG) | Image (frame 80) | movie |
---|---|---|
7.07 | ![]() | Halpha_b7.gif |
21.21 | ![]() | Halpha_b21.gif |
70.71 | ![]() | Halpha_b70.gif |
212.13 | ![]() | Halpha_b212.gif |
707.11 | ![]() | Halpha_b707.gif |
Documentation
Had a documentation meeting last Friday (ehansen01252013) to discuss new group organization and updates that still need to be done for the wiki. Will make a more detailed table of things that need to be done and who should be doing them.
Documentation Meeting - 01/25/2013
Last Time
- Jonathan's post from last time: johannjc01092013
User Guide
- Erica's new user's guide: UserGuide
Teams and Leaders
- new group organization: Teams
Blog Categories
- cleaned up some of the redundant/unnecessary blog categories. Using these will be nice for team meetings. Is there a way to have a drop down menu of the most commonly used categories (say top 5) when creating a new post?
Discussion Plugin
- discussion plugin is not working at the moment, but I had also cleaned that up and there is a newer forum (created by Erica) for discussing wiki documentation.
To-Do List
Task | Member |
---|---|
Update objects pages | Eddie |
Update specifics within User Guide | Eddie, Erica |
Update multi-physics pages | Shule, Eddie |
Debug Team Page | Baowei |
Education Team Page | Erica |
Development Team Page | Jonathan |
Organize Developer Guide | Jonathan |
Meeting Update 01/21/2013 - Eddie
- I finished the script that automatically runs astrobear through the entire parameter space for my 1D radiative shocks AND creates the plots. It is kind of awesome that I can do 87 runs and 2 plots per run in about an hour.
- The plots of the higher shock velocities do not look as good. It might be that it takes a couple of more cells to resolve the shock, but I want to investigate a bit to make sure the shock is steady.
- I will finish coding and start running the 2.5D pulsed jet today and tomorrow. Expect images later this week.
- I did a little bit of work on documentation for the wiki. Mostly just moving/renaming a few pages and creating teams and leaders. Should have a documentation meeting later this week.
Teams and Leaders
This table also has its own page at Teams. There are now team pages for each team where the leaders can organize useful long-term information. Information for meetings will probably still be handled by the blog, but it's up to the leader to decide how he/she wants to use their team page.
Team | Leader | Description |
---|---|---|
Documentation | Eddie | review updated wiki pages, discuss outdated sections and organizational issues |
Debug | Baowei/Martin | tickets: review closed, check progress of active, assign new. tests: check current, discuss need for new |
Development | Jonathan/Shule | check status of enhancements in progress, discuss new ideas and what changed in recent code revisions |
Education | Erica | code mini-lessons, seminal paper discussions, community outreach |
Meeting Update 01/14/13 - Eddie
1D MHD Radiative Shocks
I made some modifications to my problem module and generated some plots with my gnuplot script. Every run that I do will have these two plots (or something very close to these). The first is of the hydro variables: log(number density), log(temperature), and ionization fraction. The second is emission lines. Pat gave me some useful feedback, so I just have to tweak these a bit.
I already made the hydro plot look much better. I instead plotted computational number density and temperature, and I made the y-axis a log scale.
Wiki Documentation
We had a useful meeting last week about wiki organization and documentation. Erica has already made significant progress on the User Guide. I have only done a few minor things, but I will be in charge of updating some introductory things within the User Guide, as well as object pages and multiphysics pages.
We briefly discussed the possibility of having leaders and/or teams. We might not necessarily need teams if everyone is part of every team. However, I think it would be useful to have a leader to run a meeting maybe every other week or once a month. This format would also work for debugging, code development, public relations, etc.
Meeting Update 01/07/13 - Eddie
progress since last meeting
- completed about half of the runs for 1D rad shocks with B field
- started writing module for a 2D/2.5D pulsed jet
to do list
- write a gnuplot script to create plots for the 1D rad shocks. This will make plotting for all the runs easy, and also allow Anna to generate comparable plots with her data.
- use new Al cooling table to run lab simulation of jet (Francisco)
on the back burner
- implement cylindrical coordinates
Meeting Update 12/18/2012 - Eddie
- did a few preliminary runs of 1D radiative shocks with magnetic fields. All three runs are with shock velocity 50 km/s and preshock density 100 cm-3. In uG, black is 1, blue is 30, and red is 100.
- working on a script to automatically do all the runs that Pat wants
- constructed an Al cooling table to use for lab jet sims for Francisco, and compared it to the plot in Post et al., 1977…
Post Al Cooling Rate | My Constructed Cooling Table |
---|---|
![]() | ![]() |
Meeting Update 12/03/2012 - Eddie
- made progress on understanding cylindrical source terms and algorithm in the Skinner & Ostriker paper on Athena in cylindrical coordinates. ehansen11272012
- my token is activated now…how are we splitting up our resources?
- currently rerunning some mach stem simulations for longer time to see if they reach a steady state with mach reflection
This Week
- finish rerunning mach stems
- add aluminum cooling table that Francisco gave me and run his jet set up…will probably need kraken for this
- start running 1D radiative shocks with magnetic fields using parameters that Pat sent
different versions of the cylindrical source terms
I believe I have discovered the difference between the source terms I was deriving and those presented in Skinner & Ostriker. It is emphasized that the source terms presented in that paper are used only for the reconstruction step. Those source terms are derived from simplified versions of the MHD equations which make use of div B = 0 (see example at bottom of this post).
If you leave the MHD equations in their full flux-conserving form, then you get the source terms that I previously posted. Below, I have summarized all the possible source terms from these two methods. "full" corresponds to using the complete flux version of the MHD equations, and "divergence-free" corresponds to what you get when you first expand the MHD equations and take out div B.
Conservative Variables
From this point on, these source terms are written to conveniently compare with cylindrical.f90 in astrobear. The "actual" source terms are these values multiplied by a factor of -1/r. Also, the magnetic terms usually come with a factor of 1/4pi, but all these values are in computational units.
variable | full source term | divergence-free source term |
---|---|---|
0 | ||
0 | ||
Primitive Variables
Momentum —→ velocity which is not as simple as just dividing the momentum source terms by the density. Those have to take into account a product rule that introduces another source term. The density and magnetic field source terms remain the same. Energy —→ thermal pressure which also has to be derived on its own.
variable | full source term | divergence-free source term |
---|---|---|
0 | ||
0 | ||
0 | ||
Now the only part that is inconsistent with Skinner & Ostriker is the divergence-free version of the Br source term.
MHD Contributions in the Momentum/Velocity Equation
In the momentum equation, there is a term that looks like:
If you brute force this equation and take the gradient of this tensor, you get the magnetic contributions for the momentum/velocity source terms as shown in the "full" columns above. However, you can also use a tensor identity to rewrite this expression as follows: Then, you could get rid of the div B term since div B should = 0. However, there are source terms associated with div B, so your overall source terms have now changed. Using this method results in the "divergence-free" columns above.Similarly, the MHD contributions for energy/pressure and B also change if you use other calculus identities to expand the magnetic terms and get rid of any div B terms.
cylindrical source terms
I have been double checking the cylindrical source terms in astrobear, because the Athena paper (link) that I was referencing had some differences. I derived the source terms myself…this time using the full 3D MHD equations, and my results were different from Skinner and Ostriker. Astrobear has two cylindrical source term routines for both the conservative and primitive formulations. I get the same results as the conservative routine. However, the primitive routine has some errors, but those are easy to fix. Below is a table showing what I believe the source terms should be…
Conservative
From this point on, these source terms are written to conveniently compare with cylindrical.f90 in astrobear. The "actual" source terms are these values multiplied by a factor of -1/r. Also, the magnetic terms usually come with a factor of 1/4pi, but all these values are in computational units.
variable | source term |
---|---|
0 | |
0 | |
Primitive
Momentum —→ velocity which is not as simple as just dividing the momentum source terms by the density. Those have to be derived separately (see Jonathan's comment). The density and magnetic field source terms remain the same. Energy —→ thermal pressure which also has to be derived on its own.
variable | source term |
---|---|
0 | |
0 | |
So in conclusion, I guess it's possible that I am not understanding some of Skinner and Ostriker's assumptions and/or techniques, but I think they might have some errors in their source terms.
Meeting Update 11/19/2012 - Eddie
- Exchanged some emails with Francisco, and now I have some jet parameters to try out. It seemed a little odd to me that they used this aluminum cooling table (Post, et al 1977) to estimate their cooling parameter, but in simulations using Gorgon they used something else. I will give this Al cooling table a try, and see what we get.
- Read the mach stem thesis. Completed half of the initial mach stem runs, and created a new page for documenting everything on this project: MachStems. Hopefully, Kris is ready to start looking at these.
- Met with Jonathan about implementing cylindrical coordinates. I expect to make significant progress on this within the next month or so.
- All that token nonsense is done on my end, just waiting for activation now.
- Have not heard from Pat or Anna in a while about the radiative shocks.
Meeting Update 11/12/2012 - Eddie
- got my NICS token in the mail, still need to activate it
- read up on cylindrical coordinates in Athena, http://adsabs.harvard.edu/abs/2010ApJS..188..290S. Now I need to find and alter the applicable routines in astrobear.
- I am running the mach stem simulations now, at a higher resolution. Images/movies will be posted soon. I did all the analytic calculations, so I have a set of runs planned.
I was thinking I'd do for each gamma a run at Rcrit, and a run ± 0.2 Rcrit…so with 4 different gammas, that's a total of 12 runs.
gamma | Rcrit/Rclump |
---|---|
1.6667 | 1.4679 |
1.4 | 1.3932 |
1.2 | 1.3229 |
1.1 | 1.2766 |
Relating this to the real situation of two identical clumps side by side is simple. The distance to the reflecting wall is half of the separation distance if we take the separation distance to be the distance between the clump centers. So the critical separation distance is twice Rcrit.
- I will also start setting up my module to do the problem with interacting bow shocks where one is moving.
Meeting Update 11/05/2012 - Eddie
- still waiting to hear from Anna (rad shock emissions), Francisco (cooling jets), and now Kris (mach stems)
- did some analytic work on the mach stem project following Kris' page: ClumpClump/MachStems
My problem set up will now allow me to test different gammas and different separation distances. ehansen10312012
- going to start working on implementing cylindrical reconstruction for 2.5D ehansen10202012
preliminary mach stem runs
This is a preliminary run for the mach stem problem with a single stationary clump. This is a slice of a 3D simulation. The left boundary is reflecting to simulate two bow shocks interacting from two identical clumps. I am setting up refinement objects around the clump and in a small rectangular region along the left edge where the mach stem should form. Here are some of the parameters:
Object | Density (particles/cc) | Temperature (K) | Velocity (km/s) | Size (100 AU) |
---|---|---|---|---|
Ambient | 100 | 1000 | 0 | 10 x 20 x 10 |
Wind | 1000 | 1000 | 50 | —- |
Clump | 1e6 | 0.1 | 0 | radius = 1 |
I need to do some more detailed calculations to set the rectangular refinement region more appropriately, but I think this illustrates the basic idea of what I am doing.
UPDATE
I did some calculations and found that the critical angle occurs at angles larger than i thought, which corresponds to regions along the bow shock that are closer to the head. So the above image is simulating way too much of the bow shock. To get a mach stem, the reflecting wall should be much closer, so we can actually simulate this at a much higher resolution than the previous run. I also took advantage of some symmetry and am only simulating half a clump.
Meeting Update 10/29/2012 - Eddie
- old tickets: #82, #121
- finally got relative error plots for primitive variables of the radiative shock models ehansen10232012
- started setting up refinement criteria for mach stem problem, will have images up later this week
Comparing Shock Models
I regenerated the images for comparing my shock models (dashed) with Anna's (solld). I am comparing density (blue), temperature (red), and ionization fraction (black). I also made plots showing the percent relative error. The error does not look that good even though the plots of the primitive variables look okay. It just goes to show you how deceiving plots can be, and this is why quantifying differences and errors is so important.
Primitive Variables | Percent Relative Error |
---|---|
![]() | ![]() |
![]() | ![]() |
![]() | ![]() |
UPDATE 1 (10/23 1:00 PM)
I discovered something that might be screwing up the relative error plots. Anna's data does not have a fixed step size while mine does. I set it up so that we have the same number of data points, but Anna's data points are more concentrated where gradients are higher whereas my simulations are a fixed grid. Therefore, visit could be doing the relative error based on the point number rather than the actual distance behind the shock. Perhaps if I use gnuplot it will do some type of interpolation to get a more accurate error plot.
UPDATE 2 (10/23 9:00 PM)
As I previously mentioned, visit could not correctly interpolate the data to produce accurate relative error plots. gnuplot can interpolate data but only as a plotting option. I couldn't figure out how to do operations on the interpolated curves. So I did something different…I fitted polynomials to the data points and then found the relative errors between the polynomial curves. To achieve better fits, I used 10th order polynomials. I have only done this for vs = 30 km/s, but I will do the rest tomorrow. I will also post images of the fits vs. the data so that you're convinced that the fits are okay.
UPDATE 3 (10/25 11:00 AM)
The above method worked pretty well for vs = 30 km/s. The fits get a little oscillatory for temperature and ionization fraction, but they are not too bad. However, the fits for the higher shock velocities did get pretty bad (especially for temperature). So I will have to come up with something else. I think my only option at this point is to write my own little fortran routine to interpolate the data and calculate the errors for me.
UPDATE 4 (10/25 3:30 PM)
Finally got accurate relative error plots. I wrote a little fortran program that resamples/interpolates the data, then calculates relative errors and outputs to a curve file to easily create the plots in visit. I also redid the primitive variable plots using gnuplot instead of visit, so that I could have better control of the tick marks.
Meeting Update 10/22/2012 - Eddie
- learned about geometrical source terms, primarily 2.5D (ehansen10202012)
- learned how to use valgrind and wrote a script to automate testing for memory leaks. This can become part of the test suite in the near future. (ehansen10182012)
- updated the 3DCoolingJets page with some new images/movies. I still need to play around with the volume renderings. The strong cooling case still won't run in 3D, even with the memory leaks patched. Getting this run to work is not crucial, so I will put this on the back burner for now, especially since Shule wants to use the Kraken token this week. I emailed Francisco about getting the parameters to set up a run to simulate his lab experiment.
- found a typo in my H-alpha emission routine. Fixing it helped a little bit, but my data still does not match Anna's data. I am waiting to get a copy of the shock code that she is using so that I can figure out the source of our discrepancies. The following image is from before the typo fix, but this is essentially where I'm at now. v30_p2_f0.01.png
- started setting up a physically realistic run for the Mach Stem project. I still need to play around with the parameters and refinement, but I should have images later this week.
A Brief Lesson in 2.5D
A couple of weeks ago, I learned how to run a 2.5D (axisymmetric cylindrical coordinates) simulation in astrobear. More recently, I learned about geometric source terms in general, and how astrobear actually accomplishes its transformation to 2.5D. This blog post is a brief summary of what I learned and is intended as an introduction for the naive user…
The 2.5D Source Term for Mass Conservation
To illustrate the concept of geometrical source terms, it suffices to look at mass conservation in 2D and 2.5D. As you already know, the mass conservation equation can be written as:
First, we move the flux term to the other side and use a product rule: Now, we can simplify the right hand side using 2D Cartesian coordinates: Now, we can simplify eq. (1) again but this time in axisymmetric Cylindrical coordinates. Axisymmetry just means that there is no dependence on phi, so those terms from the del operator go to 0. So eq. (1) now becomes: This can be simplified further with a product rule: If you compare equations (2) and (3), the meaning of a geometrical source term is realized. The flux terms are essentially the same (just in different directions), but switching to cylindrical coordinates introduces an extra source term Similar derivations can be done for the other conservation equations to get more geometrical source terms. Those can get slightly more complicated when you allow for rotation (velocity in phi direction), and if you have magnetic fields. You could also derive geometrical source terms for other coordinate systems.
How AstroBEAR Handles 2.5D
Switching on 2.5D is quite simple; you just set iCylindrical in physics.data. 0 is off, 1 is for without rotation (aka angular momentum) , and 2 is for with rotation. When iCylindrical is 1 or 2, astrobear treats the x-axis as the r-axis and the y-axis as the cylindrical z-axis. When you look back at equations (2) and (3), you can see why mapping x —> r and y —> z makes a lot of sense. This information is also explained in PhysicsDataExplained.
Now, when astrobear is actually solving the hydro (or MHD) equations, it uses a Runge-Kutta method. Those methods can get somewhat involved, but in our case, they are basically used to solve equations of this form:
where change is a time derivative, fluxes are spatial derivatives, and sources are everything else. These methods do not care if the flux terms are in the x direction or r direction; they are just generic flux terms. So the only thing that astrobear has to do differently for 2.5D is add the geometrical source terms to the dqdt-array (the array that represents the time derivative of q). This step is done in the function SrcDerivs in source_control.f90, and the actual cylindrical source terms are defined in cylindrical.f90. If you look at cylindrical.f90 you will see all of the source terms for all of the variables in the dqdt-array (including the one I derived above). Note that source terms also have to be added to any tracers that might be present, and tracers are handled just like density.
Cylindrical Coordinates in Athena
This is the paper that Jonathan was referring to. It is about implementing cylindrical coordinates in the Athena code. The paper mentions the geometrical source terms that I discussed, but includes MHD and rotation (also, in primitive form rather than conservative). The bulk of the paper is focused on how to do the reconstruction in cylindrical coordinates without changing the Riemann solvers. This definitely looks like something we could do for astrobear. http://adsabs.harvard.edu/abs/2011ASPC..444..260S
valgrind script
I made a fancy bash script to run valgrind. Originally, I was just making this so that I could easily compile and run valgrind on any problem module since I am now the "valgrind czar". However, I have since added a lot of functionality so that this could perhaps become part of the test suite.
Basically, the script will compile astrobear, run valgrind, and then check the valgrind output for memory leaks specific to astrobear. The script will print the tracebacks for all the memory leaks that it finds to a file named astrobear_leaks.out. If that file does contain leak information, then a message will tell the user to check that file. If no leaks are found, that file should be empty, and a message will tell the user that no leaks were found. I implemented many checks and error messages to make the script fairly fool proof. There are also a handful of options the user can specify…
Usage: ./valgrind_leaks [OPTIONS] -m module
specified module must be in the modules directory
Options: -d code_dir
specify code directory (default is present working directory)
-sc —skipclean
skip clean, but still recompiles (use only if you have previously compiled with proper FFLAGS for valgrind, and you have not made changes to Makefile.inc or makefile)
-sm —skipmake
skip compiling completely if you just want to rerun
-sr —skiprun
skip running if you already ran valgrind and just want to check output
-np
choose number of processors to run on (default is 1)
The script is in /home/ehansen/astrobear_GVpp if anyone is curious. If anyone has any questions or comments, just post a comment to this blog post or email me.
update on cooling jets
It was much easier to visualize the strong cooling simulations on kraken. For future reference, anyone that wants to run visit in parallel on kraken should see this web page. Also, in order for visit to read chombos in parallel you have to do this first:
module swap hdf5 hdf5-parallel
So the new images/movies are from the strong cooling case (2D and 2.5D) and can be viewed at 3DCoolingJets. Still working on the 3D version…it appears to be running farther than before, but I have to dig deeper to see if I can find any more memory leaks. At least I am getting good valgrind practice.
Meeting Update 10/15/2012 - Eddie
Cooling Jets
My priority over the last week has been the cooling jet sims, and I am close to being done with all the runs that I had set out to do. Jonathan found a memory leak that was causing the strong cooling runs to fail. With this patched, I was able to run the 2D and 2.5D versions of the strong cooling case. The 3D version is still having some issues, but we might be able to get this one working too. All images and movies can be found at 3DCoolingJets. The volume renderings still need some work.
Cooling Diagnostics
I made cooling length and cooling time fields in fields.f90 in order to calculate the cooling strength chi (cooling length / jet radius). I had previously calculated chi by running several 1D radiative shock simulations which allowed me to come up with an analytic formula for cooling length in terms of shock velocity and ambient density (see ehansen06122012). With the new cooling diagnostics, I can see more quantitatively the effect of cooling near the jet head (see ehansen10122012). I also confirmed, at least for the weak cooling case, that my analytic method for chi is consistent with the cooling diagnostics method.
Comparing Radiative Shock Models
I have been working on comparing my data with Anna's data (Pat's student). The comparison of emission data did not look promising, so I took a step back and compared the hydrodynamic quantities density, temperature, and ionization fraction. I compared 3 different runs which vary only in shock velocity. They all have a preshock density of 100 cm-3 and a preshock ionization fraction of 0.01.
Configuration Script
I worked on writing a configuration script for astrobear. I made one that works as a bash script, but it is somewhat limited. It worked on grass, and it allowed me to compile and run astrobear without having any modules loaded or a Makefile.inc. However, we probably do not want to use this since it is specific to how our library paths (HDF5, HYPRE, etc.) are named. The best configuration script would be an AutoConf script. I am reading up on this and even started writing one, but this will take some time because it is like learning another language. Depending on when we need this to be done, it might make more sense to pass this project to someone else?
Mach Stems
I did not get around to doing much for this project, but I did notice that Kris has a nice project page that I will be referencing. ClumpClump/MachStems
plotting local chi ("cooling strength")
I was a bit skeptical that calculating cooling lengths from 1D radiative shock simulations would directly translate to these 3D jet simulations, so I implemented a field variable to calculate the cooling length on the fly for each cell. Doing it this way results in chi's that are very different. Here is a pseudocolor plot of chi from what I have been calling weak cooling.
I had estimated a chi of approximately 0.4 using my 1D radiative shock simulations, but as you can see here chi can locally be much different. However, you sort of have to take the ambient with a grain of salt. You get extremely low chi's in those regions because nothing is moving, and the material has not been shocked yet. Here is a lineout of chi along the center of the jet.
You can see that the ambient is very low, and then chi does some interesting things as it goes through the bow shock, and then it becomes constant once you enter the jet beam. For most of the jet chi > 1. Here is a zoom in of the chi lineout at the head of the jet.
In this specific region chi does fall between 0.1 and 1, so perhaps a chi of 0.4 is a good approximation after all.
UPDATE
The above plots are not the original post; they are more accurate now.
Also, since plotting hte cooling length of the ambient is kind of strange, I also decided to plot cooling time. Now the ambient looks exactly the same as the jet beam…
the actual 2.5D sims for the cooling jets
So I made a silly mistake last week. The images I showed on Monday's meeting are literally 2D, not 2.5D. I altered my problem.f90 and the .data files accordingly to run actual 2.5D sims. Here are the images for the weak cooling case.
The 2D and 2.5D versions of the moderate cooling case are running right now, and when those are done I will run the 2D and 2.5D versions of the strong cooling case. The 3D versions of both the moderate cooling and strong cooling case will have to be run on kraken, but as mentioned in previous posts I do not know if I will be able to get the 3D strong cooling case to work.
I quickly learned how to set up a 2.5D simulation, and there wasn't much on the wiki about it, so I plan on adding some details.
UPDATE
I added a small section about cylindrical symmetry at the bottom of PhysicsDataExplained.
Meeting Update 10/08/2012 - Eddie
- made the cooling jet simulations work in 2D ehansen10032012
- I didn't like the way the moderate cooling case looked. The jet head doesn't form a good bow shock; it just breaks up into the ambient. I did find a small typo in the problem module, so I'm rerunning that.
- still can't get the strong cooling case to run on kraken. It stops after only one frame dump, so it could be a memory thing. I'll look into this more this week.
- got the problem module set up and ready to go for the mach stem research ehansen10052012
Mach stem simulation set up
I have a module written and ready for the mach stem research that I'll be doing with Kris. This run was sort of just a "proof of concept" to make sure I could do runs like this and that it works. To start, I just need to tweak the parameters and do a set of runs with different gammas.
Cooling Jets - 3D vs. 2D
The cooling jets that I have been working on have all been 3D runs, but the very strong cooling cases require a lot of resolution. In order to get around this, I modified a few things in Martin's jet module to make it work in 2D. I ran the weak cooling case in both 2D and 3D to compare.
2D | ![]() |
---|---|
3D Slice | ![]() |
As expected, they are not exactly the same, so we just have to decide if 2D is good enough. For the very strong cooling case it may be our only option. The 2D run above did run extremely fast (less than a minute on grass).
plasma beta plot for 1D MHD radiative shock
Here is the plot of beta for my 1D MHD radiative shock. I also put density on the same plot, and some constant lines to show where beta = 1.0. So you can see that the density increase starts to slow down when beta goes to 1.0. But I would say that the density profile does not go relatively flat until a lower beta (approx. 0.4 or so).
Meeting Update 09/24/2012 - Eddie
1D MHD Radiative Shocks
The MHD version of my RadShock problem appears to be working, but I am not sure if the output is physically accurate. Red is hydro, blue is MHD.
I am now waiting for Anna (Pat's student) to send me her data so I can compare these MHD runs. I will also be comparing the emission lines.
Post-Processing Emissions
I finished implementing some emission lines into astrobear's post-processing routines. We now have the capability to get [O I] (6300), [N II] (6583), [S II] (6716), [S II] (6731), and H-alpha. Other lines for these species could easily be added, and we also have code for [O II] and [N I] if we want to be able to use those species in the future. Here's an example of what the output looks like for one of my 1D radiative shocks (hydro):
Jonathan also changed a few things in astrobear's post-processing routines, and I updated a bunch of problem modules to conform with these changes. Since AstroBEAR 2.0 is about to be released, I guess we would want this stuff in our new development version (aka AstroBEAR 3.0)?
3D Cooling Jets
Moderate cooling run resubmitted on Kraken, should finish this week. Not sure if we can get the resolution needed for the strong cooling case…
Mach Stems
I will start working on this project tomorrow.
Meeting Update 09/172012 - Eddie
- 3D Cooling Jets - I tested my weak cooling run on alfalfa (8 procs) and bluegene (256 procs). Bluegene should be significantly faster, but for whatever reason it runs less efficiently. I moved to Kraken and ran the moderate cooling case. This run is also fairly inefficient (35%), but since I can use so many procs (4092) it ran halfway in 17 hours. Not sure why it did not run for the full 24 hours. I have not had a chance to visualize the data yet. So it looks like I should be able to finish this on Kraken, but the strong cooling case will probably not be possible at 10 cells per cooling length.
- Line Emissions - I started coding the routines required for emissions post-processing. I am essentially taking the code from bear2fix and putting it into astrobear's processing routines. I have most of the infrastructure set up, and I will be testing SII with my RadShock module. After I have coded everything we will be able to get synthetic observations for OI, OII, NI, NII, SII, and H-alpha.
- Mach Stems - A new project that Adam, Pat, Kris and I talked about last week. With all my other things going on, I am guessing that I will not even get started on this until sometime next week.
refinement object in jet sims acting funny
I've been trying to use Jonathan's new refinement object in the 3D cooling jet sims. It's set up so that a refinement shape (in my case a rectangular prism) moves across the grid at a certain speed. Everywhere inside this shape is refined. The shape moves at roughly the same speed as the head of the jet, so that the head of the jet is always refined. However, based on what I see from the mesh in visit, the sim is not always refining everywhere within my shape.
Here are the relevant movies:
I've also tried turning on gradient based refinement within the shape but this looks even worse:
UPDATE
Sorry, the above movies were sliced at the wrong spot. Here is what they should look like:
Refinement in Shape | Density | Density w/ Mesh |
---|---|---|
Gradient-based | rho_grad1.gif | rho_gradmesh1.gif |
These updated movies are on a log scale.
New redirect macro
I installed a new macro called redirect which is very convenient for setting up redirect pages. You just put this at the top of a wiki page:
[[redirect(OtherWikiPage)]]
So when you try to access that wiki page you'll be redirected to OtherWikiPage.
I mainly wanted this for the AstroBearTesting page. If you recall, Jonathan had updated the Tests macro to display links to each test's wiki page. This link is based on the test's name, so for example the link for the FieldLoopAdvection test takes you to TestSuite/FieldLoopAdvection. However, this test also has a restart test which is called FieldLoopRestart. FieldLoopRestart does not have its own test page on the wiki so if you click on TestSuite/FieldLoopRestart you just get redirected to the FieldLoopAdvection test page.
This will be helpful for other tests that have multiple versions but just one test page.
Meeting Update 08/07/2012 - Eddie
The quarter grid was pretty simple to set up…I changed the domain and boundary conditions in global.data. Also, a couple things for the shape refinement. Run A took just under an hour on 8 procs on alfalfa (those images are below). Runs B and C will need bluegene due to their very high resolution.
The resolution had to be altered a little bit, so just as a reminder here are the parameters for my runs:
Run | Density | Chi | Base Grid | AMR Levels | Cells per Jet Radius | Cells per Cooling Length |
A | 10 | 0.358 | 48 x 12 x 12 | 3 | 32 | 11.47 |
B | 100 | 0.0358 | 56 x 14 x 14 | 6 | 298.67 | 10.70 |
C | 1000 | 0.00358 | 68 x 17 x 17 | 9 | 2901.33 | 10.40 |
Technically, the cells per jet radius increased for B and C in order to keep the proportions the same from sim to sim, but the overall domain is smaller so these should run faster now.
I need to spend most of my time studying for the prelim now, so I don't expect to be getting much research done for the rest of August.
Meeting Update 07/24/2012 - Eddie
Used Jonathan's new refinement objects for my 3D jet simulations. Here's a comparison of how it was refining before to how it refines with the shape object that I defined in the problem module:
Old Method | New Method | |
Density | ![]() | ![]() |
Mesh | ![]() | ![]() |
I'm also working on setting up a quarter grid simulation which should significantly speed things up.
In a meeting last week, we decided to keep the density contrast (eta) equal to 1 for the runs we want to do. Now, increases in jet and ambient density result in stronger cooling throughout the simulation. Here is a table for proposed runs:
Run | Density | Chi | Base Grid | AMR Levels | Cells per Jet Radius |
A | 10 | 0.358 | 48 x 24 x 24 | 3 | 32 |
B | 100 | 0.0358 | 54 x 27 x 27 | 6 | 288 |
C | 1000 | 0.00358 | 66 x 33 x 33 | 9 | 2816 |
Resolution is based on ~10 cells per cooling length.
using particle refinement for jet sims
So Martin showed me how I can set up a particle to have better refinement for my 3D cooling jet sims. The old method was just refining along momentum gradients. Now I am just refining a box around the jet head. Run C, which is the density matched (eta=1) jet is running way faster. It took bluegene 3 days to do 80 frames on 256 procs with the old refining method. With the new refining method, it has taken grass less than 2 hours on 8 procs. Here is a snapshot of the mesh with the new particle refinement:
The next step is to only simulate a quarter of the grid since everything is symmetric. This will make the runs even faster.
Meeting Update 07/17/2012 - Eddie
I have some partial runs of my 3D jet simulations. They are all 256 processor jobs on bluegene with rev 950. Here is that chart of all the runs I want to do:
run | n jet | eta | chi | required x res | base grid | levels | actual cells per cooling length | cells per jet radius |
---|---|---|---|---|---|---|---|---|
A | 0.1 | 0.01 | 0.0015 | 78,335 | … | … | … | … |
B | 1 | 0.1 | 0.034 | 3495 | 56 x 28 x 28 | 6 | 10.15 | 298.67 |
C | 10 | 1 | 0.358 | 335 | 48 x 24 x 24 | 3 | 11.46 | 32 |
D | 100 | 10 | 1.367 | 88 | 48 x 24 x 24 | 3 | 43.74 | 32 |
E | 1000 | 100 | 2.428 | 50 | 48 x 24 x 24 | 3 | 77.70 | 32 |
Run A was not even attempted, run B was taking too long so I have to figure out a better way to refine. Martin showed me how I can use a particle to refine only the jet head. I need to play around with this to make sure I get it right for each simulation. In the meantime, I have submitted restarts for runs C, D, and E to finish them up.
It's kind of tricky to compare the different runs on the same color scale since they are all in slightly different density regimes. Here are snapshots of C, D, and E at the same time:
Run | Chi | Image |
---|---|---|
C | 0.358 | ![]() |
D | 1.367 | ![]() |
E | 2.428 | ![]() |
Run C ran the longest, so here is a video (the color scale here is different from the color scale above):
The majority of my time is spent on studying for the prelim now.
Meeting Update 07/10/2012 - Eddie
- My jet sims failed to run on bluegene because mpi did not like my astrobear executable. I was having some trouble compiling in bluegene, but changing some lines in my .bashrc seemed to do the trick. The jobs have been resubmitted, so now I wait.
clarifications on resolution for my 3D cooling jet runs
Here is an updated table from my previous blog post. Again, run A was not submitted to bluegene since it requires a lot of resolution. The required x res is the fixed resolution that I would need to guarantee 10 cells per cooling length.
run | n jet | eta | chi | required x res | base grid | levels | actual cells per cooling length | cells per jet radius |
---|---|---|---|---|---|---|---|---|
A | 0.1 | 0.01 | 0.0015 | 78,335 | … | … | … | … |
B | 1 | 0.1 | 0.034 | 3495 | 56 x 28 x 28 | 6 | 10.15 | 298.67 |
C | 10 | 1 | 0.358 | 335 | 42 x 21 x 21 | 3 | 10.02 | 28 |
D | 100 | 10 | 1.367 | 88 | 44 x 22 x 22 | 1 | 10.02 | 7.33 |
E | 1000 | 100 | 2.428 | 50 | 50 x 25 x 25 | 0 | 10.12 | 4.17 |
… | … | … | … | … | … | … | … | … |
C' | 10 | 1 | 0.358 | 335 | 48 x 24 x 24 | 3 | 11.46 | 32 |
D' | 100 | 10 | 1.367 | 88 | 48 x 24 x 24 | 3 | 43.74 | 32 |
E' | 1000 | 100 | 2.428 | 50 | 48 x 24 x 24 | 3 | 77.70 | 32 |
I also figured out that if I want at least 32 cells per jet radius, then I need a minimum fixed x resolution of 348 which corresponds to a base grid of 48 x 24 x 24 with 3 levels of AMR. Unfortunately, when I set up these runs I was only considering cooling length to set the resolution, so runs C, D, and E are below this minimum resolution based on jet radius although C is very close.
The primed runs are new ones that meet the minimum resolution requirement. This essentially puts "extra" resolution on the cooling length. Fortunately, the old runs were still waiting in the queue, so I updated their global.data files to have the higher resolution. So C, D, and E were successfully replaced by C', D', and E'. I'll have to take a look at the results when I get back from vacation.
Meeting Update 06/26/2012 - Eddie
I ran some very low res sims for the 3D cooling jets to make sure the problem worked and I had the dimensions right. The left image is what it looked like at frame 70 (100 frames total). I made the domain a bit bigger in all directions so that frame 100 now looks like the right image.
![]() | ![]() |
movie |
---|
Since these test sims looked okay, I submitted 4 jobs to bluegene this morning (Monday, 6/25). The wait in bluegene looks like it's about 3 days, so I guess I'll have to run the analysis when I get back from vacation.
With the ambient density held constant at 10 cm-3, I change the jet density to change the cooling strength chi. When chi is really small, this implies very strong cooling, and it becomes difficult to resolve the cooling region. All runs have a jet velocity of 250 km/s and a jet radius of 2e16 cm. Here are the varying parameters for 5 different runs:
n jet | eta | chi | required x res |
---|---|---|---|
0.1 | 0.01 | 0.0015 | 78,335 |
1 | 0.1 | 0.034 | 3495 |
10 | 1 | 0.358 | 335 |
100 | 10 | 1.367 | 88 |
1000 | 100 | 2.428 | 50 |
The required resolution in the x-direction is based on a desired resolution of 10 cells per cooling length. You can see that the low chi runs require a lot of resolution, but there are ways around this. I can probably come up with a clever way to only track the jet head and region immediately behind it. For now, I did not submit a job for the first run in the table (njet = 0.1), and I threw some AMR levels at the other runs.
So here is my plan for when I get back:
- analyze these runs that I have submitted
- code something to force refinement only where I need it (this will help with low chi sims)
- run some sims with Francisco's parameters (Adam, can you forward this info to me please?)
updated documentation for my projects
My project pages were getting a little crowded, so I separated them and organized them in a logical manner on my page.
The work that I've done on 1D radiative shocks is now on a separate page. It has relevant equations and plots for the different types of cooling. The last plot on that page is an interesting one that I haven't shown before.
This is where I explain non-equilibrium cooling. I go into detail about what was changed from AstroBEAR 1.0 to 2.0. A lot of this information corresponds to ticket #147.
This Z cooling page is now mainly just about the implementation into astrobear. I also created ticket #224 for this, and it is part of the milestone called Implementing Multi-Physics.
This page contains remnants from my 1D pulsed jet simulations that I attempted last fall. I'm not sure if I'll be picking this project back up, or just moving on to something else for now.
This is the new page where I will post things for my 3D cooling jet simulations. I will also put the stuff about cooling length and cooling strength here, but the page is blank at the moment.
Meeting Update 06/19/2012 - Eddie
- some issues with Z cooling resolved, everything works now ehansen06082012.
- compared my shock with Pat's shock ehansen06112012.
- figured out cooling length and cooling strength for 3D jet simulations which are based on the Blondin paper ehansen06122012.
- testing Martin's jet module locally, will move to bluegene for initial set of runs
- also started reviewing for the evil prelim
pre-run analysis for 3D cooling jets
Cooling Length
I've been following this Blondin paper to prepare for the 3D cooling jet simulations. http://adsabs.harvard.edu/abs/1990ApJ...360..370B. Section III is especially important because this is where the cooling length is defined. The analytic expression is there, but for analysis Blondin uses the numerical approximation from 1D simulations. I am therefore using my 1D radiative shock module to come up with my own formula for the cooling length.
Blondin's cooling length formula is going to be different because he did not use the DM curve. His cooling curve is based on a nonequilibrium ionization calculation by Kafatos (1973). Hartigan, Raymond, and Hartmann (1987) also have a different cooling curve and therefore a different cooling length formula. Blondin's is proportional to vs4 and Hartigan's is proportional to vs4.67.
The goal is to get a formula such that
where dcool is the cooling length, na is the ambient number density, and vs is the shock velocity. b and a are the free parameters to be solved for. The formulas from Blondin and Hartigan do not match my data very well.Source | b (1016) | a | Sum of Error2 |
---|---|---|---|
Blondin | 4.5 | 4 | 125.96 |
Hartigan | 1.80 | 4.67 | 121.39 |
Mine | 3.51 | 3.2 | 10.42 |
Here, vs is in units of 100 km/s to give dcool in cm. Obviously, my formula is going to do the best because it's based on my data.
It is also important to note that the cooling lengths that I am calculating use a floor temperature of 8000 K just like Blondin. Hartigan's cooling length allows for cooling down to 1000 K.
Cooling Strength
Now, the cooling strength can be defined as:
where rj is the jet radius, vj is the jet velocity, and is the density contrast (jet/ambient). The idea is to keep rj, vj, and na fixed. Then, you change the jet density and therefore changing in order to test different cooling strengths.For some reason, I cannot reproduce the chi's that are reported in the Blondin paper. I've looked over my work several times, and I cannot find any errors. Maybe the paper has a calculational error or a typo somewhere. When I calculate chi for their "standard run" using their formula I get 0.27, but they claim that it should be 0.55. Not sure whether I should ignore the discrepancy and just do my own calculations, or if this needs to be resolved before moving on.
UPDATE
After doing more calculations for different sets of parameters, I realized that all of my cooling strengths were exactly ½ of the cooling strengths in the Blondin paper. It seems that Blondin either forgot a factor of 2 in their formula for chi, or all of their cooling strengths should really be ½ of what they reported.
Z cooling plots
Here are the other two important plots for comparing my shock module with Pat's model…
The overall shape/trend is good. Not sure why Pat's ionization fraction is a bit higher.
Density varies greatly towards the end of the cooling region. This could be because Pat's model keeps cooling past 2000 K where mine does not. Z cooling stops at 2000 K and NEQ takes over. The cooling rate for NEQ at those temperatures is very low.
Also, here is temperature plotted on a log scale. It's easier to see the differences at the lower temperatures this way.
ZCooling Works!
Thanks to Jonathan, ZCooling is finally working. He found the last tiny error in the cooling routines, and it made the shock steady. My shock is now even closer to resembling Pat's, and I believe I understand why it's different. The issues I was having with implementing a variable mu are fixed now too. So this is the full result…Zcooling, varying mu, MHD should work but these runs are not MHD. Basically, the difference between NEQ and Z cooling is that Z will give you much more cooling below about 10,000 K or so.
Pat's is in black, mine is red
There are two discrepancies here:
- The big one is the difference between our immediate post-shock temperatures. I believe this can be accounted for by mu. Pat's model uses a bunch of heavier species, and I calculated his mu to be about 1.26 where mine is approximately 1.0. If my simulations had an initial mu of 1.26 my post-shock temperature would be much closer to Pat's. I should be able to resolve most of this by implementing He into my module.
- My model has a bit of a jump near 2000 K. There is also a less noticeable one at 16,500 K. This is due to the way I am using Z cooling. NEQ cooling is a combination of a few things: hydrogen excitation, ionization, recombination, and metal excitation. Z cooling is supposed to replace the metal excitation component. Right now, Z cooling and NEQ metal excitation switch on or off at the upper and lower limits of the Z cooling table. This transition just needs to be "smoothed out".
UPDATE
I fixed the two discrepancies. I could not get He ionization/recombination to work, but just adding the right amount of neutral He increased mu and brought my solution much closer to Pat's. I also "smoothed out" my solution with what I call a parabolic weight function. I don't know if that's what the math world calls it, but it worked. I will play around a little more with the He ionization/recombination, but this is already looking good enough. I'll also post the density and ionization fraction plots soon. Again, my run is red and Pat's is black.
There is probably a good reason why my run seems to flatten out at 2000 K and Pat's keeps cooling a little more. The Z cooling table that I have stops at 2000 K, and so NEQ cooling takes over below that point. I suspect that whatever Pat uses is better at cooling at such extremely low temperatures.
Meeting Update 06/06/2012 - Eddie
First, see my recent blog post ehansen06052012.
NEQ cooling with a varying mu is steady now, but the shock for Z cooling is still moving. I haven't been able to find the error yet, but I'm still working on it. Here's a temperature movie of both NEQ and Z. NEQ is dashed black line, Z is solid red line. NEQ_Z_temp.gif
I also started reading up on the jet simulations with cooling, so that I know what to do for these upcoming 3D runs. I have Martin's set up that he used for CRL618 and the Blondin paper http://adsabs.harvard.edu/abs/1990ApJ...360..370B.. Basically, the idea here is to run 3D jets with different cooling strengths. The cooling strength is defined in the aforementioned Blondin paper. For their different runs, they have some cooling curve (similar to DM) and they change the jet density in order to change the cooling strength. My question is…do I want to replicate the stuff in this paper by using their parameters, or should I start with the default parameters that Martin has used in his CRL618 runs?
Issues with NEQ and Z Cooling
I haven't been able to get steady shocks for these types of cooling yet. I originally thought that it had to do with mu (mean molecular weight). I need to calculate mu for every cell, because it is always changing. Heavier species will increase mu, and ionization will decrease mu. However, I just realized this morning that the shocks are not steady for a different reason. I went back to the beginning and just set mu = 1, and the shocks move. I'm guessing that there's a scaling issue somewhere since that has often been the case when I am close but not quite steady. There could also be a typo in the way temperature is being calculated. I'll look into this further and update this post later today.
Why mu is important
The following equation can be derived from the Rankine-Hugoniot conditions in the strong shock limit:
This shows that post-shock temperature is proportional to mu. Temperature in general is actually always proportional to mu because of the ideal gas law. We are just so used to having mu = 1.0 for an ideal hydrogen gas. For my pure hydrogen runs, mu is just 1.0, but Pat's runs have other species. The aforementioned post-shock temperature is immediately after the shock interface and before ionization occurs. Therefore, since the gas has no pre-ionization, the mu for Pat's runs is higher.
Let's look at the temperature for both runs (Pat's is in black, mine is in red):
You can get a rough estimate of what mu should be in Pat's runs by just looking at the difference in temperature. I suspect that Pat has an initial mu of approximately 1.22.
Here is how mu is calculated in the cooling routines:
mu = (nvec(iH2)*muH2 + (nvec(iH)+nvec(iHII))*muH + (nvec(iHe)+nvec(iHeII)+nvec(iHeIII))*muHe)/npart
nvec is the number density vector, and it contains the number densities for all the different species. npart is the total number density of particles which includes all the electrons as well. This is why mu significantly decreases with ionization. For example, a fully ionized gas of HII would have mu = 0.5 instead of 1.0 for a neutral gas of all H. In other words, you double the number of particles without adding any mass (electrons are treated as mass-less). Recombination has the opposite effect as it reduces the number of electrons.
So in general, mu is important in calculating temperature. Both for the immediate post-shock temperature, and the decreasing temperature within the cooling region.
Here is what mu looks like throughout the problem domain:
Notice that mu starts off a little less than 1.0. That is because there is a trace amount of pre-shock ionization. Mu then drops rapidly (to about 0.65) after the shock due to ionization,and gradually increases due to recombination.
Meeting Update 05/29/2012 - Eddie
Not much to report this week. Last week, I spent the majority of my time wrestling with my radiative shock module and the Z cooling routines. I'm working on implementing a mu (mean molecular weight) that varies from cell to cell. Ionization decreases mu, recombination and the presence of heavier elements such as He increases mu. However, I'm not even messing with He yet. I have been able to set up an initial profile with a variable mu, but have not been able to keep the shocks steady. I'm not exactly sure what is going on, so I'll have to keep working on it.
In other news, I closed on my house last Friday, so I've been extremely busy moving, cleaning, fixing things, etc. Turns out that buying/owning a house is a ton of work, but the "16-hour-manual-labor-work-days" will be over soon.
comparing my rad shock with Pat's
I translated my data so that it starts at the shock interface just like Pat's data. Also, I believe Pat used a small magnetic field, so I added that as well. The field is weak, so it doesn't really do much. These runs are supposed to be using Z cooling which includes Pat's tables. However, they don't look much different from the previous NEQ runs.
red is my run, black is Pat's…
Temperature looks okay, cooling length is close, Pat's postshock temperature is higher for some reason. This could be due to a different mean molecular weight. In the strong shock limit, postshock temperature is proportional to mean molecular weight. So if Pat's run has some helium or other species which make the mean molecular weight closer to 1.22 instead of 1.0 then that would explain this discrepancy.
Ionization Fraction looks good. The initial increase in ionization agrees very well. The fact that my data does not get as high as Pat's can be explained by the difference in postshock temperature.
Our simulations differ mostly in the density plot. I think this is just a code thing. They both agree very well in regards to the postshock value (almost 4 times preshock which makes sense). My code stops compressing when it stops cooling, and cooling stops once the temperature goes back to the ambient/preshock temperature. Pat's run keeps cooling past the ambient/preshock temperature, so it keeps compressing.
Meeting Update 05/22/2012 - Eddie
Took me a little bit longer than I expected to sort through Pat's data and plot it nicely in Visit with my data. I can compare density, temperature, and ionization fraction. So far I've only looked at temperature because it doesn't look quite right. They are certainly within the same ballpark though…
red is my simulation, black is Pat's data
Meeting Update 05/15/2012 - Eddie
Nonequilibrium cooling (NEQ) aka Modified DM is working with my radiative shock module. This type of cooling is more efficient than regular DM, so the problem had to be physically scaled down a little bit in order to see the cooling region.
Here is the expected temperature profile NEQ_temp.gif
Here are the profiles for various densities (total, H, and HII): NEQ_rhos.gif
The problem was that I didn't realize that the source term modules currently use q in primitive form. The ionization routines assume q is in conservative form. So all I had to do was a few little conversions and everything worked. The routines will need some more modification and testing to use other available species: H2, He, HeII, HeIII. For now I'll just stick with H and HII, and I'll move onto Z cooling to see what problems that creates.
UPDATE
I adjusted how ionization/recombination is handled, so that when cooling turns off the species densities can still change. This leads to more physically accurate results: NEQ_rhos1.gif
Meeting Update 05/08/2012 - Eddie
Sorry, not much to look at today, just a brief summary of what I'm working on:
- Still working on getting ionization to work in my radiative shock module. I feel like I'm getting closer though. I'm no longer getting any compiler or runtime errors. The output just looks off. It may just have to do with a postshock jump condition for ionization? I will contact Pat on this.
- I'm also working on a revision of the code that will take care of the BScale typo (#198) and protections for initialization (#171). The first ticket is trivial. The second ticket came back with an error when Baowei tried to run the test suite, so I'll take another look at that.
- Lastly, I'm going to start running 3D jet simulations with different cooling strengths with Martin's CRL618 setup.
UPDATE
It appears that NEQ cooling (modified DM) yields much stronger cooling than just DM. I shrunk the problem domain by a factor of 40 to get a shock profile similar to the DM simulations. However, the shock is still not steady. The simulation seems to develop a forward and reverse shock. So something is obviously not right, I'll keep working on it.
Meeting Update 04/24/2012 - Eddie
I successfully added MHD to the 1D radiative shock simulations. The magnetic field is only in the y-direction (perpendicular to the fluid flow). This was initially trickier than I thought it would be. Having magnetic fields changes the shock jump equations and the ODE within the cooling region. I made a lot of changes to my problem.f90 so that the module as a whole could easily switch between hydro and MHD. As expected, the presence of the magnetic field lengthens the cooling region. Here is a movie that compares the temperature profile for hydro and MHD (both simulations use DM cooling). DM_compare_temp.gif
I'm still working out some kinks in using the next type of cooling (NEQ cooling). The tricky part here is now I have to create initialization profiles for hydrogen H and ionized hydrogen HII (I'm not even considering molecular hydrogen, or helium and its ionized species). The cooling routines calculate these ionization and recombination rates, and these can be used to update my initial quantities to create an initial profile. These rates have to be self-consistent with the cooling rate that NEQ cooling gives. So the problem essentially becomes solving two different ODEs simultaneously. I think I have this done correctly, and I am now just running into issues with the way the main cooling routine interfaces with NEQ cooling.
Meeting Update 04/10/2012 - Eddie
DM cooling added to Radiative Shock module
Due to the initial conditions of the problem, the immediate post-shock temperature is near the peak of the DM cooling curve. As a result, the simulations that use DM cooling cool much faster than the analytic simulations where the rate ~ T2. I found that, on average, the DM rate is approximately 250 times greater than the analytic rate within the cooling region. So for the DM simulations, I scaled the length down by a factor of 250. This resulted in a cooling region that was about 10% of the problem domain as it was in the analytic case.
However, even with the same resolution analytic and adiabatic simulations(~40 cells across the cooling region), the DM simulation was not very steady. I ran the DM simulations 10x longer than the analytic and adiabatic simulations that we looked at last week. I also looked at the analytic simulation at this longer runtime, and it was still steady. The snapshots below are all at time 0.12 (the final runtime of the analytic and adiabatic runs). Here is the DM 400 cell simulation:
I therefore doubled the resolution to ~80 cells per cooling region. Here is the 800 cell simulation:
I again doubled the resolution to ~160 cells per cooling region. Here is the 1600 cell simulation:
It looks like I am again running into resolution issues. However, I think that the interpolation of the DM cooling curve might be to blame?
Update
Decreasing the ambient velocity (aka shock velocity) solved the problem. The ambient conditions used for analytic cooling only work for analytic cooling. They lead to an immediate post-shock temperature that is beyond the peak of the DM cooling curve. Therefore, with DM cooling, these radiative shocks would be unstable. Reducing the shock velocity reduces the post-shock temperature. I found that changing the length scale as previously mentioned, and changing the velocity from 100 km/s to 80 km/s keeps the post-shock temperature below the critical value for stability (which is approximately 2e5 K). Here is the simulation for DM cooling (400 cells):
Meeting Update 04/03/2012 - Eddie
The Radiative Shock simulation works!
I've updated my wiki page with a description of the problem, some relevant equations, plots, and movies. Scroll to bottom of zcooling for details.
Radiative Shock Simulation
In this problem, you set the ambient conditions and use shock jump equations to solve for initial post-shock conditions. Then, the initial post-shock conditions become your initial boundary conditions for the ODE inside the cooling region. I've been using Mathematica to solve this ODE for the pressure profile. The profiles for other important hydro quantities follow from this. The analytic form that Mathematica gives is not pretty, but when plotted it looks qualitatively correct.
For a while, I was struggling to get my profiles to look like those in Delamarter '01. I realized that the solutions for the pressure, temperature, etc. profiles in the cooling region are strongly dependent on gamma. So I made gamma a free parameter again, instead of just using 5/3. One problem is that Delamarter makes no mention of what value he used for gamma. Also, he does not give values for the initial post-shock quantities, so there is no way to calculate what he used for gamma. The best I can do is guess by looking at the plots. I think I get close to the initial post-shock values with a gamma of 1.25. The profiles look qualitatively correct, but my cooling region is about 10 cells (2.5*1016 cm) smaller than that of the Delamarter simulation.
Another issue is that my solutions are not steady on the timescales stated in Delamarter. It should be relatively steady for at least 4000 years, but my solutions are not at all.
I see four possible reasons for these discrepancies:
1) There is something wrong in AstroBEAR
2) There is something wrong with the Delamarter simulations/calculations
3) The analytic formulas from Mathematica are wrong+
4) I am still not using gamma correctly++
+Note that I solved the differential equations in Mathematica analytically and numerically…the results are the same.
++Another thing that I haven't tried yet is using two different gammas. One gamma for the cooling region, and a different gamma elsewhere.
I have not yet tried using my own numerical ODE solver, so maybe that's the next step. I'll post some equations and plots before our next group meeting.
Meeting Update 03/13/2012 - Eddie
Cooling Curves
I made a blog post last week with a plot of cooling curves to compare Pat's tables with Dalgarno-McCray. It's interesting how different the curve can be based on the electron density and ionization fraction.
My Library
I went through a bunch of old emails from Adam and Pat, and I made a list of papers that I want to take notes on to become more familiar with them (for my library on the wiki). It seems like a long list that could get longer, so I will probably focus more on the key papers. Adam, are there any papers in this list that you feel shouldn't be there?
Cooling Project
I updated my version of astrobear to try the changes that Jonathan had made with how the pressure is handled by the cooling source term. It changed a little bit, but there are still discrepancies between the cooling and adiabatic simulations. It appears as though it helped, the post-shock temperature in the cooling cases went up a little bit, but it's hard to tell due to the oscillatory nature of looking at frames that are several time steps apart.
So I've moved on to writing a 1D Radiative Shock module, and I will most likely have results from that next week. Just as a reminder, this follows what was done in AstroBEAR 1.0 (1D Radiative Shock test page) and Ch. 4 of Delamarter '01 (pdf). This module could also become part of the test suite if we want.
Cooling Curves
I made this plot to look at how the Z cooling curve compares to the DM curve. Z cooling which, as a reminder, uses Pat's tables is highly dependent on the electron density ne and ionization fraction Xh. Here, I only plotted 2 possibilities: low ne, low Xh and high ne, high Xh. By "low" and "high" I mean the minimum and maximum of Pat's table.
Also, it is important to note that outside of the temperature range of Pat's table, the cooling routines currently revert back to neqcooling. It's more difficult to make a plot of this since there are no tables for this type of cooling. The cooling rates from ionization and recombination are calculated on the fly for neqcooling.
Meeting Update 03/07/2012 - Eddie
finished resolution study
I've worked on completing the results for my 1D jet runs. The zcooling page now has data tables for the post-shock temperature for Z, DM, and no cooling. Here is a new plot that contains data from all the runs:
lineouts to look at shock structure
I've made some lineouts of density, velocity, pressure, and temperature to get a better understanding of how the simulation is interpreting the shock structure. Every image is taken at the middle of the simulation for the highest resolution run. The solid black line is density, dashed blue is velociity, dotted green is pressure, and dash/dot red is temperature. Everything is in computational units, and it's plotted on a log scale. The structure is more difficult to see in the cooling cases, so I included some magnified images for those two cases.
No Cooling
DM
Z
try new changeset, run old astrobear setup
This week I plan on using the newest changeset that Jonathan just pushed. It might have some affect on my simulations. Also, I will be making a new simulation that follows what was shown on the old astrobear website…http://www.pas.rochester.edu/~bearclaw/delamarterplots.html.
important papers
I've also been doing more reading. I think I will make a brief "library" of the most important papers on my user page similar to what Erica did. u/ehansen
Meeting Update 02/21/2012 - Eddie
I have been studying the resolution dependencies on the cooling routines. It seems as though I need a lot of cells per cooling length in order to achieve post-shock temperatures close to analytic predictions. However, I never quite reach the post-shock temperature that I expect. So I don't know if I am still under-resolved. This seems unlikely, because the post-shock temperature appears to be converging on some value, just not the value I expect. Maybe I am missing a fundamental difference between cooling and not cooling that I haven't accounted for.
I have kept diligent notes on this for the past couple of weeks. You can see more details on the zcooling page. The bottom subsection titled 1D Jet Simulations might be moved to a more permanent project page in the future.
Meeting Update 02/07/2012 - Eddie
My 1D Jet simulation now runs literally in 1D with Jonathan's recent 1D implementation. This makes the simulations much faster ~1 minute on alfalfa. I have run the simulations with the new cooling routines both with and without magnetic fields. I'm now working on getting the emission maps for H-alpha and SII. However, the emission routines in bear2fix are proving to be more troublesome than we initially thought.
I get results for H-alpha that look ok. I'm not sure if these are exactly right, but at a glance they make some sense. I am not getting anything for SII. So far I have traced the problem all the way back to how it calculates the temperature. It seems that some of the chombo information is not being read in properly. This is likely due to the chombo being in 2D, but the simulation and many of the arrays being allocated for 1D. Or there is some other bug that I can't locate. All I know for sure at this point is that TempScale is wrong. I have set TempScale to be 1000, but bear2fix thinks it is some crazy number ~10-4. This makes input temperatures way too low to get any emission from anything.
Meeting Update 01/24/2012 - Eddie
- Google calendar up and running. So far it only has the Astro lunches and colloquia, our weekly group meeting, PAS colloquia, and journal club. Any one in our group can create/manage the events on the calendar. What else do we want on this calendar? The AST462 course (not really a group thing, but it does affect Adam, Erica, and I), other weekly meetings, upcoming conferences? Leave comments! Let's use this blog as a discussion as it was originally intended!
- Speaking of blog comments, did anyone ever look into email notifications, or do we even want them? It seems that there is a plugin for this: http://trac-hacks.org/wiki/FullBlogNotificationPlugin. Should I take a stab at this or leave it for someone who is more familiar with trac and how plugins are handled?
- Zcooling routines have been implemented. I'm now in the process of testing them with my 1DJet simulation. I'll try out Jonathan's new 1D implementation as well. The documentation on Zcooling is coming along. See zcooling for more details. The table reformatting stuff will be important for anyone who gets this code in the future. The section on the interpolation that I used will be added soon.
- I also started documentation on the Wind object, since I made some changes to it in order to use it in my 1DJet module. See the wiki page WindObjects.
UPDATE : Jonathan did look into email notification for blog comments a while ago. The Announcer plugin is supposed to support email notification for blogs in the same way it does for tickets. However, Jonathan could not get this to work. I looked into it today as well and found the same things that Jonathan did. The solution is to use the separate FullBlogNotification plugin. This is a standalone plugin that does its own notification without going through Announcer.
The way notifications work is the wiki sends an email to username@pas.rochester.edu. Announcer allows each user to specify an email address in their preferences on the wiki if they don't want to use this pas address. In short, this feature will not work for blog comments.
Meeting Update 01/17/2012 - Eddie
I have a different simulation set up that is now more suitable for what we will want for the HST proposal and further research. It is effectively a 1D pulsed jet running through an ambient medium. The jet is pulsed via velocity perturbations. I added perturbation parameters to the wind object, so anyone that wants to simulate something with a pulsed inflow might be able to use that. It currently supports square wave and sine wave perturbations. I plan on adding a random perturbation option as well. The 1D jet simulation also has a magnetic field perpendicular to the propagation, and cooling can be turned on. Here's an image of what I have so far on a run with no cooling:
I've been reading up on interpolation schemes. Tomorrow I will dig a little deeper into the cooling routines, and really start making progress on adding Pat's cooling tables. Hopefully, next week we'll have this new form of cooling implemented and the aforementioned simulations running.
Meeting Update 01/10/2012 - Eddie
The nonequilibrium cooling routines are coming along. This has also been called the modified Dalgarno-McCray cooling. The files have been ported over from astrobear 1, and the necessary changes have been made to make it work in astrobear 2.
I used Bin's radiative instability module as a quick check to see if NEQ cooling was working. You can tell from the output that it is definitely cooling, and it is cooling more than DM cooling. I started a new wiki page on this: u/neqcooling. So far this page just contains some visit plots from the radiative instability module.
I am now working with a different module to just run a shock through an ambient material. This set up will make it easier to see the post-shock material, because the radiative instability module will often produce oscillations. After I determine that this is accurate, I can start working on the 6th type of cooling which will include Pat's tables, which he has already sent to me.
Non-equilibrium Cooling
The non-equilibrium cooling routines appear to be working. I made a new page where I will document things. So far this page only has some initial plots to look at for the 1/6/12 conference call. u/neqcooling
Meeting Update 12/13/2011 - Eddie
Work on porting the non-equilibrium cooling routines to astrobear2.0 is under way. I've made some of the obvious changes that were required and am now working through the gritty details.
In regards to the test suite, I've been looking into the OrbitingParticles problem again because the images don't look right, and it's not running correctly any more. When I ran it on alfalfa doing a separate compile and run it worked fine. However, it did not run correctly with the buildproblem command. Hopefully, I'll have this resolved before Friday's demonstration, but we still have a handful of other tests that are working just fine if I don't.
I also just started converting my own hydro code to 2D, but I doubt this will be done any time soon. This is just a side project that probably won't receive much attention for a while.
Meeting Update 12/06/2011 - Eddie
Not much to say about the new cooling procedures. I'm waiting to hear back from Kris on that.
I have done quite a bit with the test suite. See ehansen12022011. I also updated all of the tests' wiki pages with Matt's new plugins.
Several test suite updates (astrobear changeset 701 + 702)
I finished fixing the RadiativeInstability files so that they now only produce 1 chombo at the correct computational time. Then, I updated the reference chombos that are in the /data/tests directory on Clover. I also noticed some errors while doing the post processing on this problem. I fixed the bear2fix.data files in astrobear and the clover directory.
Another thing that I noticed while doing the chombo comparison for RadiativeInstability was a relative error of infinity! I realized that the y-momentum is zero for this problem. Therefore, I rewrote bear2fix's i_compare.f90 to avoid this divide by zero issue. It now replaces an infinite error with one of three things:
1) zero error if the new chombo also has zero momentum component
2) an error equal to the new chombo's momentum component if the new momentum component is very small (which would lead to an error within the 'successful' tolerance)
3) an error relative to the average of the reference chombo's momentum components (which could lead to a 'success' or 'fail')
Lastly, the OrbitingParticles problem was added to the test suite. Producing a useful image was the most difficult part here. I "tricked" bear2fix into plotting the sink particle potential. To do this I told bear2fix to plot By. The particle potential is stored in the same location in the q array that By would be if there were magnetic fields.
All the relevant aforementioned changes were pushed to scrambler_devel.
Test Suite update
Ruka's Bondi Hoyle problem has been successfully implemented to the test suite complete with bear2fix.data files and images. I also updated the test suite page, so that it now contains information on what our testing procedures actually do. The OrbittingParticles test will be coming soon; I just need to resolve some runtime Hypre errors. Also, Matt is working on the wiki integration which should also be finished soon.
comparing chombos
The bear2fix routine that compares chombos was altered to be more easily used with the test suite. The errors are now consistent with what is described on the Test Suite page. Basically, the relative errors are now computed as the error relative to the reference chombo as opposed to relative to some average. Due to some foresight on Matt's part (when considering wiki integration), the data in the output files is now formatted as floats with one number on each line.
Meeting Update 11/15/2011 - Eddie
I have completed updating my hydro code with a second order MUSCL-Hancock scheme. See the comparison plots on my buildcode page (scroll to the bottom). Also, I finished updating all of the content on the page which now includes short descriptions of the approximate Riemann solvers and the MUSCL-Hancock scheme.
The last few chapters of Toro deal with source terms and higher dimensions. Should I implement either of these into my code, or should I now move on to the hand-full of AstroBEAR development projects that I have on my plate?
As far as testing goes, I'm working on the OrbitingParticles problem. This one is a bit tricky because I want to produce an image of something visually insightful which is no problem with Visit. However, we want the post-processing scripts to be easily used for all tests which requires us to produce images with bear2fix. So I may have to add some code to bear2fix, but I have to think about this some more. Any ideas on what might be the easiest yet insightful thing to plot for this problem?
Some small test suite updates
I updated the tests directory on Clover with the necessary reference chombos and bear2fix.data files for the tests that we currently have implemented. So, as Matt said yesterday, we are now prepared to run weekly tests in order to compare chombos and create images.
I also slightly updated the Test Suite page. The list of tests in the table now has everything that is in the test suite, and the ones we are currently working on including who is working on them.
Meeting Update 11/08/2011 - Eddie
Test Suite Progress
The tests directory was updated and pushed into the development repository on Clover. See ehansen11042011. A makefile for alfalfa was created so that Matt's buildproblem script will work as is on that machine. The Orbiting Particles problem will also be added soon.
Approximate Riemann Solvers
I finished implementing an HLLC solver and a Roe solver into my hydro code. See my buildcode page for details. My next task in this area is to implement a Roe solver into astrobear.
Higher Order Schemes
I have just started my reading on higher order schemes. I plan on using the basic 1D MUSCL-Hancock scheme for my code. Is there a particular limited slope or slope limiter that I should use? Toro has a section on extending this scheme to higher dimensions as an unsplit finite volume scheme. I'm assuming this is the one we want for astrobear? Also, it is my understanding that astrobear already has routines to calculate slope limiters, so I won't have to code that part of the MUSCL scheme correct?
Tests directory updated
The tests directory in astrobear has been updated. It now contains 7 tests: Rayleigh-Taylor, Uniform Collapse, Field Loop Advection, a Field Loop restart test, and 3 Radiative Instability tests with different parameters.
So if you have properly set up your PATH to use Matt's buildproblem script, you can run (from the code directory):
buildproblem -t -np 4
This will compile and run all the aforementioned tests.
Roe Solver implemented
I'm on a roll today. The Roe Solver has been successfully implemented into my 1D Hydro code. As with the HLLC solver, I will not post all of the plots for the various tests. However, I did make a density plot for test 1 because there are some interesting things to note. I first ran the test without an entropy fix, and then ran it again with an entropy fix. Comparing those two plots is probably the most interesting thing you can do with the Roe solver. Again, you can see the plots on my buildcode page all the way at the bottom.
I have yet to add documentation to my buildcode page about the approximate solvers, but will do so in the near future. I do not plan on writing many details, but rather just a brief overview of the methodology behind each solver.
Now that I have an understanding of the Roe solver, I can start coding it into AstroBEAR. On second thought, this might be a bit more complicated. I will need to learn more about extending this to 3D and different equations of state.
So for now, I will start reading about higher order Godunov schemes to put into my hydro code.
HLLC solver implemented
I have successfully implemented an HLLC solver into my 1D hydro code. I decided to go with HLLC instead of HLL since it does a better job of resolving contacts. All the plots from the test runs are consistent with Toro. I posted the test 1 plot for density to show how the HLLC solver can actually sometimes do better than the exact solver. More specifically, it does a better job of resolving that entropy "glitch" inside the left rarefaction. The comparison plots are on my buildcode page at the very bottom.
I'm now working on a Roe solver, and after that is finished I will add some documentation about these two approximate solvers to my buildcode page.
Meeting Update 11/01/2011 - Eddie
- I finished my first hydro code last week. Just as a reminder, it is one dimensional, uses an exact Riemann solver, and uses a first order Godunov upwind scheme for updating.
- I ran the same tests that Toro runs, and fixed a couple of errors. All the plots for these tests can be found on my build code page.
- Also, this page is getting more and more detail on it every week. There's still quite a bit of information that I want to put on it that I've already learned. Just haven't gotten around to it yet.
- The next step is to put some approximate Riemann solvers into my code. I've already begun reading about them, and will probably have this step done by the end of the week. Should I systematically go through Toro and put all approximate solvers in my code, or are there just a couple specific ones that I should use?
Hydro code finished!
I have finished building my first 1D Hydro code. It uses an exact Riemann solver, and a first order Godunov upwind scheme for updating. At first it was not giving the correct solution, and after staring at it for a while Jonathan found a typo in one line in which the wrong gamma constant was being used. Now it seems to work great. I will run the 5 different tests in Toro, and make plots to show in next week's meeting. Also, I'm working on a separate program to compute the exact solution to directly compare it to the numerical one.
I have definitely learned a lot through this process, and my understanding of the numerical methods is much clearer now. I'm currently working on my buildcode page which documents the important aspects of a 1D Hydro code.
Meeting Update 10/26/2011 - Eddie
Following the examples in Toro, I finished writing an exact Riemann solver and a 1D code that solves the inviscid Burger's equation. It uses a 1st order Godunov upwind scheme. I'm still trying to test it and also write a code that can solve more general problems, but I think it would be good to start looking at approximate Riemann solvers now.
Meeting Update 10/18/2011 - Eddie
- put together files for restart test
- read Toro Ch.6 - The Method of Godunov for Non-linear Systems
- solved Riemann problem analytically
- looking at exact solver in astrobear, will start coding tomorrow
field loop restart test
In addition to the normal tests that we are going to run in the test suite, Martin thought it would be a good idea to have one restart test. The 2D field loop advection problem was chosen for this. I completed all the runs to make sure it worked, copied the necessary files over to the test directory in grass, and updated the wiki page: TestSuite/FieldLoop2D. This test will be run just like all the other tests in the test suite. The only difference is that this test comes with an out directory that already has a chombo, so that astrobear knows when/where to start from.
Meeting Update 10/11/2011 - Eddie
Agenda:
- Test suite preparations done ehansen10062011
- Read Pat's review article on measurements of physical conditions in jets, sent him an email…hopefully start learning specifics about his cooling table soon
- Read Toro Ch. 5 - Notions on Numerical Methods, gathering thoughts/ideas for Wednesday's code meeting
- Galactic Dynamics course project ideas? ehansen10102011b
- Group organization ideas? ehansen10102011a
Galactic Dynamics Course Project
I have an opportunity to do some research/learning that is relevant to both AstroBEAR and a course that I am taking. Eric Blackman is teaching Galactic Dynamics this semester and I need to do a presentation and write a paper. The topic is supposed to be "controversial" and probably somehow related to galactic dynamics. As an example, Blackman mentioned the Dark Matter vs. MOND debate. I'd like to do something that will help me understand some physical process that is either already in AstroBEAR or maybe something that might go into AstroBEAR in the distant future. My favorite idea so far is to do something on AGN jets. Does anyone have any suggestions for a more specific AGN jet topic? Or any different paper ideas on a different topic related to galactic dynamics?
Group Organization
I think it's important that we clarify what google groups we need and what teams we should have. We also need to discuss who should lead each team, how the team should function, etc. In particular I'm thinking of the Development, Debug, Testing, Documentation, and any other teams. We'll talk more about this in tomorrow's group meeting, but this blog post is intended to generate a discussion so please comment!
Test pages updated
I've finished updating some test pages and getting the bear2fix.data files ready. The following test pages have been updated, complete with images from the reference chombo: TestSuite/RayleighTaylor2D,TestSuite/UniformCollapse,TestSuite/FieldLoop2D. All bear2fix.data files necessary to reproduce images from a new chombo, and those necessary for chombo comparison have been copied to Martin's test directory in grass.
The only slightly tricky one was UniformCollapse since it is a 3D problem. I had to make a very minor change in bear2fix that makes a density plot in x and y by averaging over z. So to reproduce the image for that problem, you'll need the current version of bear2fix which I already pushed to the devel repo.
Meeting Update 10/03/2011 - Eddie
Activity last week:
- Read Toro Ch. 4 - The Riemann Problem for the Euler Equations
- Read Dopita Ch. 3 - Collisional Excitations
- Created Development Procedure page
- Updated test suite pages
- Gathered .data files for RT and UniformCollapse tests
Development Procedure
I just finished altering astrobear and adding my module files, Since the procedure I went through is pretty standard, I thought it would be useful for others (specifically new users) to know the correct way to do this. Here is the new wiki page: Development Procedure. It is linked under the Tutorials page and also the Code page. Please let me know if there are any errors or if you think more detail should be added.
Meeting Update 09/27/2011 - Eddie
Mostly just been reading this past week. From Toro I have finished Ch. 2 (Notions on Hyperbolic Partial Differential Equations) and Ch. 3 (Some Properties of the Euler Equations). I am currently on Ch. 4 (The Riemmann Problem for the Euler Equations). I see that Toro will start getting into some Fortran examples in this chapter, so pretty soon we should decide on what kind of 1D Hydro code I should write.
Also, I made some changes to astrobear when writing my RT module to make uniform gravity work. I plan on committing these changes to the repo, but first need to run tests on bluehive. However, mercurial is currently not working on bluehive, so I am waiting for CRC to fix that.
Meeting Update 09/20/2011 - Eddie
Progress made last week:
- New RT page that summarizes everything I did - ehansen09142011
- More detail on AMR control - ehansen09162011
With RT finished, I have begun reading more. So far I have read Dopita Chs.1 (What is the Diffuse Universe?) and 2 (Line Emission Processes), Toro Ch.1 (The Equations of Fluid Dynamics), and Hartigan's review article on the Measurement of Physical Conditions in Stellar Jets
Playing with AMR
For much of the latter part of this week, I have been using my RT module to learn how to control mesh refinement. I've added a few things to the ControllingRefinement page that should help new users. There's a new subsection on how changing qTolerance and DesiredFillRatios affects the mesh with some snapshot examples. There's also more in the ProblemSetErrFlag section now with a sample code. Since a lot of this information is useful on several wiki pages, I've placed links where appropriate.
RT finished, see wiki page
I've conquered 2D Hydro RT. All necessary simulations and calculations are complete. See my new wikipage for details on the process that I went through u/ehansen/RT. Is this an appropriate/good spot for a page like this?
Martin, let me know which files you want for the test suite and I'll get them organized/ready for you.
Tracking the RT Interface
At first, I was skeptical that my bear2fix routine was really tracking the interface of the RT instability. I plotted the bear2fix values over the visit density plot to check this. It appears that it tracks the interface fairly well. Certainly good enough to get a decent value for the growth rate. Interface Movie
Meeting Update 09/13/2011
I have been working on the RT growth rate calculation. I am currently trying to track the interface between the heavy and light fluids by tracking the position of the density gradient. The simulation growth rate is calculated through a line of best fit method. This method is very sensitive to which data points you use. There is one specific range of the data where the simulation growth rate is calculated to be less than 0.25% relative error with respect to the analytic growth rate. This region is from t=0.5 to t=1.5 (i.e. a short time after the perturbation and before the the instability becomes nonlinear).
Also, my simulation has some improper nesting which may be responsible for the asymmetry that we saw last week. This probably doesn't affect the growth rate calculation, but nonetheless I'm playing around with the qtolerance and fillratios to get rid of that.
Meeting Update 09/06/2011
Uniform gravity is now working in my RT module. Here is a movie of the latest simulation. It is a 50 x 150 grid with 2 levels of AMR.
RT Movie (Test Suite Simulation)
The only thing left to do is run bear2fix to calculate the growth rate, and check it with the analytic value. If the values reasonably agree, this will be all set to go into the test suite.
my schedule this semester
Since I won't be at the meeting tomorrow on account of studying for the prelim, here's what my schedule looks like this semester…
Galactic Dynamics - TR 2:00 - 3:15 Intro to Plasma - TR 3:25 - 4:40 Hydro Stability - MW 4:50 - 6:05 (If I do take this course it will only be an audit, so consider this time slot as open if needed)
Meeting Update 08/15/2011
Still studying for the prelim like a madman. In regards to research, I need to figure out how to turn on uniform gravity for my RT module. Shule pointed me in the right direction last week, but I'm not sure if it's working yet…update on this to come later this week.
Meeting Update 08/08/2011
I've really been focusing on studying for the prelim for the past couple of weeks, so I have not done much in terms of research. Perhaps sometime this week Shule can show me how he calculated growth rates for the RT Instability. Also, once Martin finishes the scripts for the testing suite, I will be able to do more towards that project. However, the prelim will remain my top priority for the month of August.
RT Update
I ran a higher resolution simulation with my module, all parameters the same as before. The grid is now 50 x 150 (half of the first run), but AMR is turned on with 3 levels. So the effective resolution is 4x the first run? RTmovie (Higher Resolution)
Module Tutorial Updated
I've finished updating the Wiki tutorial page on writing a problem module. The biggest changes were in the initializing a grid routine. I put in more info on the q-array, which USE statements are required, a heavily-commented sample module, and a brief bulleted list of what I felt was most important.
The q-array might have more details that I missed. Also, I did not add anything about the aux-array since I don't know much about it. For the same reason, I didn't touch the Additional Physics sections.
RT Instability appears to be working
I have a seemingly successful simulation of the RT instability. The simulation is in 2D with a fixed grid of 100 x 300. The next step is to crank up the resolution, turn AMR on and then calculate the growth rate to check if it is indeed working accurately. RTmovie
Field Loop Advection Practice
So for the past couple weeks the new users' group has been playing around with the Field Loop Advection problem. I've looked at a lot of the data files, just trying to understand what each one is for. In my most recent run of the problem I changed it to 2D and ran for 20 frames instead of the default 10. We've also been doing some intro to visit. For this run, I created this simple movie that shows the magnetic field vectors: Bvec_movie
I realize this entry is not too exciting, but it's also helpful just for getting used to editing the Wiki, posting blogs, linking different pages, etc…