3D Isotropic Turbulence

Here is a test of the problem module, 'Isotropic Turbulence' in 3D. The 2D version of this code ran fine (see my previous blog post here: erica02232015).

In 3D there might be a bug? Instead of isotropic sources, I am seeing columns of turbulent sources. A slice through the box in the xy plane, would resemble the 2D simulation.

Contour of rho, momentum, and pseudocolor of rho:

movie​

movie​

movie​

• OutflowWind
  1. the blowout found in 2D run blog:bliu03232015 was caused by too-short outflow_duration. Results after fixing that:

​movie

2. ​gamma=5/3; ​gamma=5/3,zoomed

• Worked with users: #438

• Working on 3D anisotropic conductivity solver

Colliding Flows

• High Resolution (~800 ppi for color, ~300 ppi for b&w) Figures made: CDMs (0, 15, 30, 60), InvBeta (0, 15, 30, 60), B vs. n (0, 15, 30, 60), Spectra (0, 15, 30, 60), Mass History (0), Model of Flows
• HR Figures to be made: Mass History (15, 30, 60), Histograms for sinks (We are thinking about these three figures I made: ​1, ​2, ​3)
• Querying for the mass history is a pain in the ass in VisIt and could take over a week. I am just going to make a post processing executable and run quick jobs on the simulations to pull the information of the rho in the box and simulation out.
• In the past few weeks I have learned a significant amount about what goes into preparing figures for publication. Not a trivial task. Visit exports pseudo-vectorized images. Exporting as an EPS is an option in Visit, but when you zoom in on the printed figure there is no scalability. So I saved visualizations done in Visit as high resolution PNGs and converted them to PDF (ImageMagick rasterizes, but that won't matter if we save as a PNG anyway). Some other softwares can export actual vectorized images. These you don't want to touch with ImageMagick (use ghostscript instead). So we have a bunch of really heavy figures. I wrote a script on how to compress them, see ​here. Erica wants me to document all of the stuff I learned.

Shape

Meeting Marcus Freeman, an RIT graduate student, on the week of the 6th to work on Shape. Will post about the outcome of that meeting afterword. Before then I should probably work on regridding the simulation so it has the necessary variables we would like?

Turbulence Simulations

​GIF
Met with Jonathan last week. He handed over the problem module and .data files. Suggested I make a 2D simulation first, so I did this over the weekend. Looks as expected: clump and wind. Hot desaturated really makes the shock stand out. There is a weak magnetic field (beta = 100) and MHD is turned on.

Space is minimal on Bamboo and Grass, so I ran this simulation using 4 processors on Clover (which is only being used about 56%). It took all weekend and has 6 more frames to go. I am going to start moving over to a more powerful machine, like BS or BH2. The configure file is having issues on BS (see Baowei's e-mail).

Here are some of the parameters in the problem.data we'll be considering for the study:

&ProblemData
spacing=.2           ! Lattice constant
thickness=.03        ! Radius of wire
beta=100             ! magnetic beta
mach=20              ! mach number
screen_x=.25         ! location of screen in x
rho_wire=1000d0       ! peak density of screen
rho_wind=1d0     ! density of wind
rho_amb=.01          ! ambient density
/

We could vary the beta and mach, or the spacing and thickness of the wires. I suppose we'll see once I get some jobs set up and the code made on BH2.

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.

In order to create a 1-D spherical solver, I write a python version of Riemann solver. The long term goal is to add cooling, thermal conduction and radiation to examine theories and provide an accurate benchmark for AstroBEAR.

Simulations are run with 500 points.

Below are results from Roe solver: (each picture consist of rho, v, P)

Below are results from HLLC solver:

I am developing TVD-MUSCL scheme now, Toro's book is wrong on this Chapter. I have succeeded in writing TVD-MUSCL 2 years ago.

Binary seperation = 5 AU

Secondary mass= 0.15 SM

Secondary accrete gas around it

Isothermal temp = 1000 K

Primary's parameters:

Primary mass = 1 SM

Base density is fixed at 1.18E-5 g/cc throughout the primary (even in equatorial outflow)

Alpha of the secondary is 1. Alpha of the primary varies as below:

The velocity at the base of the wind does vary near the equatorial plane. In order to mimic the equatorial outflow induced by the in fall of planet (more massive than Jupiter), its open angle is (1/9)*Pi. Its velocity surges from 0.04 km/s to 2.4 km/s. Then equatorial outflow last only for 0.6 yr and start at frame= 150, i.e. time= 16.5 yr

Density plot + vector velocity of mid-plane:

density.gif​

Z-Angular momentum + density contour of mid-plane:

The whole angular momentum is negative, i.e. anti-clockwise. Which means, angular momentum has been transferred from stars to gas.

Lz.gif​

This simulation is basically based on Parker's model. The disk is not obvious and the pulsation is also not clear. The pulsation is subsonic and smeared out by the high density gas around the primary.

Which is supposed to be seen if the simulation is correct? Disk or spiral shock?

Magnetic field strength

To get the magnetic field strength of the flows for the paper in gauss, I think I can do this,

Take the computational magnetic pressure and multiply by pscale to get magnetic pressure in cgs,

B(cgs)²/8pi = (B(cu)²/2)*pscale

Then multiply by 8pi and take the square-root to get B (cgs).

Now, I am pretty sure that the unit of gauss describes 'B', but wikipedia says the unit is for the flux density, which I suppose one might interpret as energy density aka magnetic pressure. If that were the case, one would not take the square root..

Considering just B, I get,

B = 1.3 microgauss

This seems reasonable, considering typical estimates for molecular clouds and ISM are ~ 10 microgauss..

Ram pressure

In my physics.data, xmu = 1.27, which seems to be the mean molecular weight in the code. Mean molecular weight is the mean mass of a particle in atomic mass units, or amu. One amu is technically 1/12 the mass of Carbon-12, and so is approximately the mass of 1 nucleon, or equivalently, one proton (m_H).

rho = n * xmu * amu

Xmu seems to be used to set the scales for velscale (which is a scaled sound speed). It seems then that this means the particles in the gas are effectively *mostly* hydrogen with a small fraction (25%?) being 2x the mass of hydrogen, ie Helium.

Pram = n*amu*xmu*v²

Am working on compiling my dev branch for 2.5 self gravity, and thinking of a test problem (the former seems to be going smoothly, the latter will discuss in meeting today). Maybe a few days out from being done with this ticket. Next up on my plate might be 1D gravity.

I also suggest making sure the code can compile on all local machines in the next fix of the code.

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.

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)

1. Outflowwind: density blowout in 2D For a 2D run I did last week, the density blows out after running some time..

​movie

nDim     = 2                            ! number of dimensions for this problem (1-3)
GmX      = 100,400,0                    ! Base grid resolution [x,y,z]
MaxLevel = 1 !5                         ! Maximum level for this simulation (0 is fixed grid)
LastStaticLevel = -1                    ! Use static AMR for levels through LastStaticLevel [-1]
GxBounds = 0d0,-200d0,0d0,100d0,200d0,0.d0      ! Problem boundaries in computational units,format:

2. Working with SUNY user and on anisotropic conductivity solver.

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.

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.

The gamma = 5/3 runs:

movie​

The gamma = 1.4 runs:

d = 4, 4.25, 4.5, 5.5 movie​

d = 3, 4, 5.5, 6.5 movie​

d = 3, 4 movie​

d = 5.5 movie​

I've been continuing to streamline the chombo io process and now have parallel write working. And it was 2x as fast using 8 processors on a 200 MB chombo file.

I took luminosity into consideration.

The primary's luminosity is large so the photons on the dust drive the dust and gas away from the primary. The secondary is assumed to be dim - although it actually maybe not given that it is a 0.5M_{\odot} WD. In all, the secondary star has full gravity on gas while the primary has no gravity on gas as soon as the dust form at 4 AU.

I also assumed that, there is a magic that the dust evolves such that the photon-dust interaction changes. The scattering and absorption of photon diminish at 12 AU away from the primary. This is only a test in order to see what will happen.

The primary is pulsating periodically with a period of 36 yrs and pulsation last for 1.1 yr. The pulsation is similar to the last post.

Temperature is isothermal at 2000 K. Separation is 7 AU. At this separation and temperature, the accretion mode is probably WRLOF or RLOF.

Total simulation bos is 320AU\times320AU\times8AU

density.mov​

It looks like a steady state disk as gas is continuously coming out from the primary.

Science Meeting

Projected streams

Here is another image showing projected streamlines. This time, we are looking at the final frame of the 0, 30, and 60 case from left to right. We can see from this how straight the interface became in the 60 case — it looks like the 0 case on the far left. (Thought: does the interface evolve faster in the 60 case than the 30 case — perhaps with enough time, the 30 case would also straighten out).

Projected streams vs sliced streams

Last time we saw this image of the projected streams,

This time, here is the same image but made of slices rather than projections (that is slices in density and streamlines down the middle of the box).

Coding meeting

Simulation Proposal

Added tracer to the outflow object:

​movie ​zoomed in movie

So we have been seeing lots of cool instabilities in the simulations of the planet wind/stellar wind interactions. So far all these sims have been run as isothermal with

\gamma = 1.01

This means we expect cooling to play and important role for both the shocked planetary wind (pw) and shocked stellar wind (sw). I wanted to check these assumptions and so I did a brief calculation.

We are taking about Hot Jupiter's here so the orbital radius is small

a = 0.025 AU

Stellar Wind

At that radius the stellar wind is still acceleration and is likely subsonic (M<1). That means no shock. If I assume its a solar type star with a corona the parameters are likely

V_{sw} = 500 km/s

T_{sw} = 10^5 K

and using

\rho_{sw} = \dot{M}_{sw}/( 4 \pi V_{sw} a^2)

gives me a number density of

n_{sw} = 10 g/cm^3

To get a cooling time scale I just use

t_{cl} = E_t/\dot{E_t}

and I approximate the cooling rate

\dot{E_t} \sim n^2(10^{-22} erg cm^3/s)

Putting in #s I get

t_{cl} \sim 10^9$ s

To see if the wind actually cools I need a dynamical time to compare it to which I chose to be the time crossing time of the shocked flow past the planet (I use Jupiter's radius for the planet and assume the interaction happens at 10 R_j)

t_{pc} = 10 R_p/V_{sw}

This gives me

t_{pc} \sim 10^7 $

Thus

t_{pc}/t_{cl} \sim 100 $

and the stellar wind will not cool.

Planetary Wind

For the planetary wind shock things look a bit better but not too much.

This is density [ [Image(movie10019.png,width=300)]]

attachment:movie1.gif​

This is contours of velocity and energy

Working on the paper, and see last week's meeting post:

​https://astrobear.pas.rochester.edu/trac/blog/author/erica

• OutflowWind
  1. latest results for 2D high-res and 3D low-res after sonic surface fixing are here: blog:bliu02282015_2
  2. got 3D high-res data, haven't analyzed yet.

• Ablative RT
  1. meet with LLE next week?

• XSEDE Allocation
  1. expires in June 30 2015
  2. Current resources: Stampede ~700,000 SUs(32%) left, Gordon 520,000 SUs (74%) left, Trestles 100,000 SUs (100%) left. Details can be found at: wiki:ProjectRuns

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.

movie​

I think there is a good story here, and this is my interpretation:

1. Bow shocks intersect.

2. Bow shocks get "puffed" up due to increase in pressure in intersection region.

3. The intersection angle crosses the minimum critical angle and a Mach stem forms.

4. The Mach stem grows until the intersection angle hits the maximum critical angle.

5. 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:

1. When the clumps are too close together, the Mach stem will grow beyond the maximum critical angle and disappear.

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

3. 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:

movie​

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.

movie​

FYI: I made my own Library (copying Erica). I wrote an outline about the 2007 Shape paper, ​here. In case anyone is interested about the point of Shape.

Since our group is trying to utilize the Shape software for our simulations, I figure I would start to document helpful components to visualizing and modeling in Shape. Here I will briefly explain how one can make a movie, or animation in Shape. Some days toward the end of the day, I am going to go through the templates on their website. With Shape however, I have two main objectives:

1. Learn how to convert our hdf5 simulation data to ASCII, so that it can be fed into Shape.

2. Learn how to visualize this external data in Shape. What limits are there and what can Shape tell us about our data compared to other visualization softwares?

So if the template I check out is useful, I'll make a blog post about it explaining the objective of the template. Eventually I'll use these to make a series of wikipages documenting Shape for our purposes. However I am yet to be a "Shape Master." Otherwise I would have done 1 & 2 already. This template is quite simple: how do you make a movie in Shape?

Results:

(Make sure you check out the .gifs! They are really pretty.)

Image 1	Image 2
​BW GIF ​Red/Blue GIF	​GIF
Image 3	Image 4
​GIF	​GIF

This data is from one of the ​Shape templates, titled Animation 1 - Rotation. It is a bipolar nebula rotating around all three x-y-z Cartesian axes. Image 1 is a 2D projection, comparable to our column density map movies. However the object is not evolving in time. Image 2 is the line profile for the velocity vs. pixel intensity as the object rotates. Image 3 is the mesh used for the object. Image 4 is its PV-Diagram.

Guide:

Image 5

This screen-shot was taken post-rendering. Click the big button with the Shape-S in the top, left-hand corner to view the data as seen here. Prior to rendering, the PV-diagrams will be black. Note that all of the data parameters are already in the template. You can adjust the colors if you want too. For instance I have also made a movie color coded to indicate Doppler-shifting, i.e. Red/Blue (see above). Other options include grey scale, rainbow (color), red/blue, gradient, and spectrum. Click the movie-film looking button in the same row as the renderer to get to the animation module.

Image 6

Image 6 is what the animation module looks like. Instead of clicking the rendering button, you'll click the animation button (which is circled and marked as the second step). First you want to adjust what format your output will be. Note the timeline at the bottom of the GUI. As Shape makes your animation, it will tell you how far it is in the process. On the right hand side of the Animation module is the variable tree. Below the variable tree, in the table, is what is referred to as the variable stack. Essentially the functionality of this part of the animation module in v. 5 seems not much different from v. 4. Y

Here is a table of screen shots of all the parameters listed in the animation module.

Image 7. General

Image 8. Variable

Image 9. Output

Image 7 (in order from top to bottom):

• Name of your file output
• Number of frames (where you want them to start and end in the simulation)
• Current frame tells you when Shape is in spitting your .png output (there is also a timeline at the bottom of the GUI that starts from first frame and ends at the number of the last frame) relative to the simulation
• Start and end times of your simulation for some given time units (years, days, hours, minutes, seconds).
• Like current frame, it tells you what current time your simulation is being animated at in your chosen units.
• Distribute & Fields: Using particles for rendering or output velocity vector information may require Shape to redistribute the particles in every frame. This is done by default. If you are not using particles, then disabling the distribution can save processing time.
• Render: Render while you animate. This is essentially required.

Image 8. (Currently investigating)

Image 9. Above in the results section of this post, you can see all of the options visualized that are listed here: 2D Maps, PV Diagrams, 3D Mesh, Hydro, Plots (Images) (Note: the numbers to the right indicate LxW size of the image, so these plots are 512x512), Plots (Ascii), Math Variables, and Time Units. You can denote the image type, and indicate where you want the images to be saved (a working directory).