So Eric and I were discussing the specific angular momenta in the gas yesterday and you can calculate the velocity of the wind at any given point (assuming the secondary does not influence it's motion) by solving for the characteristic that leaves from the surface of the primary with a velocity vector pointed at the position. It is difficult to solve exactly and some approximations must be made… If we assume that the wind velocity is much greater than the orbital velocity, and that the distance to a given position is much greater than the primary's orbital radius, then we can easily esimate the time delay (t-t₀) between the wind leaving the primary at t=t₀ and arriving at a position x at time t as t-t₀=|x|/v_wind. We can then determine the location of the primary at this time, and then assuming that the primary's radius is << |x|, we can estimate the direction of the wind as being the same as the vector going from the primary towards |x| or v(x,t) x (x-x_p(t₀)) = 0. Then to calculate the magnitude of the wind we know that it will be constrained between v_w +- v_orbit depending on the degree to which v_p(t_0) is aligned with (x-x_p(t₀))

See the following illustration

The upshot, is that the secondary tends to be in a region where the specific angular momenta from the primary is positive (which is retrograde)…

Here is the 'corrected' approximate solution plotted (here the orbit is clockwise and the z-axis if pointing out of the plane)

Click here​ for the original version. Note where |x| = 1 - there is no difference

And here is another more exact solution where the 'time delay' is not calculated based on the distance from the origin, but by the distance to the current location of the primary. The exact solution would require solving a transcendental equation to find the distance to the position of the primary at the retarted time…

and here is the calculations that went into the above approximation

And here is the exact solution calculated from the simulation

Here the yellow corresponds to positive specific angular momenta - and the system as a whole has negative angular momenta - so yellow ⇒ retrograde disk

Here is a few frames later after the secondary has begun to torque the wind…

The secondary should likely give a net torque that is prograde (so more blue then yellow) but positive specific angular momenta does not automatically form a disk since the disk has to be orbiting the system's center of mass as well as orbiting the secondary. So does the infalling material have to have a net specific angular momenta that is greater than the secondary to form a prograde disk?

And here is a zoomed out version of the approximate solution… Note the whole thing rotates with the same pattern speed - so just spin your monitor to see a movie

Also was thinking about what the flow looked like from the point of view of the secondary - I don't see any shear in the flow, but there is of course divergence. But in the co-rotating frame, the secondary sees a nearly uniform flow subject to a Coriolis force as well as a gravitational force… something like

So setting up the initial solution can be accomplished by the following:

First start with the assumption that the retarted time is the current time

t_r=t

DO

Calculate the displacement vector from the primary at the retarded time

\vec{d} = \vec{x}-\vec{X}_p(t_r)

Calculate the wind normal so that

(v_w \hat{n} + \vec{V}_p(t_r)) \times \vec{d} = 0

\hat{n} \times \vec{d} = \frac{-1}{v_w}\vec{V}_p(t_r) \times \vec{d}

It is possible to solve \hat{a} \times \vec{b} = \vec{c} for the unit vector \hat{a}

\hat{a}=\frac{\vec{b} \times \vec{c} + \sqrt{\vec{b} \cdot \vec{b} - \frac{(\vec{b} \times \vec{c}) \cdot (\vec{b} \times \vec{c})}{\vec{b} \cdot \vec{b}}} \vec{b}}{\vec{b} \cdot \vec{b}}

Calculate the wind velocity from the primary

v=|v_w \hat{n}+\vec{V}_p(t_r)|

Update the retarded time using the new distance and wind speed

t_r=t-\frac{d}{v}

END DO

The only problem occurs when there are multiple solutions for the retarded time…

This will occur once we reach distances of order v_w/\Omega

If we switch to a rotating frame that rotates counter to the orbit so the angular speed is \Omega, then

\vec{v}_r=\vec{v}-\Omega \times \vec{x}

and

\alpha=\alpha_0+\Omega t

so that

\Omega t_r+\alpha = \alpha_0 - \Omega \frac{d}{v_w}

Resources:
​http://www.rcsg.rice.edu/home/
So the rice machine seems to require a netid which is only for students or staff members at rice.
​http://www.rcsg.rice.edu/apply/
I guess we need to contact one of Pat's students or Pat himself for setting up a guest account.

clump sims:
Z Polarized Clumps with different beta:
Density:
​http://www.pas.rochester.edu/~shuleli/SingleModeMovies/betacom.gif
Field Energy:
​http://www.pas.rochester.edu/~shuleli/SingleModeMovies/fmcom.gif
Field Components:
​http://www.pas.rochester.edu/~shuleli/SingleModeMovies/bcomp.gif

X Polarized Clumps with low beta:
​http://www.pas.rochester.edu/~shuleli/SingleModeMovies/xplb.gif
Field Components:

Uniform Bx with reflective boundaries:
​http://www.pas.rochester.edu/~shuleli/SingleModeMovies/refbx.gif
Uniform Bx and By with reflective boundaries (snapshot at different time):

Mostly working on turbulence stuff and developing routines for creating random fields using FFT's as well as spectra…

Created Layout objects which allow for a distributed fixed subgrid across all of the processors. There are also routines for loading data from the AMR mesh into the layout and extracting data from the layout into the AMR mesh - as well as transferring data from one layout to another

Created PFFTPlan objects which allow for easily performing parallel FFT's of data. It creates nDim layouts where in the 1st layout there are no divisions in x, in the second layout there are no divisions in y and so on…

Creating spectra processing objects which use PFFTPlans to load various fields from the AMR mesh into a local data array - then perform a parallel FFT - and then bin the results into radial bins and output the results in a curve file. Supports scalar fields (density spectra) as well as vector fields(velocity spectra…) Vector fields output the spectra of each component as well as the spectra of the total, the solenoidal, and the diverging components.

Plans can be used to generate initial conditions much faster than summing up sines and cosines… And layouts should allow to take the data in the base grid - and remap it on to a subsection of a finer level (for using forced turbulence results as initial conditions for a clouds velocity field…)…

movie​

Surprise​

1. Tickets:
All active tickets were processed—either assign to new owners or closed. Details can be found at ​https://clover.pas.rochester.edu/trac/astrobear/blog/bliu03252012.
I'm testing to use Trac tickets as a way of managing my projects.

2. Optimization: New Grid-Generating Algorithm
I was working on comparing the new grid-generating algorithm with the old one using real AstroBEAR modules but ran into bug-fixing like the Ticket 184: ​https://clover.pas.rochester.edu/trac/astrobear/ticket/184. The primary results were promising: the time used from the new algorithm so far was comparable with the best case of the old one. Will work on pictures…

I'll post images soon.

I'll submit the magnetic tower paper to ApJ and arxiv this Friday.

Binary

• The 40AU run with the parameters discussed last Friday is running. It's now at its early phase and seems to be doing well.

I'm working on plots of mean angular momentum vs time. I've used the data from the simulation that showed a tilt of 90deg (see last week's post). I've only had time for a quick inspection. The facts that jz is negative (because particles orbit clockwise) and larger (about 300 times) than the other ang. mom. components, are expected. Trying to understand the gradients in plots and the relative values of the momenta.

AGN jet truncation by stellar winds. The 1st test of Model 2 is still running. It's a little more than half way though. The RG's wind looks square-ish because of (i) the decreasing radial resolution of this test (see ​https://clover.pas.rochester.edu/trac/astrobear/blog/agnJetTruncation. This will be fixed for production runs) (ii) the very high density contrast with the ambient medium. The stellar wind's momentum in this case is rather strong and able to truncate the very light AGN jet, despite the fact that the latter is about 700 times faster. I plan to start running Model 3, milder stellar wind, by weeks end or early next week.

PN wind. I'm working on the referee's suggestions. Have not finished yet. See progress in Figures 2, e.g. , 3, 4 and 5, as well as the new section 4.2 in ​http://www.pas.rochester.edu/~martinhe/PNpaper-mnras2.pdf. Plan to send it to Adam, Eric and Raghvendra soon.

CRL618. The jet version of the runs with a toroidal AGB and higher densities, with and without ambient rings, are running. Should finish by week end.

Quick look

Here are the results of my sim that was set up with same params as Bannerjee's:

Density: ​http://www.pas.rochester.edu/~erica/BannerjeeRunRho.png

Radial Velocity: ​http://www.pas.rochester.edu/~erica/BannerjeeRun_Vrad.png

Bannerjee paper link:

​http://adsabs.harvard.edu/abs/2004MNRAS.355..248B

Jonathan and I went through all the active tickets we have. Several tickets were closed and most of still active ones were assigned to new owners.

Here's a summary of what we have

A more detailed report can be found here ​https://clover.pas.rochester.edu/trac/astrobear/report/8

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*10¹⁶ 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.

Hypre

Sorted out several issues with bluegene and hypre - although I am still a little perplexed by hypre's fortran interface… Was able to use it successfully on bluegene - but when I tried the same on alfalfa - I ran into weird mpi errors… I think it has something to with the assumed length of integers passed through the hypre interface??? Our own interface works correctly - so this is probably not worth pursuing… Also discovered that when running with AMR on 100's of processors, that hypre's memory requirement may be quite large unless you use a version of hypre configured with "—with-no-global-partition". See Hypre Memory usage

Restarts

I modified the restart IO to be much faster on bluegene… Restart time was reduced from 3 hours to 3 minutes! See #180. I still need to check this into the repository…

Colliding Flows

These runs are going well on kure. Need to start thinking about analysis… Here's a movie​ of a clumpy flow. I have a similar setup but with a smooth flow running on kure now.

Gravo Turb

After fixing hypre issues, runs were able to proceed for about .25 free fall times. Then some combination of protections/resteps/etc… caused runs to 'lock up'… Maybe initial setup was too ambitious - or refinement criterion was not sensitive enough to avoid large errors in energy/gravitational potential at coarse fine boundaries… Looking at a similar setup in 2D on alfalfa to try and determine what may have gone wrong… Here's a movie​

Collected some code improvement ideas on to DevelopmentProjects

This week is really hectic with my stat mech midterm and other homework assignments. I am running the Bannerjee set up now though, but will likely not be able to post plots until after my exam on Friday..

Run is located at:

/alfalfadata/erica/PROJECTS/BE_stability/12FEB2012/Bannerjee_stuff

Params:

• R = 1.62 PC
• Rho_central = 3.35*10^(-21) (g/cc)
• Xi = 6.5
• Box size = (25 pc)³
• Collapse initiated by reducing thermal energy of clump in problem.f90
• Close or exact matches of Bannerjee's reported params (Rho_c, R, Xi, isothermal sound speed, Mass clump), except my Temp seems to be 2times greater than Bannerjee's — 40K instead of 20K. Have a hunch it may do with the fact that I am using a Xmu = 2.. Should this be something to rectify before comparing results?

New Algorithm for Patch Refinement

New algorithm is good at doing refinement for line-shape areas (​https://clover.pas.rochester.edu/trac/astrobear/blog/bliu03112012), but not good at ring-shape areas (​https://clover.pas.rochester.edu/trac/astrobear/blog/bliu03112012#comment-5) where the old algorithm probably works better.

Here's the result for 2D-search (checking splitting point along both x and y)

Now working on a new recursive 3D algorithm.

Recursive Inflection Algorithm for Patch Refinement

This new(er) algorithm combines the old algorithm and the splitting-cost-check algorithm, and so the good parts of both algorithm. The following are results for some ErrFlag patterns:

Clumps

For uniform field clumps, the data are on I've discussed with Ruka on the 2D images and line plots we want to have in the paper. He now has the data and instructions to make the images.
Here's his trial animation:

​http://www.pas.rochester.edu/~smurugan/SYMP.gif

Done some "single mode" magnetic clump runs. Here I define the "polarization" as the direction at which the poloidal field is pointing. We need to discuss two distinct polarizations: aligned with the shock and perpendicular to the shock. Since the poloidal and toroidal components can have different ratios inside the clump, we actually have one more free parameter. In the paper the shock is directed on X direction. So X polarization is the aligned case, Y or Z polarization is the perpendicular case.

Here I show the animation when the clump is polarized on Z direction, poloidal vs toroidal ratio is 1, volume averaged beta = 3.15:

​http://www.pas.rochester.edu/~shuleli/SingleMode/zpol.gif

The poloidal vs toroidal field strength is shown in the following animation:

​http://www.pas.rochester.edu/~shuleli/SingleMode/ztvp.gif

Next is the X polarized clump:

​http://www.pas.rochester.edu/~shuleli/SingleMode/xpol.gif

Comparing the density evolution, the Z polarized clump is delayed:

The higher resolution runs are about to start with volume averaged magnetic field strength equivalent to beta = 0.5 and 8, corresponding to the uniform field case.

BE Dynamo Generation

Now the subcycled runs on 1 processor. Here I show the animation of viscous MHD KH instability under when the Braginskii viscosity cycles twice per hydro step:

​http://www.pas.rochester.edu/~shuleli/MKH/MKH22.gif

Binary

See ​https://clover.pas.rochester.edu/trac/astrobear/blog/binaryBondi for images and movies.

• In general, for times <~ 1orbit, the 40 AU runs form strongly tilted disks with angular momenta vectors which are
  • contained very close to the orbital plane
  • point at the winds incoming direction; ~10^o from the primary's position.
• This tilt seem independent of:
  • the soft radius
  • whether the gas initial condition is a light static ambient medium or the AGB wind
  • q, for I tried for q=m₁/m₂=1.5,2
  • the secondary's spin
• The formation of this structure is ~4 times faster for r_soft=0 than for r_soft=2dx
• I see a less pronounced tilt in a 11au run (see ​https://clover.pas.rochester.edu/trac/astrobear/blog/binaryBondi 20 mar '12)
• The tilt doesn't happen in runs with:
  • a<~5au separation (as in M&M)
  • more resolution (because the 11-40au grids are larger) to see the Hill radii. I've not used more resolution yet, for these sims running times are rather large and wanted to see if I could get disks in sensible walltimes.
  • the 2 conditions above + v_wind=20km/s
• I do not fully understand why such a dramatic tilt in the 40au case, yet the wind velocity, the shear and the wind density near the secondary's position, are all higher in the 11au sims than in the 40au ones. Also, the No. of cells per r_Hill scales down with the size of the grid.
• As reported by email, one of my 40au test produced a disk on the orbital plane, but for:
  • an initial AGB wind with a density perturbation (which are expected in AGB stars; i.e. early pulsations)
  • boundary inflow, which was not intended.

I found the tilt again once I removed both the boundary inflow and the wind density perturbation. I'd like to try the perturbation again (without the inflow).

Next steps:

• 40 au runs with higher resolutions
• Longer runs
• these runs should take about 1.5-2 weeks to run, even with the sandwich grids I'm using which, BTW, are enough to resolve Hill radii. I've not seen significantly faster runs in bluegene so far. I suspect, however, based on previous tests, that the bluegene sims will stop due to the error reported in ticket 121.

CRL618, several new plots and movies at ​https://clover.pas.rochester.edu/trac/astrobear/blog/crl618Figures#comment-1

AGN jet truncation. Model 2 is running now, see ​https://clover.pas.rochester.edu/trac/astrobear/blog/agnJetTruncation

Magnetic tower. Waiting for comments from Eric, Pat, Sergey and Jerry. Will submit to ApJ next friday.

PN winds. I'm working on the referee's suggestions, see Figures 2,3 and 4 in ​http://www.pas.rochester.edu/~martinhe/PNpaper-mnras2.pdf

Color figures for the printed version?

Clumps:
The 4 uniform field runs are finished (fixed grid). I will meet with Ruka next Monday to discuss the analysis to do on the data.
The single mode runs will start shortly.

RT stuff:
Ported the implicit thermal diffusion code over to the new version of the code. I have been doing RT tests using Riccardo's data, which generates the same outputs as we've seen before. (spike growing at the interface)
I'll either use my own RT profile or invent some new quantitative tests to see whether the algorithm or the data is the reason in the next few days.

Subcyclings:
I have been working on the subcycling part of the code with isotropic thermal diffusion. The code is currently running fine but does not looking correct. I think we can get it done faster than we thought it would be, except for the persistent send/recv part. It will also be helpful for the RT stuff if we can compare the subcycled thermal diffusion vs the implicit thermal diffusion.

Wiki Page for Ticket Assignment Procedure

​https://clover.pas.rochester.edu/trac/astrobear/wiki/TicketAssignmentPage

Current Active Tickets

Testing Results for new patch refinement algorithm

​https://clover.pas.rochester.edu/trac/astrobear/blog/bliu03112012

Binary.

I'm continuing working on the 40AU case. High res sims have produced a "disk-like" structure, the angular momentum vector of which seems to point towards the primary. It doesn't make sense to me. I'm investigating this but suspect it may be related to the fact that the secondary pulls before the it is inside the wind (I used this setup for it worked well for our previous runs with Bondi accretion and smaller separations). I want to restart the sim (i) once the secondary is inside the wind (ii) with a secondary which starts pulling gas after restart. The restart procedure has changed in the latest code revisions, so I'm updating my module. IMAGE: 2d density maps showing a part of the orbital plane close to the secondary. Left: rsoft=0; Right: rsoft=dx. Times are different in each panel.

CRL618.

The tracer field that represents the DM cooling seems to be working well. This sim ran 40% and is scheduled to restart in bgenen soon.

As discussed last time, I increased the densities of both the clump and the ambient medium by a factor of 10 in order to get stronger cooling and thus thinner lobes. In addition, the AGB wind is toroidal. Here's and image from the low res test, the high res one should run soon.

PN paper. I've heard back from MNRAS: "…Minor revision of your manuscript is requested before it is reconsidered for publication."

Magnetic tower. Waiting for comments from Adam, Eric and Pat to submit. I had a discussion with Manoj and he's very interested in out magnetic towers to model some of their recent Hershel observations. He scheduled a talk for me and will ask Adam to give one, in a meeting that he, Dan and other people will have in Rochester in about 3 weeks.

AGN jet truncation. As I said last time, the RG wind looks strange for it is subsonic. New runs:

1. smaller grid and jet radius by a factor of 2 + more refinement levels,
2. same as 1 but for a jet kinetic luminosity of 10⁴⁵ers/s (instead of 4x10⁴²),
3. same as 1 and 2 but for a magnetized jet. Should have an effect on the RG wind.

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.

Here are the Bannerjee plots for my sub-critical sphere (xi=4) that collapsed when the ambient medium was in pressure equilibrium with the edge of the sphere and at the same density. Entire box is isothermal.

Here is density: ​http://www.pas.rochester.edu/~erica/density.png

Here is radial velocity: ​http://www.pas.rochester.edu/~erica/radialVel.png

The scaled times from the Bannerjee paper run from ~3.7-7, so my time samples are slightly different. The radial velocity trend looks nice. As for the density, it seems that the ambient medium falls inward and piles up at the edge of the sphere at earlier times. By the last time curve, the curve has the right shape as the curves in Bannerjee's paper. One big difference I notice is that with my run the ambient seems to have set up some sort of equilibrium with the sphere. In contrast, in Bannerjee's sim the ambient stays uniformly flat throughout the collapse. I see that in that run, the ambient is 100 times lighter than that sphere's outer edge and the box's boundaries are very far from the sphere. There are other differences between the simulations as well, such as there the sphere was a critical sphere with a 10% density enhancement to initialize the collapse. Also, in my simul, a sink particle formed between t/to = 3.84 and 5.12. Is it no longer correct to compare sims after sink formation?

What next steps should be taken?

The testing results of the old/current algorithm for refinement patches which includes a bug can be found here:
​https://clover.pas.rochester.edu/trac/astrobear/blog/bliu03052012

The new algorithm calculates the splitting cost at each position (along one direction for the time being) and finds the optimal (with lowest cost) position. The following shows the testing results of the new algorithm. It clearly shows that the new algorithm fixes the bug in the old one.

1. 16X16 Diagonal patch

2. 32X32 Diagonal patch

3. Random patch 1 =

4. Random patch 2

5. Random patch 3

6. Random patch 4

The code dies with the following message Out of memory trying to allocate 412094160 bytes

• nBoxes is the number of grids passed into hypre for that level's solve
• nCells it the number of cells inside of those grids
• x_min is the lowest index in the x-direction of all the boxes passed into hypre
• x_max is the highest index in the x-direction…
• y_min is the lowest index in the y-direction etc…
• MinBoxSize is the number of cells in the smallest box
• MaxBoxSize is the number of cells in the largest box
• Table₁ , Table₂ , & Table₃ are the sizes used to calculate the index table
• Table size is the product of those sizes.

So hypre was working as indicated assuming that it was compiled without the '-with-no-global-partition' option. This compiler option is strongly recommended on BlueGene for amr problems using hypre.

Here is a graphical example of hypre's global partition index table… This would results in a 9x6 table with an entry for every black dot. The horizontal lines are along the top and bottom of every grid, and the vertical lines are along the left and right edge.

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.

Stability of Bok Globules, theory and numerics:

​http://adsabs.harvard.edu/abs/2007PhDT........19R

Isothermal spheres embedded in homogenous medium: ​http://adsabs.harvard.edu/abs/1986A%26A...165....1U ​http://adsabs.harvard.edu/abs/1980A%26A....83...95S ​http://adsabs.harvard.edu/abs/2011MNRAS.414.2728Y ​http://adsabs.harvard.edu/abs/1968MNRAS.140..109Y ​http://adsabs.harvard.edu/abs/1970MNRAS.151...81H ​http://adsabs.harvard.edu/abs/2000PASJ...52..217H ​http://adsabs.harvard.edu/abs/2004astro.ph..2021K ​http://adsabs.harvard.edu/abs/1996ApJ...469..194M ​http://adsabs.harvard.edu/abs/1979A%26AS...38..295C ​http://adsabs.harvard.edu/abs/1980MNRAS.190..497K ​http://adsabs.harvard.edu/abs/2000prpl.conf...59A

Virial Thm analysis:

​http://adsabs.harvard.edu/abs/2007IAUS..237..410D

Numerics: ​http://adsabs.harvard.edu/abs/1984Ap%26SS.100..447J ​http://adsabs.harvard.edu/abs/1990ASSL..162..333B

Stability in rotating isothermal spheres: ​http://adsabs.harvard.edu/abs/1989ApJ...341..220T ​http://adsabs.harvard.edu/abs/1988PThPS..96...50K ​http://adsabs.harvard.edu/abs/1984ApJ...286..529T ​http://adsabs.harvard.edu/abs/2007arXiv0706.3205C ​http://adsabs.harvard.edu/abs/1985Ap%26SS.113..125C ​http://adsabs.harvard.edu/abs/1988PThPS..96...50K

MHD stability analysis: ​http://adsabs.harvard.edu/abs/1982PASAu...4..371P ​http://adsabs.harvard.edu/abs/1997RoAJ....7..143B ​http://adsabs.harvard.edu/abs/1976ApJ...207..141M ​http://adsabs.harvard.edu/abs/2005MNRAS.359.1083G

Stability Self Gravity: ​http://adsabs.harvard.edu/abs/2002A%26A...386..732C ​http://adsabs.harvard.edu/abs/2003astro.ph..7532F

Equation of state: ​http://adsabs.harvard.edu/abs/1986Ap%26SS.124..295M

Boyle's law and Gravitational instability: ​http://adsabs.harvard.edu/abs/1956MNRAS.116..351B ​http://adsabs.harvard.edu/abs/2001A%26A...375.1091L

Grav. Collapse of isothermal sphere: ​http://adsabs.harvard.edu/abs/1992AAS...181.5801F ​http://adsabs.harvard.edu/abs/1977ApJ...218..834H

Additional physics: ​http://adsabs.harvard.edu/abs/2001PhRvD..63h4005S ​http://adsabs.harvard.edu/abs/2011arXiv1105.5128J ​http://adsabs.harvard.edu/abs/2011A%26A...535A..49S ​http://adsabs.harvard.edu/abs/1996hell.conf..191B ​http://adsabs.harvard.edu/abs/1992Ap%26SS.193..173Y ​http://adsabs.harvard.edu/abs/1975MNRAS.172..441Y ​http://adsabs.harvard.edu/abs/1975MNRAS.171...85Y ​http://adsabs.harvard.edu/abs/1968MNRAS.140..109Y ​http://adsabs.harvard.edu/abs/1958MNRAS.118..523B ​http://adsabs.harvard.edu/abs/1984ApJ...276..737S ​http://adsabs.harvard.edu/abs/1988MNRAS.234..821O

• Patched ​Back Links Macro
  • See ProcessingObjects for an example using [[BackLinksMenu]]
  • While patching figured out how to modify wiki macros and recompile them etc… so I modified BackLinks to automatically be embedded in a collapsible structure instead of appearing next to the [[PageOutline]]
  • After you modify the source code in BackLinksMacro.py you have to:
    • compile the source python code to generate a .pyc file.
    • restart the web server to refresh the cached versions of the plugin

      johannjc@clover plugins > emacs BackLinksMenu.py -nw
      johannjc@clover plugins > python
      Python 2.6.6 (r266:84292, Sep 15 2010, 16:22:56)
      [GCC 4.4.5] on linux2
      Type "help", "copyright", "credits" or "license" for more information.
      >>> import py_compile
      >>> py_compile.compile("BackLinksMenu.py")
      >>> exit()
      johannjc@clover plugins > sudo apachectl restart
      

  • Also some plugins are packaged as egg files in which case you have to first unzip them - make the changes, compile, and then rezip them. Some egg files are not rezippable, but should have a setup.py script that you can pass to easy_install that should create the egg.
• I couldn't find an easy way to add [[BackLinksMenu]] to every wiki page, but did learn about page templates.
  • Template pages can be created in the Template sub directory. Currently there are three template pages
    • PageTemplates/Templates - for any random page
    • PageTemplates/TestTemplate - for test problems
    • PageTemplates/ProjectPage - for new project pages
  • When you create a NewPage you can select any of the existing templates as a starting point.
  • All of these templates have the [[PageOutline]] and [[BackLinksMenu]] macros already.

• Reopened ticket #127 after running higher res versions of the field loop advection.
• Added references to IICool to external links page as well as a link to the old astrobear site at the bottom of WikiStart
• Fixed carbuncles appearing in RadiativeInstability test problems #178 and updated the bear2fix.data.img file to show more than just the coarsest grid and modified the problem module to trigger refinement before the first step near the right edge.
• Made various code improvements
  • Added CellsPerJeansLength variable which allows users to adjust the jeans-length based refinement. Was previously a parameter = 4. Can be adjusted inside of a users problem module init - not in physics.data
• Duplicated the TestMacro that Matt wrote and had it use a directory on clover instead of the test repository when generating images etc… This enables folks to immediately see the output of the current test on the web before the code gets committed etc… Also added the test results from the last version of the code checked in to the official repository to AstroBearTesting.
• Pondered the Fedderrath particle creation tests
• Met with Eric to discuss GravoTurbulent sim params. Also have a 1024 core reservation that just started this afternoon.
• Running colliding flow sims on kure.

​movie

clump project

So I did some cooling shock problem using the clump settings, just removed the clumps. The shock front looked fine when hviscosity and periodic boundaries are used.
So the clump runs are being restarted from frame 0 using those settings. The following animations show Bx Strong, Bx Weak, By Strong, By Weak respectively. They are not finished yet, but are running. The final frame should be 60. So By runs are almost finished, Bx runs are half way through.
​Bx Strong

​Bx Weak

​By Strong

​By Weak

cloud collapse project

Eric and I thought it will be very interesting to do the BE collapse problem as mentioned in last week's journal club:

a paper about dynamo generations at various scales in the cloud collapse problem and the impact of a realistic prandtl number.

The runs will be similar to those of the Federrah 2011 paper, but put in physical viscosity, that is, Braginskii when gyrofrequency >> collision frequency, Spitzer when the opposite is true. The resistivity will be Spitzer, which we already have.

Based on these thoughts, I implemented Braginskii viscosity into the code. The numerical method is based on Parrish et al. 2012 paper:

​http://arxiv.org/abs/1201.0754

They did the HBI/MTI simulations with explicit Braginskii viscosity.

I then went on to test 2D Braginskii viscosity in a MHD KH instablitiy problem using James Stone's parameters:

​MHD Kelvin Helmholtz Instability with Braginskii type Viscosity
​MHD Kelvin Helmholtz Instability with Braginskii type Viscosity, Comparison

It is best to have the code switching between Spitzer and Braginskii viscosity based on the relation between local gyrofrequency and collision frequency. Unfortunately we do not have that yet. Not far away though.

code dev

A problem we discussed before: when considering the energy equation for the operator split induction equation, do we need to apply the induction flux to the energy equation explicitly as well? I now doubt to explicitly calculate the energy equation is the correct way. So I added some additional calculations to apply energy flux induced by multiphysical processes in both the braginskii viscosity and resistivity solvers.

I think we are at a point where sub-cycling is becoming a bottle-neck. Some thoughts are needed. The problem is:

How to do subcyclings in AMR without a significant increase in the size of ghost zones and at the same time maintain accuracy.

Just talked with Eric and we decided to try the following set of runs…

Constant params

• rho=100 particles/cc
• Mach = 5
• cells per jeans length 25 (or higher if computationally feasible)

Variable params

where

• II is Inoue and Inutsuka
• ISO is isothermal at 41 K
• All runs will start with windowed velocity field of fully developed turbulence.

To achieve a virial parameter of .2, we need to make the density higher or make the box larger… Adjusting the box size is the easiest thing as it will keep the temperature and the velocities the same. Before we were using R = 50 pc and a Mach number of 8.5

We also discussed a steady state gravo-turbulence where mass is taken from small scale collapsed structures and put back in at the large scale with or without additional velocity forcing.

Fully Developed Turbulence

Note for the isothermal run, only the mach number of the turbulence is important - everything else is scalable

For the IICool simulation the mach number, the mean density, and the box size all become important considerations…

So now we have turbulence initial conditions that we need to map onto a cloud

The Total Mass of the system in cu is 0.42836E+07 - and since the ambient has a volume of 100^3-(4/3*pi*203^{) and density of 1, this gives a cloud mass of approximately 3.3e6 or a mean density of 51 though it should be closer to 100.. It is likely less because of the softening… If we consider the radius to only be 80% then the mean density would be 99 particles/cc. There is still a net velocity of order 9 cu and a velocity dispersion (mass weighted) of 33. Given that the sound speed is of order 4.5 this is a mach number of 8}

For the isothermal run - the density can be arbitrarily scaled without effecting the dynamics… For the IICooling run - this is not the case. If should be possible to adjust the window to raise the min density…

Due to an accounting error… the mach number ended up not being the same for the two runs…

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

Particle creation requirements

1. violates the truelove criterion
2. is on the highest level of refinement,
3. surrounding flow is converging in each direction,
4. has a central gravitational potential minimum,
5. is Jeans-unstable,
6. is bound, and
7. is not within racc of an existing sink particle

1. and 2. are the same as Krumholz

3. seems reasonable to avoid particle formation in shearing layers

4. seems like it would favor particle creation in the center of the grid…

Consider a 1D system with two sink 'planes'. The Greens function is just \frac{|x|}{2} and let's consider three points at x=-1,0, & 1. The potential will then be \frac{|x+1|+|x|+|x-1|}{2}

As you can see only the central point would form a sink particle even though there are 3 equal masses. The paper is not entirely clear about the local gravitational potential and could possibly be referring to a gas potential calculated using only the local mass as opposed to the local gas potential of the entire domain.

Indeed, for check # 5 it calculates the self-energy of the control volume and does calculate a potential using only the local contributions… So they probably are referring to this same local-mass potential with regards to check 4. However, any cell surrounded by a uniform region that violates the true love criterion will have a minimum for the potential at the center of the region because the control volume will be centered on that cell. One could use the periodic approach where the mean density of the region is first subtracted off - in which case the potential would be zero everywhere… and would require the density to be the largest at the center… However, a better check - and one that wouldn't be dependent on the centering of the control volume might just be to make sure that the density is the highest in that cell compared to the surrounding cells in the control volume. This will keep particles from forming in very close proximity. We should update the code with one of these two approaches.

5. This requires that the surrounding control volume be as dense as the central cell - otherwise it will violate the true-love criterion without forming a particle… we have…

E_{grav}=-\frac{3GM^2}{5R_c}=\frac{16}{15}\pi^2G\rho^2R_c^5

E_{th}=\frac{3}{2}M\sigma^2=\frac{6\pi}{5}\rho R_c^3 c_s^2 where \rho \sigma^2=P

Jeans density criterion is that \rho > \frac{\sqrt{P \pi \gamma/G}}{R_c} which is equivalent to \rho > \frac{\sqrt{\rho c_s^2 \pi/G}}{R_c} \rightarrow R_c > \sqrt{\frac{\pi c_s^2}{\rho G}} = \lambda_J

Jeans energy criterion is that

2E_{th} < |E_{grav}| \rightarrow \frac{12\pi}{5}\rho R_c^3 c_s^2 < \frac{16}{15}\pi^2G\rho^2R_c^5 \rightarrow \sqrt{\frac{9}{4\pi^2}\frac{\pi c_s^2}{\rho G}} < R_c or that \rightarrow R_c > .4775 \lambda_J or in other words, if a region violates the density criterion, we would expect \frac{2E_th}{E_{grav}} = .4775^2

I've checked in the code that when the density is uniform and just above the jeans criterion that the ratio is .48^2 which is consistent with the various assumptions and discretizations… at least when gamma=5/3. When gamma = 1 this changes slightly to closer to .7^2.

So if there are only a few cells inside of the control volume near the jeans density, it may not satisfy the energy constraints initially… But it seems unlikely that it would artificially fragment without having enough mass present to satisfy this criterion.

After criterion 4, criterion 6 tends to be the most restrictive - as a collapsing region will not form a particle until enough kinetic energy has been converted into thermal energy at the center, and then cooled.

-Binary. I've ran 3 tests with a=40au. I have not seen the formation of a disk after 1 orbit. I suspect this happens because I'm using a ~4-times smaller resolution than before, because of the larger value of a. Here's an image of the run:

Left: Zoom in to see the secondary, showing the log density (cu) and the velocity field in Mach units. Material is pulled towards the particle but does not go about it. Right: vel field on the entire orbital plane. AGB is in the upper left while the secondary is on the bottom right. The latter does't seem to affect the wind strongly. Either the wind is too fast, 15km/s, too dilute, M'=1x10-5Msun/yr, or the gravitational pull of the secondary (1Msun, rsoft=2cells, 1283 + 4 particle refs.) is weak.

I'm now running a higher res version of it.

-Magnetic tower:.

• Andrea has replied, he's happy with the paper.
• I've calculated the volume time average mean flux ratios, <Q> in the following way. However, Eric and I agree that this volume calculation is not well defined for the fluxes; they should be taken on horizontal slices through the jet beam, thus we believe that Figure 7 in the paper (logarithmic flux rations color maps) does convey the fact that our towers remain, as they should, PDF from base to head, throughout the sims.

The volume calculations:

(1/N_cells) SUM_x,y,z [ f_P(x,y,z,t)/f_k(x,y,z,t) ], (1)

where x,y,z are cells inside boxes which contain:

• (i) the jet (RIGHT PANEL; x,y ⇐ |r_jet|; z⇐ height_cavity)
• (ii) the whole magnetic cavity (LEFT PANEL; x,y ⇐width_cavity; z⇐ height_cavity).

The time averages of the above:

SUM_t { (1/N_cells) SUM_x,y,z [ f_P(x,y,z,t)/f_k(x,y,z,t) ] }, (2)

where t = 42,84,118. These yield, <Q>_jet~ 21; <Q>_cavity~14.

-CRL618:.

• I've implanted a tracer field which takes the values of the energy loss due to cooling. Will give us a better idea on the emission of our nebulea.
• I've been doing low-res tests with a toroidal AGB (ambient medium) distribution using eq (13) from our PN paper, but with other, more torus-like, alpha and beta values. I've added the ambient shells to this as well. I see (i) concave, rather than a round, bow shocks, (ii) lobes with very similar, or smaller (fatter), aspect ratios than the ones produced with a spherical ambient. Not sure this is helping to get CRL618-looking objects.

-AGN jet truncation. Eric and I have discussed the relative speeds, densities and temperatures of the ambient, red giant (RG) stellar wind and the AGN jet. Here's our latest test run:

​https://clover.pas.rochester.edu/trac/astrobear/blog/agnJetTruncation

We study the 3D interaction of stellar winds from red giant (RG) stars and AGN jets. We're based on

¹​http://adsabs.harvard.edu/abs/2006MNRAS.371.1717H

Note this blog is in chronological order.

6 june '12

Latest sim ​http://www.pas.rochester.edu/~martinhe/2011/agn/hd-strongJet-edge2.gif shows that the red giant's (RG) wind, which crosses the edge of the jet's beam, strongly affects the collimation of the AGN jet. Note that in the previous sim ( ​http://www.pas.rochester.edu/~martinhe/2011/agn-6apr-dens2.gif) the RG's orbital path intersected the jet's axis, affecting the jet's beam more drastically.

9 Apr '12

The other parameters are as in Run 3 (below) except for:

• the jet is launched earlier
• the RG's initial position is farther form the jet's axis

No. density movie: ​http://www.pas.rochester.edu/~martinhe/2011/agn-6apr-dens2.gif

I see that a reconfinement shock forms in the jet, the height of which increases in time and, at this point, is at a height close to the RG path. These shocks are expected⁺ and seen in many other agn jet sims. This process enhances the jet's kinetic energy flux at that height, hence affecting the winds' interaction.

⁺​http://adsabs.harvard.edu/abs/1991MNRAS.250..581F

2 Apr '12

Because of the opening angle of 5^o , the jet has a density contrast of 0.0001 with respect to the ambient medium at the position of the RG, 1100pc from the origin

No. density movie: ​http://www.pas.rochester.edu/~martinhe/2011/agn-dens-2apr.gif

14 March '12

^* Because of the opening angle of 5^o , the jet has a density contrast of 0.001 with respect to the ambient medium at the position of the RG.

The jet will be launched once the star has reached the middle of the grid in the horizontal direction, not the other way around.

13 March '12

Meeting with Eric and Martín. We found, by taking the momentum ratio of the jet over the star on the interception surface, that the back force of the jet upon the star is very weak; several orders of magnitude. The passage of the star though the jet beam is ~ 10⁵ yr.

New runs:

• Run 2: as Run 1 (below) but (i) 8 times more resolution by a combination of grid size, jet radius and refinement level. (ii) The jet will be launched once the star has reached the middle of the grid, not the other way around.
• Run 3: as Run 2 but with L_jet-kinetic=10⁴⁵ erg/s and the corresponding jet density for vel_jet=.5c,
• Run 2.1: as Run 2 but with a magnetized jet with beta=1, the geometry of which will be decided later.
• Run 3.1: as Run 3 but with the magnetized jet of Run 2.1.

3 March '12

^* Because of the opening angle of 5^o , the jet has a density contrast of 0.001 with respect to the ambient medium at the position of the RG.

The sRG enters the edge of the grid with a horizontal velocity of 200km/s once the jet's head has reached a third of the grid along the vertical direction.

Run time=2.3days, 24procs, afrank, bluehive

Log No. density movie. The jet's bow shock quickly leaves the grid and affects the RG's wind spherical structure (bright small structure on the left). The simulation takes place inside the cocoon then. The jet is not truncated; it's quite strong relative to the wind, in agreement with¹. The interaction of the wind and the jet's beam edge affects the structure of the cocoon asymmetrically. A boundary effect rises affecting the jet beam structure near the base. It relates to the grid structure, easy to solve.

​http://www.pas.rochester.edu/~martinhe/2011/agn6mar.gif

Mach No. The shape of the RG wind is amorphous, not spherical, despite that fact that I do inject a spherical wind, for it is subsonic in relation to the cocoon material.

​http://www.pas.rochester.edu/~martinhe/2011/agn8mar.gif

Jonathan and I were trying to do some optimization on the refinement patches, as the smaller the patches are the more resources (memory and computing time) it needs:

​https://clover.pas.rochester.edu/trac/astrobear/blog/bliu02272012

The current algorithm works as the following: if the filling ratio is less than the desired ratio, it cuts the refinement patch into smaller ones until the filling ratio is larger than the desired ratio according to the inflections of ErrFlags.

In Figure 1 filling ratio of the patch (big box) is 40%. When the desired filling ratio is larger than 40%, AMR goes to smaller patches. We can see the running time increases when Desired filling ratio goes beyond 40% in Figure 2.

Figure 1

Figure 2

For refinement patch in Figure 3, however, the current algorithm couldn't find smaller patches even when the filling ratio is less than the designed one.

Figure 3

Figure 4