I did some simulations of the circumnuclear disk with a larger inner cavity. My previous simulations, which are shown in the first figure below, have an inner cavity with a radius of 1 pc. The new simulations have an inner cavity of 2 pc and are shown in the second figure.

The plots show radial profiles for the surface density for Model 1 (without magnetic field but with outflow, dotted lines), Model 2 (with magnetic fields and with outflow, dashed lines) and Model 3 (without magnetic fields and without central outflow, solid lines). The vertical line marks the location of the disk's initial inner rim. The time is given in units of the orbital time-scale at 1 pc, which is 4.5e4 yr. The profiles have been calculated by using radial bins of the surface density with a bin size of 0.02 pc.

In the simulations with a small cavity we found that for Model 1 the disk was moving rapidly inwards, and that for Model 2 the inflow was to some extent suppressed by the magnetic field.

With a larger cavity the situation seems different, now for Model 2 the disk material moves inwards much faster than for Model 1, only after a very long time (note that the last subplot of the second figure corresponds to 22 orbital timescales and not 13 orbital timescales as in the first figure) the surface density of Model 2 is again lower than the surface density of Model 1.

According to these new results, is our previous statement that the magnetic field suppresses the inflow of material towards the black hole still valid?

To measure the magnitude of numerical viscosity in disk simulations with AstroBEAR I have performed calculations of a spreading ring.

A short reminder on the spreading ring problem:

If we assume an accretion disk that is rotationally symmetric and geometrically thin and that its angular frequency does not change with time we can derive the following equation that describes the time dependent evolution of such an accretion disk:

\frac{\partial \Sigma}{\partial t} + \frac{1}{r} \frac{\partial}{\partial r} \left[ \frac{\frac{\partial}{\partial r} \left( \nu \Sigma r^3 \frac{\partial \omega}{\partial r} \right)}{\frac{\partial}{\partial r} \left( r^2 \omega\right)} \right] = 0

Here \Sigma is the accretion disk's surface density, \omega its angular frequency and \nu the viscosity of the gas. For a Keplerian gravitational potential, constant viscosity and an initial condition in the form of a delta peak at position r_0

\Sigma \sim \delta (r-r_0)

this equation has the following analytical solution (eq. 1):

\Sigma (r,t) \sim \frac{1}{12 \pi r_0^2 x^{1/4} \tau} \text{exp} \left[ - \frac{(1+x^2)}{\tau} \right] I_{1/4} (\frac{2 x}{12 \tau})

with

x=r/r_0

\tau=\frac{\nu t}{r_0^2} (eq. 2)

and the modified Bessel function of the first kind I_{1/4}.

As a delta peak is numerically difficult to handle I use this analytical solution as initial condition, with r_0 = 2\,\text{pc} and initial values for \tau of \tau_{\text{i}}=0.01 and \tau_{\text{i}}=0.001, respectively.

The following movie shows the general behavior of such a ring. The numerical and physical parameters are the same than those of my simulations of the CND unless stated otherwise, e.g., the resolution at the location of the ring is about 0.04 pc. However, I switched of cooling and magnetic fields for the spreading ring calculations.

First animation: surface density for $\tau_{\text{i}} = 0.001$​

The following two figures show the radial surface density profiles, the first one for \tau_{\text{i}}=0.001 and the second one for \tau_{\text{i}}=0.01. For each snapshot of the simulation I fitted a curve according to eq. 1 to the surface density profiles, which are also shown in these figures.

This fit allows to determine the corresponding value of \tau. Thus we have \tau as a function of time, determining the slope of this function gives the viscosity \nu (see eq. 2).

I define the parameter

\alpha_{\text{num}}=\frac{\nu}{\nu_{\alpha}}

where \nu_{\alpha} = h \cdot c_{\text{s}} is the maximal alpha viscosity, with h=0.2 \text{pc} and c_{\text{s}}=1300 \frac{\text{m}}{\text{s}} I use typical values of the CND.

For the \tau_{\text{i}}=0.001 simulation I get a value of

\alpha_{\text{num}}=(3.93 \pm 0.13) \cdot 10^{-2}

and for the \tau_{\text{i}}=0.01 simulation a value of

\alpha_{\text{num}}=(2.86 \pm 0.08) \cdot 10^{-2}

Considering that in the literature a value of \alpha=0.1 is often used to account for viscous disk evolution, this result is not extraordinarily good, but I think it shows that our simulations are not dominated by numerical viscosity.

I furthermore did a spreading ring calculation with central outflow, the following movie shows that the material of the ring is slowly blown away by the wind, contrary to the simulations of the CND where the material of the disk was moving towards the central black hole due to its angular momentum loss.

Second animation: surface density for $\tau = 0.001$, with central outflow​

The following figure shows the corresponding radial surface density profiles. Because the development of the ring is completely different from the analytical solution of the spreading ring, I do not show any fits to eq. 1 here.

Some feedback on Erica's paper

Shape

I can visualize the data and have it imported into the 3D Module. However I have yet to be able to get it to be viewable in the Render Module. I can rotate the simulation around and such. Martin and I have used the same set of data. However you'll notice that it looks a bit different than mine. I think it depends on how many lines you want to import. However there seems to be 'two' ways to import it. I want to meet with Martin once more, and neaten this up.

Marissa's Shape

Martin's Shape

• Would anyone be interested in me giving a tutorial on how to use SHAPE?
• Once I figure this out I'll write a tutorial here on the wiki.
• Need to get a hold of Bruce's data and regrid it to ascii.
• Any chance SHAPE can take HDF5 files yet?

Wire Simulations

Made an executable that makes CDMs. Here is a preliminary taste from the first submission of the jobs. It goes to frame 22.

CDM down-x	CDM down-y	CDM down-z
​GIF	​GIF	​GIF

• Visualizing the actual hdf5 files on Clover and BH2 using the contour operator kills visit. Not sure what to do. Clover's visit didn't even like visualizing the bov files, so I did these on Bamboo, and even Bamboo had to wait between each image it spit out.
• Updated the ​Run Statistics page. Currently have 114,520 sus left.
• hydro is at frame 54, and mhd is at frame 36.
• Can no longer write to BH2. Should I delete the 2D simulations from before?

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.

g1.4_d5.5_highcontrast.gif​

With cooling, clump separation = 3 rclump:

cooling_d3_rho.gif​

New gamma = 1.2 runs, with d = 3, 4, 5:

g1.2_d3d4d5_rho.gif​

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

movie​

movie​

movie​

Here is a snapshot of frame 51 of 200 of the hydro run and frame 34/200 of the MHD run.

Monopole expansion works in 2.5D

Here is a comparison of phi for the same 2.5D and 3D simulations of a clump at the origin:

The boundary conditions in the 3D case are reflect on the x,y,&z axes, and in the 2.5D are reflect-cylindrical on the r and z axis. For the other box sides, the only possible BC is multipole, which has monopole, dipole, and quadrapole terms possible in 2.5D now (but only monopole in 2D). Both runs use the same solver.

Here is a line out of phi:

You can see they are very close. The spreadsheet curve shows their differences, the highest being inside of the clump, closest near the boundary, and then approaching zero just outside of the clump.

For reference, here is a comparison of the 2D simulation with the 2.5D:

Here, the 2D simulation setup is exactly the same as the 2.5d case, except the boundary condition is plain reflecting, and the solver is structPCG instead of structGMRes (which gave a segfault..). We can see that the 2D self gravity did a worse job of reproducing the full 3D spherical behavior, as expected. Most notably here is the puffier appearance of phi. This is because of the lack of 1/r terms in the stencil in xy space compared to rz space. In rz space, phi is more sharply peaked.

Next steps

This simulation (clump at origin) only produces a monopole term in the multipole expansion, and so we would like to also check the behavior of higher orders terms using different mass distributions.

A recent paper used PLUTO to carry forward full MHD orbital dynamics of planetary winds interacting with a stellar wind.

This paper is exactly where we want to be going. So for Jonathan, Baowei and I reviewing this in detail is a must. The authors are Matsakos, Uribe & Königl

​http://adsabs.harvard.edu/abs/2015arXiv150303551M

Here is my literature list, Adam (still have to fill in the CF section, have read a lot of papers there, and need to grab them):

​https://astrobear.pas.rochester.edu/trac/wiki/u/EricasLibrary

I will have some results up to share on my 2.5D gravity ticket for tomorrow's meeting (last week was working on that)

Here is a new git page I made: ​https://astrobear.pas.rochester.edu/trac/wiki/u/erica/GitRepos

Conferences:

IAUS 315: FROM INTERSTELLAR CLOUDS TO STAR-FORMING GALAXIES: UNIVERSAL PROCESSES? ​August 2015, Honolulu

Conditions and Impact of Star Formation ​Sept 2015, Switzerland

Dynamics of the Interstellar Medium and Star Formation ​Sept 2015, Heidelberg

Feedback in the Magellanic Clouds ​Oct., Maryland

Astrobiology and Planetary Atmospheres ​Sept 2015, Chile

Astro hack week ​Oct 2015, NYC

Star Formation 2016 ​August, 2016, Exeter

• 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:

movie​

• 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 (10⁶ instead of 10⁴) 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?

• Planetary Wind
  1. Refinement outside the outflow. Resolution (16 zone per radius) large window with size of 200. ​movie up to 29 frames; ​movie with rho_min fixed

• PN results from Bruce.

;

• XSEDE allocation
  1. Gordon & Comet almost gone. Stampede ~250,000 or 11.3% left. Details in wiki:ProjectRuns

The motivation is to set up a central installation for the code so a user won't need to re-compile the code every time he switches the problem module.. — Feel free to put your ideas here…

I. Binary Folder

1. Compile all problem modules and generate an executable file for each module.
2. Create a bin folder which contains all executable files (can be links)
3. Sample Data files folder which contains all the data files for each module

II. AstroBEAR Library

1. compile the compute engine as a library
2. each module problem we have now or the user-developed problem module will link the library and make

Doing reading for next project on outflows.

Did CIRC poster session.

Met with Alyssa Goodman.

Worked on my ticket for 2.5 self gravity

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:

MachStems.

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

Wire Simulations

• Thanks to Baowei we got the code to make on Stampede. Submitting jobs there now for the run directories Jonathan forwarded to me.
• Made a ​stampede statistics page so we can pay attention to our resources.
• Writing up a study on the cpu usage needed for these, can find it there soon. Going to also do a write up on what is different about the production runs on the ​studies page too.

Shape

• Met with Martin on Friday. Showed him what I've been doing. I forwarded him the data I've been using. He downloaded Shape and is going to play with it himself. We plan to meet again tomorrow.
• Martin says that in the past the software usually gets rid of the mesh when your data has been properly imported into Shape. This means there might be something weird with how I am importing it, or visualizing it in the 3d module.
• Scroll to the bottom of the .dat file to see how many lines there are. One can set this to the lines per entry..
• Find the max and minimum position values for the whole data set, take the average of these and use it for the "center."

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.

movie​

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.

I ran the rotating variant of the TrueLove Collapse problem. Here is an animation with density contours and the particles along with their spin vectors.

Here is the angular momentum of the gas, the particles, and the internal spin of the particles in the X, Y, and Z directions (top to bottom)

And here is a closeup of Jz of the gas, the gas+particle, and the gas+particle+spin

I'm a little surprised by the ratio of spin to angular momentum, but I guess in a binary system, there is a lot more angular momentum in the orbit of the stars, then in the rotation of the stars themselves. It is also surprising that the Z component of angular momentum of the particles becomes positive at times - but perhaps that has more to do with the choice of the center of the origin…

In order to understand the context of the discussion I am going to lay out here, please view the following video made by the makers of SHAPE v. 5 ​Shape Version 5 - External Data Visualization. Back in March I met with Marcus at RIT, and we worked on this together. What I will be detailing here are comparable results to what we saw. However we were really too excited that we saw something, and thus we weren't super skeptical of the results.

Attached to this blog post I have included two .dat files, or ascii files, from Martin's directory — in case anyone wants to download Shape, follow the video, and try to find it out. They are titled SHAPE00000.dat and SHAPE00022.dat.

Video Summary

1. Open Shape. Go to the 3D Module. Insert a "mesh" by clicking either sphere, torus, etc. The point of this is to "hold" the external data on some shape.

2. Rename your Shape, and choose what file format you want. We are using ascii.

3. Input properties dialog.

Questions & Concerns:

1. What is dNum ? The narrator in the video says that it takes a "slice" of the data. Does this mean that you only consider a certain number of rows in your data? Why 256, why 1? Especially when you are only inputting a single frame of your data. There is also this parameter n in the input properties dialog. Wouldn't these be the same thing based on what the narrator has said?

2. The narrator says that n is the "lines per entry" in the data. However wouldn't we just have 1 row based on how people typically use file I/O? When you expand the input parameters dialog box, the data you can see will adjust with the box.

3. For any data set that I obtain I will need to know the format string, or the label for all of the columns in the data. I didn't ask Martin what his labels are, however Marcus and myself made an educated guess about what they could be.

4. The narrator says that one needs to put a dummy variable at the end of their string. Why do we need a dummy variable if we know all of the columns?

5. The Center… how do you know how your data is centered? Then there are extents. Marcus and I sort of played around with this. You can also adjust the radius of the mesh object from Step 1 above to see what makes sense for how your data looks.

4. Delete the modifiers already there in the 3D Module (Density, Temperature).

5. Add the emissions species by accessing the Physics Module.

Questions & Concerns:

1. In doing this I have always followed the instructions given by the narrator in the tutorial. I don't think that this is part of the problem when trying to visualize external data in Shape. I say this because we are actually able to "see" something, which I'll illustrate at the end of the blog post.

2. When do you know what species to use? One will probably need to know if their simulation requires a particular emission.

6. Apply the species modifer.

7. Now go to the Render module.

This is what the data looks like originally if you use comparable parameters from the video.

However if you change the width of the Image, it starts to look like a sphere…

And if you change the center (in the input parameters dialog) you simply just move around the mesh.

attached

A quick update. Expect more later.

• Colliding Flows: Edit 3 plots for their final version. Redid the inverse beta maps and b vs. n plots, however I will have to quickly revisit the B vs. n plots.

• Turbulence Simulations: Jonathan sent me a new problem module that utilizes a "fractal" mesh. I am going to begin to set up production runs and make the code.

-Jonathan tells me to keep an eye on cpu time for the MHD run in particular in case we need to scale down resolution. I think I am going to make a page, comparable to the one Erica made for the Colliding Flows, that tracks the statistics for each restart.

-I am assuming I will have to start running these on and XSEDE machine now. Going to test Jonathan's new problem module on BH2 first. Then work with Baowei to get code made on Stampede.

-I plan to update the page I have made of the 2D runs after I finish setting this up.

• Shape: Adam e-mailed Martin. Before I meet with Martin I am going to write a blogpost detailing my questions. Bruce also let me know he has generated some models of stellar outflows on alfalfa. Once I figure out things with Martin I'll have a go at that. Will have to regrid frames we prefer so they are accessible to Shape.

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.

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.

movie​

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.

movie​