Thoughts on binary sims
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_{0}) between the wind leaving the primary at t=t_{0} and arriving at a position x at time t as t-t_{0}=|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})) = 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_{0}))
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
DO
Calculate the displacement vector from the primary at the retarded time
Calculate the wind normal so that
It is possible to solve
for the unit vector
Calculate the wind velocity from the primary
Update the retarded time using the new distance and wind speed
END DO
The only problem occurs when there are multiple solutions for the retarded time…
This will occur once we reach distances of order
If we switch to a rotating frame that rotates counter to the orbit so the angular speed is
, then
and
so that
Meeting Update 03.26
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):
Meeting Update
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 |
Baowei's Meeting Update for Mar 27 2012
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…
Martin's update, 27 mar '12
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.
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 |
---|
Meeting Update
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:
Tickets
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
Owners | Tickets |
---|---|
johannjc | #182 Spectra processing objects #179 Current test suite does not test Isothermal solver |
shuleli | #126 strange field behavior on amr edges with uniform gravity #121 3D, AMR, large sims abort possibly due to memory problem #151 Thermal Diffusion #152 Implementing Magnetic Resistivity #153 Implement Viscosity |
bliu | #183 Optimization for Grid Generation Algorithm #173 Interpolation options can trigger protections #71 Adaptive message block sizes #174 Point gravity and outflow properties stored in chombo #176 Adding additional tests #127 Oscillations in PPM MHD #150 Implementing Self Gravity in Cylindrical Coordinates #154 Sink Particles in 2D and Cylindrical #179 Current test suite does not test Isothermal solver #169 We should post links to papers using Astrobear on the wiki |
ehansen | #147 porting over 'NEQCooling' ?#171 If the initial conditions trigger protection we should stop the run and print an appropriate message |
erica | |
No Owner | #82 Create suppression file for valgrind's I/O errors #87 Improve parallel performance on Chombo HDF5 writes #92 Incorporate MPI_PACK_SIZE() into packing algorithm #155 Roe Solver |
A more detailed report can be found here https://clover.pas.rochester.edu/trac/astrobear/report/8
Radiative Shock Simulation
In this problem, you set the ambient conditions and use shock jump equations to solve for initial post-shock conditions. Then, the initial post-shock conditions become your initial boundary conditions for the ODE inside the cooling region. I've been using Mathematica to solve this ODE for the pressure profile. The profiles for other important hydro quantities follow from this. The analytic form that Mathematica gives is not pretty, but when plotted it looks qualitatively correct.
For a while, I was struggling to get my profiles to look like those in Delamarter '01. I realized that the solutions for the pressure, temperature, etc. profiles in the cooling region are strongly dependent on gamma. So I made gamma a free parameter again, instead of just using 5/3. One problem is that Delamarter makes no mention of what value he used for gamma. Also, he does not give values for the initial post-shock quantities, so there is no way to calculate what he used for gamma. The best I can do is guess by looking at the plots. I think I get close to the initial post-shock values with a gamma of 1.25. The profiles look qualitatively correct, but my cooling region is about 10 cells (2.5*10^{16} cm) smaller than that of the Delamarter simulation.
Another issue is that my solutions are not steady on the timescales stated in Delamarter. It should be relatively steady for at least 4000 years, but my solutions are not at all.
I see four possible reasons for these discrepancies:
1) There is something wrong in AstroBEAR
2) There is something wrong with the Delamarter simulations/calculations
3) The analytic formulas from Mathematica are wrong+
4) I am still not using gamma correctly++
+Note that I solved the differential equations in Mathematica analytically and numerically…the results are the same.
++Another thing that I haven't tried yet is using two different gammas. One gamma for the cooling region, and a different gamma elsewhere.
I have not yet tried using my own numerical ODE solver, so maybe that's the next step. I'll post some equations and plots before our next group meeting.
Meeting Summary
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
Meeting
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)^{3}
- 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?
Meeting Update 03/20/2012 - Baowei
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:
Meeting Update 03.19
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:
Martin's update, 20 mar '12
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_{1}/m_{2}=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?
Meeting Update 03.13
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.
Meeting Update 03/13/2012 - Baowei
Wiki Page for Ticket Assignment Procedure
https://clover.pas.rochester.edu/trac/astrobear/wiki/TicketAssignmentPage
Current Active Tickets
Total Active Tickets | 18 (include the two just closed) |
---|---|
astrobear | 15 |
bear2fix | 1 |
wiki | 1 |
Over 6 weeks | 4 |
Over 4 weeks | 8 |
Over 2 weeks | 10 |
Testing Results for new patch refinement algorithm
https://clover.pas.rochester.edu/trac/astrobear/blog/bliu03112012
Martin's update, 13 mar '12
Binary.
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. |
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:
- smaller grid and jet radius by a factor of 2 + more refinement levels,
- same as 1 but for a jet kinetic luminosity of 10^{45}ers/s (instead of 4x10^{42}),
- same as 1 and 2 but for a magnetized jet. Should have an effect on the RG wind.
Meeting Update 03/13/2012 - Eddie
Cooling Curves
I made a blog post last week with a plot of cooling curves to compare Pat's tables with Dalgarno-McCray. It's interesting how different the curve can be based on the electron density and ionization fraction.
My Library
I went through a bunch of old emails from Adam and Pat, and I made a list of papers that I want to take notes on to become more familiar with them (for my library on the wiki). It seems like a long list that could get longer, so I will probably focus more on the key papers. Adam, are there any papers in this list that you feel shouldn't be there?
Cooling Project
I updated my version of astrobear to try the changes that Jonathan had made with how the pressure is handled by the cooling source term. It changed a little bit, but there are still discrepancies between the cooling and adiabatic simulations. It appears as though it helped, the post-shock temperature in the cooling cases went up a little bit, but it's hard to tell due to the oscillatory nature of looking at frames that are several time steps apart.
So I've moved on to writing a 1D Radiative Shock module, and I will most likely have results from that next week. Just as a reminder, this follows what was done in AstroBEAR 1.0 (1D Radiative Shock test page) and Ch. 4 of Delamarter '01 (pdf). This module could also become part of the test suite if we want.
Meeting update
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?
Optimization: New Algorithm for Patch Refinement
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
Hypre memory usage
Level | nBoxes | nCells | x_{min} | y_{min} | z_{min} | x_{max} | y_{max} | z_{max} | MinBoxSize | MaxBoxSize | Table_{1} | Table_{2} | Table_{3} | TableSize | TableSize/nBoxes | TableSize/nCells | |
0 | 1024 | 262144 | 1 | 1 | 1 | 64 | 64 | 64 | 256 | 256 | 8 | 8 | 160 | 10,240 | 10 | 0.0390625 | |
1 | 246 | 157536 | 33 | 25 | 33 | 96 | 94 | 96 | 512 | 1344 | 20 | 21 | 8 | 3,360 | 13 | 0.021 | |
2 | 869 | 1132952 | 65 | 53 | 67 | 190 | 188 | 186 | 512 | 3584 | 44 | 46 | 24 | 48,576 | 55 | .042 | |
3 | 4242 | 8477904 | 131 | 107 | 135 | 378 | 376 | 370 | 512 | 20480 | 112 | 116 | 95 | 1,234,240 | 290 | .145 | |
4 | 4201 | 24779672 | 315 | 313 | 321 | 704 | 696 | 704 | 512 | 46080 | 192 | 190 | 189 | 6,894,720 | 1641 | .27 | |
5 | 10653 | 8818504 | 657 | 627 | 659 | 1380 | 1368 | 1362 | 512 | 4608 | 354 | 369 | 337 | 44,020,962 | 4132 | 5 | |
6 | 1902 | 1404112 | 1024 | 1024 | 1024 | 2728 | 2710 | 2642 | 512 | 2744 | 465 | 449 | 444 | 92,700,540 | 48738 | 66 |
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_{1} , Table_{2} , & Table_{3} 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.
Cooling Curves
I made this plot to look at how the Z cooling curve compares to the DM curve. Z cooling which, as a reminder, uses Pat's tables is highly dependent on the electron density ne and ionization fraction Xh. Here, I only plotted 2 possibilities: low ne, low Xh and high ne, high Xh. By "low" and "high" I mean the minimum and maximum of Pat's table.
Also, it is important to note that outside of the temperature range of Pat's table, the cooling routines currently revert back to neqcooling. It's more difficult to make a plot of this since there are no tables for this type of cooling. The cooling rates from ionization and recombination are calculated on the fly for neqcooling.
Digging through some literature on the physics of embedded clouds..
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
Weekly Summary
- 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.
- See ProcessingObjects for an example using
- 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.
- Template pages can be created in the Template sub directory. Currently there are three template pages
- 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.
Meeting Update 03.07
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
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.
GravoTurbulent Params
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
Run | alpha | Gravity | EOS |
A | 1 | SG | II |
B | .2 | SG | II |
C | 1 | SG | ISO |
D | 1 | Static | II |
E | .2 | Static | II |
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
Run | 128^{3} | 256^{3} | 512^{3} | Xmu | M_turb | v_turb | alpha | alpha t_cross/v_turb | t_cross | Temp | SoundSpeed | n | lBox |
ISO | 10 | 2 | 2 | 1.4 | 12.5 | 96.8272888743251 | 2.5 | 0.01333 | 0.516383352062004 | 20 | 4.472135955 | 100 | 50 |
IICool | 10 | 2 | 2 | 1.27 | 12.5 | 66.4996345290655 | 1.1792 | 0.01333 | 0.751883831 | 16.9812533471972 | 5.31997076232524 | 100 | 50 |
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*20}3^{) 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…
Run | 128^{3} | 256^{3} | 512^{3} | Xmu | M_turb | v_turb | alpha | alpha t_cross/v_turb | t_cross | Temp | SoundSpeed | n | lBox |
ISO | 10 | 2 | 2 | 1.4 | 21.5 | 96.8272888743251 | 2.5 | 0.01333 | 0.516383352062004 | 20 | 7.746183109946008 | 100 | 50 |
IICool | 10 | 2 | 2 | 1.27 | 12.5 | 66.4996345290655 | 1.1792 | 0.01333 | 0.751883831 | 16.9812533471972 | 5.31997076232524 | 100 | 50 |
Meeting Update 03/07/2012 - Eddie
finished resolution study
I've worked on completing the results for my 1D jet runs. The zcooling page now has data tables for the post-shock temperature for Z, DM, and no cooling. Here is a new plot that contains data from all the runs:
lineouts to look at shock structure
I've made some lineouts of density, velocity, pressure, and temperature to get a better understanding of how the simulation is interpreting the shock structure. Every image is taken at the middle of the simulation for the highest resolution run. The solid black line is density, dashed blue is velociity, dotted green is pressure, and dash/dot red is temperature. Everything is in computational units, and it's plotted on a log scale. The structure is more difficult to see in the cooling cases, so I included some magnified images for those two cases.
No Cooling
DM
Z
try new changeset, run old astrobear setup
This week I plan on using the newest changeset that Jonathan just pushed. It might have some affect on my simulations. Also, I will be making a new simulation that follows what was shown on the old astrobear website…http://www.pas.rochester.edu/~bearclaw/delamarterplots.html.
important papers
I've also been doing more reading. I think I will make a brief "library" of the most important papers on my user page similar to what Erica did. u/ehansen
Pondering the Fedderath particle formation requirements
Particle creation requirements
- violates the truelove criterion
- is on the highest level of refinement,
- surrounding flow is converging in each direction,
- has a central gravitational potential minimum,
- is Jeans-unstable,
- is bound, and
- is not within racc of an existing sink particle
- and 2. are the same as Krumholz
- seems reasonable to avoid particle formation in shearing layers
- 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
and let's consider three points at x=-1,0, & 1. The potential will then beAs 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.
- 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…
where
Jeans density criterion is that
which is equivalent toJeans energy criterion is that
or that or in other words, if a region violates the density criterion, we would expect
I've checked in the code that when the density is uniform and just above the jeans criterion that the ratio is
which is consistent with the various assumptions and discretizations… at least when gamma=5/3. When gamma = 1 this changes slightly to closer to .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.
Martin's update, 6 mar '12
-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:
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
AGN jet truncation
We study the 3D interaction of stellar winds from red giant (RG) stars and AGN jets. We're based on
^{1}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
Run 4 | AGN jet (base) | RG wind | ambient medium |
---|---|---|---|
resolution [cells; 1=9.6x10^{18} cm] | 32 cells | 6 cells | 64^{3} * 2^{3} |
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
Run 3 | AGN jet (base) | RG wind | ambient medium |
---|---|---|---|
density [part./cc] | .01^{} | 50 | 1 |
vel [km/s] | 1.5e5 (~c/2) | 200 | 0 |
L_{kinetic} [erg s^{-1}] | 1.9x10^{45} | … | N/A |
Orbital velocity [km/s] | N/A | 600 | N/A |
pressure [cu] | consistent with an opening angle of 5^{o} | ←same | ←/4 |
temperature [K] | 1.25x10^{10} | 2.5x10^{6} | 3.1x10^{7} |
radius | 200pc | 10AU (large?) | … |
mass-loss [M_{sun} yr^{-1}] | 4.8x10^{-2} | 1.9x10^{-3} | 0 |
resolution [cells; 1=4.8x10^{18} cm] | 64 cells | 8 cells | 64^{3} * 2^{4} |
^{} 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
Run 2 | AGN jet (base) | RG wind | ambient medium |
---|---|---|---|
density [part./cc] | 1^{*} | 10^{7} | 1 |
vel [km/s] | 1.5e5 (~c/2) | 200 | 0 |
L_{kinetic} [erg s^{-1}] | 4.7x10^{46} | … | N/A |
Orbital velocity [km/s] | N/A | 600 | N/A |
pressure [cu] | consistent with an opening angle of 5^{o} | ←same | ←half |
temperature [K] | 10^{10} | 1.3x10^{3} | 3.1x10^{9} |
mass-loss [M_{sun} yr^{-1}] | 1.1x10^{-2} | 1.5x10^{-6} | 0 |
resolution [cells; 1=4.8x10^{19} cm] | 64 cells | 8 cells | 64^{3} * 2^{4} |
^{*} 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.
Log No. density movie (still running): http://www.pas.rochester.edu/~martinhe/2011/agn-dens-20mar.gif |
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^{5} 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^{45} 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
Run 1 | AGN jet (base) | RG wind | ambient medium |
---|---|---|---|
density [part./cc] | .1^{*} | 1 | .1 |
vel [km/s] | 1.5e5 (~c/2) | 200 | 0 |
L_{kinetic} [erg s^{-1}] | 4.2x10^{42} | … | N/A |
Orbital velocity [km/s] | N/A | 600 | N/A |
pressure [cu] | consistent with an opening angle of 5^{o} | ←same | ←half |
mass-loss [M_{sun} yr^{-1}] | 1.1x10^{-2} | 10^{-4} | 0 |
resolution [cells; 1=62.5pc] | 16 (4 ref. levels) | 4 (ditto) | 32^{3 } |
^{*} 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^{1}. 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.
Optimization: Running Time VS Desired Filling Ratio for Refinement Area
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 |
---|