I patched the BackLinksMenu macro so now pages that link to a given wiki page will be listed in a collapsible section at the top of the page provided the page has [[BackLinksMenu]] at the top. Also created some new page templates for new pages, test problem pages, and project pages. See PageTemplates/Templates, PageTemplates/TestTemplate and PageTemplates/ProjectPage

Now when you create a new page you can select one of the templates as a starting point. To create new templates just create a page in the PageTemplates directory

So I was thinking it would be useful on the wiki if every page had a sort of directory heading that showed how one would get to the page? For example see ProcessingObjects

Also created a graph showing connectivity between different wiki pages. It's a little too large to see at once…

It looks like there are a couple of useful

I've started running some hires 3D simulations of gravitationally induced turbulence. Here are a couple images of the column density at 17 Myr

as well as a ​movie

I went with a cloud at a density of 10 particles/cc and a temperature of 160 K and a radius of 20 pc for a total cloud volume of about 8500 M_Solar. The cloud material initially has a Jeans Length of 80 pc. I was able to get pressure equilibrium by using the Inuoue and Inutsuka cooling/heating curve which gives an equilibrium temperature (pressure) for a given density plotted below. x-axis is log10 number density and y-axis is log10 pressure in units of (Kelvin/cc) (ignore the horizontal line)

This curve allows for gas at 10 particles/cc and 160 K to be in pressure balance with gas at .23 particles/cc and 6784 K.

It would be possible to raise the density to 100 particles/cc and still be in pressure balance with gas at 1 particle/cc if we wanted to make the cloud more massive.

The cloud is initially kicked with a solenoidal velocity perturbation with a k^-2 spectra and a total kinetic energy that is only about 1/3 the gravitational binding energy although there is plenty of internal energy to support the cloud

Since the Kinetic energy is 1/5 of the thermal energy, the initial perturbations are sub sonic. There is still structure that initially appears, but given the scales, I believe it is driven by thermal instabilities rather than gravitational.

Here are some histograms of density at approximately 8 Myr (red line) and 16 Myr (green line).

and density weighted histograms of the Jeans Length

Eventually the central part of the cloud collapses and forms 4 sink particles at around 17 Myr, while the rest of the cloud is still fairly whispy.

Given that the resolution is 100 pc / 2048 the code should form particles when the jeans length is .2 pc which should occur at densities of 1e5 with a corresponding jeans mass of about 10 solar masses. Given that the cloud contains 8000 M_Solar eventually there should form 100s of 10 solar mass cores - provided there are mechanisms for increasing the density to 1e5. Global contraction may be necessary to provide the density enhancement - as the cooling instabilities can only compress the material to 120 particles/cc or so…

In the colliding flows, we are embedding dense small clumps in the colliding streams to provide perturbations. We want these clumps to be in pressure equilibrium with the stream, and for the stream and the clumps to be in a thermally stable region (n < 1 or n > 10). It is also good if the clumps are large enough to be resolved, but be light enough to not plow through the opposing flow. The two stable densities allowed with the minimum density contrast occurs at arround 10 and .23 particles/cc. We decided to have a mean density of 1 particle/cc and a flow radius of 20 pc which gives a Mass flux of 523 m_sun/Myr for about 15 Myr or a total mass about 8000 M_Solar which would be comparable to the mass in the spherical cloud.

It would also be interesting comparing this to an isothermal run at 160 K or perhaps much lower… I was also thinking it would be good to always refine the Jeans length by 64 cells or so - until some maximum resolution. Also I've noticed that these runs are rotationally symmetric which is due to a difference in bluegene's compiler handling certain random number generating functions differently. That issue has been fixed.

• Binary: ​https://clover.pas.rochester.edu/trac/astrobear/blog/binaryBondi

• CRL618: We need to met and discuss what Adam and Bruce chatted about last week.

Preliminary look at a 2 clumps, one after the other, run made with periodic BC. The initial dens-contrast was 100 (twice as in the previous models, in order to have sensibly dense 2nd clump). The simulations in which I placed the 2nd clump into the grid, once the 1st one had left, generated cfl and nan errors which led to infinitesimally small dt. A wider grid needed. Yet, the lobe structure generated by the ambient shells-1st clump interaction is erased by the 2nd clump. Also, 2nd clump generated shocks should have quite low temperatures, may be hard or imposible to observe. ?

MOVIE

Now working on the run to produce a cooling map.

• Magnetic tower paper. Sent though for final revision. Some plots about the flux ratio, which I have not included in the paper:

• PN paper submitted to MNRAS

• AGN jet truncation. Managed to make the code cope with real high Mach numbers, order 100, for these runs. These are needed because of the large differences of temperature and velocity between the AGB wind and and the AGN jet. Meeting with Eric this Friday, 2:30pm, to discuss this. Test runs looking good.

Clumps:

See the project page for updated runs and images:

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

I have 6 jobs with 1000 cpu, 2days each that should start within this week (3 of them are on top of the queue). So we can for sure finish the uniform field runs this week. And maybe we can try some AMR runs to see if we can find the bug.

a) 1000 cpu may run better b/c of the larger memory.

b) the newer code in which Jonathan added some fixes last weekend may work better. (probly not related to the memory bug though)

c) do simplified runs to test the behavior: hydro only, or 1 level refinement.

Just a heads up that alfalfa ran out of space - and we only have 3.7 TB of free space left on the three current raid arrays. Here is a summary of current data usage:

The following summaries our understanding to the big ratio of Data memory over Data file size for AMR inspired by the results from Jonathan's memory checking tools — a typical number for the ratio could be 30~80.

Using AMR, extra ghost data are needed to do the interpolation. These ghost data size could be big when the refined patches are small.

Take a 3D problem with two-step updates for example,
Ratio $R=2*\frac{(m_{x}+16)*(m_{y}+16)*(m_{z}+16)}{m_{x}*m_{y}*m_{z}}$
where 2 comes from the copy we save for later restart and 16 comes from the ghost data.
So when m_{x}=m_{y}=m_{z}=8, we have

R=54.

Or the Data Memory is 54 times the size of the data file.

The smaller the AMR finer patches, the bigger the ratio is. So in Figure 1 we will have much smaller ratio than in Figure 2 though the total patches sizes are same. How these AMR patches distribute depends on the specific problem and calculation.

Figure 1

Figure 2

Ran a series of tests on the BE module. I found that I can get a stable sphere when running on my old data files, with the current version of astrobear I had been using; there was an issue with my data files. I then found I could get a stable sphere when running my current data files but at a lower resolution and smaller box size. This narrowed down the change in the data files to resolution/box size causing the instability. I am now trying to rule out what causes the instability (i.e. a larger box or a higher resolution or the combo of the the 2..?). It seems that perhaps the larger box in which the sphere sits may be throwing it out of equilibrium? Perhaps the BE solution can withstand only so much material outside of the sphere that is out of equilibrium…

Also, of possible interest.. in the runs that were collapsing (but that shouldn't have been…) The set up was initially in hydrostatic equilibrium, but seemed to remain in equilibrium throughout the collapse… ?? I made the vector plots of acceleration of (-gradP-gradPhi)/rho and saw vectors in the ambient pointing inwards (the gravitational acceleration), but inside the sphere all vectors canceled out — i.e. the forces balanced out inside the sphere. This behavior remained the same as the collapse happened. This is contrary to what I would think should happen. In the event of collapse, the inward gravitational acc vectors should dominate inside of the sphere. Something to perhaps look into at another time..

Sorry, likely won't be able to post movies or figures this week with homework, but will try to have a concrete diagnosis of what went wrong before…

fixed grid sim: ​http://pas.rochester.edu/~shuleli/cloudbyhigh.gif

AMR: run into the "fail to allocation" bug similar to that of ticket 121 and ticket 129 after about 20 hours.

Resistive: Implementing the tensor form.

I have been studying the resolution dependencies on the cooling routines. It seems as though I need a lot of cells per cooling length in order to achieve post-shock temperatures close to analytic predictions. However, I never quite reach the post-shock temperature that I expect. So I don't know if I am still under-resolved. This seems unlikely, because the post-shock temperature appears to be converging on some value, just not the value I expect. Maybe I am missing a fundamental difference between cooling and not cooling that I haven't accounted for.

I have kept diligent notes on this for the past couple of weeks. You can see more details on the zcooling page. The bottom subsection titled 1D Jet Simulations might be moved to a more permanent project page in the future.

I updated the development procedure page according to Jonathan's suggestion:
​https://clover.pas.rochester.edu/trac/astrobear/wiki/DevelopmentProcedure

As Adam asked, we will have two people (Baowei and Eddie) in charge of testing the code and checking in the test-passed code to the official repository. So whenever you have code you want to check in, just notify me or Eddie to do the test. If the test passed, we will upload the results to test repository and ask you/everyone else to verify. After that, we will push the code to the official repository. If the test fails, we will point you to the reference and simulation images as well as the necessary information to reproduce the failed test and leave it to you to determine why the test failed and to fix any possible bugs.

Check out the following pages for updates.

Binary. ​https://clover.pas.rochester.edu/trac/astrobear/blog/binaryBondi The central diks test + Bondi accretion

CRL618. ​https://clover.pas.rochester.edu/trac/astrobear/blog/crl618Figures (I've added movies and updated 2d color scale density maps for the ones with higher resolution). The run with 2 sequential clumps in bluegene's queue.

PNe paper for MNRAS. Orsola is reading it. I will submit tomorrow.

Magnetic tower paper. Adam is reading it.

I've sent the teragrid proposal.

I've resumed work on the AGN jet truncation project.

• Implemented Bondi Accretion… Updated SinkParticles with description of Bondi accretion routine… Someone should update documentation on SinkParticle
• Debugging Testing script… Making sure it works and catches failures etc…

More on resistivity, conduction implementation.
See "Multiphysics in AstroBEAR" project page.

More results of higher resolution Shocked clumps, as well as the resistive shocked clumps.
See "Shocked Clumps with Magnetic Field" project page.

Dennis et al.: ​http://adsabs.harvard.edu/abs/2008ApJ...679.1327D

New site with info for the paper: ​https://clover.pas.rochester.edu/trac/astrobear/blog/crl2013

17 dec

Short AGB rings separation, v_bullet=300km/s model, inclination angle of 30^o

TESTS Integrated "Doppler-shift" image. This is the bullet model at 200yr, 30o inclination angle, and only velocities from 10-50 km/s are mapped (based on a phone discussion with Bruce; I could map any other vel components if wanted).
Ditto with the "red-blue" Shape table.

Time [yr]	50	100	200
non-integrated density maps [part cm-3] and velocity field comparison Logarithmic false gray-scale density maps and color velocity field on the plane at the middle of the computational domain. Left panels show the bullet model, while right panels show the jet model. The bullet/jet are the densest, whitest, central features. Behind the bullet, or to the right of the jet (without loss of generality), we see the cavity formed by the bullet/jet, which is separated from the ambient medium by a contact discontinuity. The ambient medium is stratified and shows a series of concentric spherical shells (see Section "Initial conditions"). The vertical short red, or blue, lines at the top of the maps show the positions at which the lineouts of Figures X Y have been taken. The bullet and the jet propagate at 300km/s away from the ambient medium's densest region (up), forming an elongated lobe. Comparing the left and the right panels in FIGURE X [gray-scale density maps], we find that the lobes formed by the bullet and the jet are quite similar. As the bullet (jet) penetrates the ambient's shells, these form regularly separated vertebrae-looking features along the lobe. The axial velocity of these features decreases with radial distance from the axis.
non-integrated density lineouts [part cm-3] comparison. See top row for lineouts' positions.
non-integrated vel. lineouts comparison. See top row for lineouts' positions.
Cooling emission, slit marked in blue, Bullet
Cooling emission, slit marked in blue, Jet
PV diagrams, Bullet
PV diagrams, Jet

5 Dec 2012, long AGB rings separation, v_bullet=200km/s model

Time [yr]	50	100	200
non-integrated density maps [part cm-3] comparison
integrated Emission & non-integrated vel. comparison
non-integrated density lineouts [part cm-3] comparison. See top row for lineouts' positions.
non-integrated vel. lineouts comparison. See top row for lineouts' positions.

2 Apr '12

• Clump/Jet to ambient density contrast of 100 (was 50)
• Toroidal AGB wind, based in Frank & Mellema '94, alpha=.7 and beta=.8
• Higher initial densities
• 128x128x192+2AMR levels (~4 days to run in 1024 bgene procs)

There are some spurious numerical artifacts ahed of the jets' head at late times.

16 mar '12

14 feb '12

Clump/jet velocity of 400 km s^-1, dens_clump/dens_amb=50, larger grid.

Lineouts are taken along the lines marked in this figure >

Clump, stratified and ringed ambient medium:	t=88yr	t=176yr	t=247yr
Log(dens) maps. MOVIE: ​http://www.pas.rochester.edu/~martinhe/2011/crl/clump-ring-dens.gif
Log(dens) axial lineouts
Axial vel lineouts
Clump, stratified ambient medium:	t=88yr	t=176yr	t=247yr
Log(dens) maps. MOVIE: ​http://www.pas.rochester.edu/~martinhe/2011/crl/clump-dens.gif
Log(dens) axial lineouts
Axial vel lineouts
Jet, stratified and ringed ambient medium:	t=88yr	t=176yr	t=247yr
Log(dens) maps. MOVIE: ​http://www.pas.rochester.edu/~martinhe/2011/crl/jet-ring-dens.gif
Log(dens) axial lineouts
Axial vel lineouts
Jet, stratified ambient medium:	t=88yr	coming soon	ditto
Log(dens) maps. MOVIE: ​http://www.pas.rochester.edu/~martinhe/2011/crl/jet-dens.gif
Log(dens) axial lineouts
Axial vel lineouts

SYNTHETIC IMAGES

Ambient =r^-2 + rings, time~176yr, the slit is r_jet/2 displaced from the symmetry axis

Inclination	Clump, log(rho2)	Clump, pv		Jet, log(rho2)	Jet, pv
90o
60o

Ambient =r^-2, time~176yr for the clump and time~88yr for the jet (still running), the slit is r_jet/2 displaced from the symmetry axis

Inclination	Clump, log(rho2)	Clump, pv		Jet, log(rho2)	Jet, pv
90o
60o

Bruce's thoughts

RED IS SHOCK-EXCITED [NII] GREEN IS HALPHA — MOSTLY COMING FROM THE STAR BLUE IS MOSTLY SCATTERED STARLIGHT ============== CO observations regarding the ambient medium

take a look at Fog 2 and the rest of Sa ́nchez Contreras, C., Sahai, R., & Gil de Paz, A. 2002, ApJ, 578, 269 (optical emission-line study of CRL618): ​http://iopscience.iop.org/0004-637X/578/1/269/pdf/55888.web.pdf)

Here's my summary of the highlights of this paper

Fig 2 shows long-slit P-V diagrams of selected emission lines. Scattered stellar Halpha is significant near the center, and the Halpha line is thermally broadened everywhere. The other lines show a linear rise in V with slit position from the star…close to a Hubble flow. But the lines don't arise from the insides of the fingers. They come from an outer sheath where the lateral expansion of the fingers produce shocks. (Inclination corrections are important in the interpretation of these spectra.) (Note: one might expect two velocity components at each position arising from the near and far sides of the laterally expanding fingers. This isn't seen, but the limited dispersion of the spectrograph could be the reason why. Our models need only explain the general trends of the P-V diagrams.

The optical fingers are relatively low- and constant-density outflows (~5000 cm-3) plowing through a denser and slower AGB wind (10^{6 cm-3; 18 km/s). The density of the AGB wind appears to decline as r}-2. The shock temperatures (shock speed) derived from emission lines range from 10,000 to 25,000K (75-200 km/s). (They do not directly measure the temperatures along the fingers and through the tips using emission line diagnostics such as [NII]5755/6584.)

The optical emission in the fingers consist of (1) scattered starlight, (2) scattered Halpha and other emission lines (e.g. [OIII]) from a small, dense (10^{6 cm-3, T~14000K) HII region within 0."4 of the core, and (3) shocked gas at the interfaces along the edges and at the tips of the fingers. (Scattered light from the HII region is seen inside the lobes but at different velocities than the intrinsic emission from the lobes.)}

The shock speeds needed to explain the observed emission-line ratios are 75-200 km/s. An independent measure of v_shock at the fingers tips comes from the width of emission lines at the tips of the fingers: 200-230 km/s. The shocks can cross the fingers laterally in a few years.

The inclination angle of the fingers is 24 ± 6 deg, so the predicted space motions of the tips are 80/tan24 = 180 km/s, in good agreement with the estimated shock speeds

The shocks along the edges of the fingers are radiative (detailed discussion in section 9 p 286). The cooling time behind the shock is weeks. So energy from an ongoing fast wind is needed in order to maintain the radiation in the finger edges

The measured density along the length of the fingers is 5000 cm-3 with only small and local variations. This constant density contrasts with the density distribution of the surrounding neutral gas.

(BB: question for the Rochester crew: according to models, will the jets accelerate as they plow through decreasing ambient density? About knots, how will the general shapes of the fingers differ in an external medium that is (1) uniform, (2) declining as r^{-2? Can we use the shapes of the fingers to deduce the density falloff in the external medium?)}

============== More recent CO observations with somewhat higher spatial resolution and a spatio-kinematic model of the ambient molecular cloud

============== (movie 1)

============== Martin, two of the most notable features of the evolving structure of CRL618 are the motions and brightness changes of the tips of the fingers. This is seen on the 3-frame movie that's attached. open it with a browser. The movie frames are from 1998, 2002, and 2009 (not equally spaced)

One thing that we have yet to explore with the models is the changes in tip brightness as the tips cross the rings in the dust distribution. One might expect them to brighten as the bullets or jets encounter rings of higher density and to fade in between. In fact, after I stare at the movie, I think that I see the reverse…the tips fade when they hit the rings. Of course, the enhanced foreground extinction will affect the optical brightness.

Too bad that we don't have an IR move of the fingertips!

movie 2

============== Balick's first thoughts on the models 92/10/12)

FILES:

• ​http://www.pas.rochester.edu/~martinhe/2011/crl/PastedGraphic-1.tiff
• ​http://www.pas.rochester.edu/~martinhe/2011/crl/CRL618residmovie.gif
• ​http://www.pas.rochester.edu/~martinhe/2011/crl/crl618.f656-606-656movie00.gif
• ​http://www.pas.rochester.edu/~martinhe/2011/crl/CRL618.f547m.2009.6-1998..pdf
• ​http://www.pas.rochester.edu/~martinhe/2011/crl/Hydro_comments.pdf

Older model info

• ​https://clover.pas.rochester.edu/trac/astrobear/blog/author/blin
• ​https://clover.pas.rochester.edu/trac/astrobear/blog/syntheticImages
• ​https://clover.pas.rochester.edu/trac/astrobear/blog/martinhe01202012

Embedded disks in the AGB wind. I've ran two sims:

1. a disk contained in the orbital plane,
   • log density movie (upper part of the orbital plane, initially the disk is on the left and the AGB on the right): ​http://www.pas.rochester.edu/~martinhe/2011/binary/binary-n-disk1-dens.gif
   • vx movie (same display setup as the dens movie above): ​http://www.pas.rochester.edu/~martinhe/2011/binary/binary-n-disk1-vx.gif

I've only glanced at these movies and data, so the following are simple observations which need some serious thought.

• The structure of the disk is quickly lost
• The disk's former gas evolves and adoptes a structure which looks like the tilted disks of the binary sims (see table below).
• There is a shock wave coming from the AGB. This happened because I didn't adjust the density to match the one corresponding to the low resolution I used for this test. Yet, seems the diks looses its structure before the shock reaches the secondary's radius. I'll fix this for future tests, if any.

2. a disk with an ang. mom. vector at an angle of 30^o with the orbital ang. mom. vector.
   • log density movie (upper part of the orbital plane, initially the disk is on the left and the AGB on the right): ​http://www.pas.rochester.edu/~martinhe/2011/binary/binary-n-disk2-dens.gif
   • vx movie (same display setup as the dens movie above): ​http://www.pas.rochester.edu/~martinhe/2011/binary/binary-n-disk2-vx.gif

29 Mar '12

a=40 AU tests.

About run speed. It is proportional to:

1. the resolution

1.1 dx_max / dx_min

2. the filling ratios. They are not crucial for these sims because I'm using particle refinements, not to the AMR (which will not really help to seed up the runs). Yet, I control the radii of cells of each level, taken form the particles' center. I've adjusted these to be small enough to aid the run speed but large enough to capture the AGB wind and the Bondi radius.
3. dt, which is proportional to cfl*dx_min/ max(v_w, c_s, v_orbit)
4. the No. of processors used, N_p
5. cluster communication and processors' clock speed
6. the advance speed (computation time per grid level).

My tests indicate that for separations ≥ 20AU, the resolution should be >~ .047AU (64³+5refs) to form disks. I consistently see that bhive's runs are faster than bgene's (even when I use four times more processor in the latter than in the former). What worries me is the advance speed. The code's output shows that the simulation spends a good deal of time in this process. e.g. test 6 (column 6 in the table below) shows advances between the levels 4, 5 and 6 (which go back and forth several times between each level 0 dt) take as long as 12 secs. So every full timestep advance takes ~minuts to happen. Why, and how can I improve this?

r_B=2Gm₂/(v_w² + c_s² + v₂²) [from ^*],

r_B^'=Gm₂/(v_w² + c_s²) ,

where r_B, m₂, v_w, c_s, and v₂ are the Bondi radius, the secondary's mass, the sound speed and the secondary's orbital velocity with respect to the center of mass.

Test	q=m1/m2	resolution	rsoft/dx	rB/dx	rB/rsoft	rB'/dx	rB'/rsoft	Tw [K]	vw [km/s]	M'w [10-5Mo/yr]	time [orb]	log(dens/cu),vel/Mach
1 (bgene, "Rsoft8")	1.2/.6=2*	642x16+6refLev, dx=.023AU	2	448.7	224.3	272.1	136	300*	15	1	.016
2 (bgene, "2", running)	1.2/.6=2*	1282x32+6refLev, dx=.011AU	2	448.7	224.3	544	272.1	300*	15	1	.002
3 (bhive)	1.2/.6=2*	642x16+7refLev (slightly larger grid), dx=.015AU	6	718	119.6	435.3	72.6	300*	15	1	.102
4 (bhive, "2", running)	1.2/.6=2*	1282x32+5refLev, dx=.023AU	3	448.7	149.6	272.1	90.7	300*	15	1	.061
5 (bgene, "3", running)	1.5/1=1.5+	1282x32+5refLev, dx=.023AU	2	693	346.6	424.1	212.1	1000	15	10	.008
6 (bgene, "4", running)	1.5/1=1.5+	1282x32+6refLev, dx=.011AU	2	1386.4	693.1	848.2	424.1	1000	15	10	.001	< tilt
AGAIN 6 (bgene, "4", running)											.0027a
7 (bhive, "3", running)	1.5/1=1.5+	642x16+5refLev, dx=.047AU	4	301.1	75.3	179	44.7	3000	9	2	.062

^* val-Borro et al.

⁺ M&M '98

^a ​http://www.pas.rochester.edu/~martinhe/2011/binary/gene-4.gif

20 mar '12

As in M&M '98, model 1. I did a similar test before (see 21 feb '12 post), but it was shorter, with less resolution and with a separation

• Temp=3000K
• AGB mass-loss = 10^-5 M_o yr^-1, from time=0
• vel_wind=9 km/s, from time=0
• a=40AU
• circular orbit
• q=1.5 (m₁=1.5; m₂=1 M_sun)
• r_soft=2
• 64x64x16cells + 5 particle grid refinements
• grid: x,y ⇐|32AU|; z⇐|8AU|.

Running the same model but with r_soft=0.

As in Val-Borro et al. '09 (​http://adsabs.harvard.edu/abs/2009ApJ...700.1148D)

• Temp=1000K
• AGB mass-loss = 10^-5 M_o yr^-1, from time=0
• vel_wind=10 km/s, from time=0
• a=40AU
• circular orbit
• q=2 (m₁=1.2; m₂=.6 M_sun)
• r_soft=2
• 64x64x16cells + 5 particle grid refinements
• grid: x,y ⇐|60AU|; z⇐|15AU|.
• Running time (24 afrank p, bluehive) = 4 days

Log(dens) on orbital plane

Here are some snapshots of the orbital plane velocity field, ​http://www.pas.rochester.edu/~martinhe/2011/binary/40au-bb5-vel.pdf

Here are some plots of the mean angular momentum in the grid as a function of time, ​http://www.pas.rochester.edu/~martinhe/2011/binary/mAngMom.png

28 feb '12

As in M&M '98 model 3:

• Temp=2400K
• AGB mass-loss = 10^-5 M_o yr^-1 (they actually use 8x10^-6)
• a=17AU
• vel_wind=15 km/s
• circular orbit
• q=1.5
• r_soft=2
• 128x128x64cells + 4 particle grid refinements
• grid: x,y ⇐|25AU|; z⇐|12.5AU|; runs faster. I've checked and there's no boundary inflow.

Density iso-contours, edge-on, time=.5 orbits, zoom in. No tilting. The captured wind does show a disk-like shape, but at this point the disk is forming.

24 feb '12

Same model as yesterday, time=2 orbits. There's some inflow from all the boundaries onto the grid (left panel, small v_z from all boundaries), so I need a larger, or a cubic, grid and/or wind objects to block the inflows. Easy to solve. The inflow makes the flow patten quite messy and destroys the spiral structure induced by the orbit (middle panel, orbital plane view). The disk is not strongly affected by the inflow, which is good (right panel, low-res perspective view of iso-contours showing that the disk is plane symmetric and has an asymmetric bow shock which faces the AGB star).

I'll stop this simulation and send another one with the same setup but: a cubic grid and outflow-only wind objects at the boundaries.

23 feb '12

A very handsome disk!!'''

• Temp=3000K (as in M&M '98)
• AGB mass-loss = 9.9x10^-6 (as in M&M '98)
• a=3.7AU < M&M '98
• vel_wind=9 km/s (as in M&M '98)
• circular orbit
• q=1.5 (as in M&M '98)
• r_soft=2
• 128x128x64cells + 4 particle grid refinements⁺
• grid: x,y ⇐|1|; z⇐|.5|. Runs faster. We can use the full cubic grid for production runs.

Next steps:

1. try M&M '98 model 4.
2. elliptical orbit

⁺ Runs faster (1 orbit per 16 running hrs) than 32x32x16+6 particle refs (based on another test using 32³ cells +5 particle refinements which produced 1 orbit per 17 running hrs; yet for the latter test I used r_soft=0 which produces high velocities in the central-most orbits of the disk, hence reducing the timstep); the more grid refinements I request, the grater the numer of sub-steps the code has to perform. This is for 48 afrank processors.

22feb '12

Time=2.7orbits, no tilt.. The secondary keeps accreting mass and capturing wind material. I still see significant density gradients in the bound gas. Mastrodemos & M do not see this, so I think this is purely a shock effect. De Val-Borro et al. do see something similar, specially for their large wind region models (​http://adsabs.harvard.edu/abs/2009ApJ...700.1148D), but this may be a 3D effect. This image shows 5 flow steam lines of the wind captured gas.

I want to run this problem again with 1 more refinement level and a r_soft=2dx which should reduce the velocity of the central most cells. I presume this will make the disk radius larger.

21 feb '12

As in M&M '98 model 1:

• Temp=3000K
• AGB mass-loss = 9.9x10^-6
• a=3.7AU
• vel_wind=9 km/s
• circular orbit
• q=1.5
• 32³ + 5 particle grid refinements
• r_soft=0

I see a disk-looking structure (with a larger radius than the one in the previous test), with a Keplerian-like vel distribution and a flow structure that's symmetric with respect to the orbital plane (i.e. no tilting), after~1.3 orbits. The simulation is still running. I want to see it at 5 orbits.

Log density in grayscale. Velocity field in color scale in Mach units.
Zoom in to the disk

Central disk + Bondi:

I do not see tilting in the central region of the disk. Note that for this setup (taken from test 2*) rsoft=4dx, hence the particularly soft-looking center of the disk. This is about 25% of the run time that I used in test 2*, but it is well within the tilting growing time.

^* ​https://clover.pas.rochester.edu/trac/astrobear/blog/martinhe02032012

20 feb '12

As in M&M '98:

• Temp=2400K
• AGB mass-loss = 9.5x10^-6
• a=12AU (they actually have 12.6 AU)
• vel_wind=25km/s
• circular orbit
• q=1.5
• 32³ + 4 particle grid refinements
• r_soft=0

I see a disk-looking structure, although smaller than I expected, after~7.5 orbits. Again, it is symmetric with respect to the orbital plane, hence Jonathan's Bondi accretion implementation see to have solved the tilting enigma. I need more resolution to capture the full formation of the disk.

15 feb '12

Same as below but 64³ + 5 particle grid refinements and temp_wind=1000K (and in an exploratory elliptical orbit)

PRELIMINARY: I see a disk-looking structure after~3.7 orbits which is symmetric with respect to the orbital plane!!

14 feb '12 Jonathan has implemented, and tested with a central disk problem setup, the Bondi accretion of sink particles. I've been trying it in one instance of my binary-formed disks problem with:

• isothermal solver, gamma= 1.001
• r_soft=0 (as suggested by Jonathan)
• a=20 AU
• q=1.5
• vel_wind=20km/s; temp_wind=2000K (line M&M '97); mass-loss_wind=10^-5 M_sun /yr
• 64³ + 3 particle grid refinements, so dx~2.1 AU
• r_Bondi=4 AU; r_Bondi/r_soft=NA

I do not see a disk even after 20 orbits. The wind captured structure does look different than the one formed without* the Bondi accretion. I.e. it shows a light core and a dense tail (see left and upper right panels), but not the other way around (as in * ). Also, the structure is symmetric with respect to the orbital plane (see lower right panel).

I'm now trying this setup, but with these changes:

• 5 particle grid refinements ( not 6 yet for it will be rather slow), so that r_Bondi=8dx (instead of 2dx)
• temp_wind=1000K, instead of 2000.

^* see test 1 in ​https://clover.pas.rochester.edu/trac/astrobear/blog/martinhe01292012

3 runs:

Test for stability - ambient 10 x lighter:

Radial velocity lineout: ​http://www.pas.rochester.edu/~erica/rhoweight10_vrad.gif

Log density lineout: ​http://www.pas.rochester.edu/~erica/rhoweight10_rhoLineout.gif

Log density: ​http://www.pas.rochester.edu/~erica/rhoweight10_rho.gif

Pressure: ​http://www.pas.rochester.edu/~erica/rhoweight10_pressure.gif

Test for stability - ambient 100 x lighter:

Radial velocity lineout: ​http://www.pas.rochester.edu/~erica/rhoweight100_vrad.gif

Log density lineout: ​http://www.pas.rochester.edu/~erica/rhoweight100_rhoLineout.gif

Log density: ​http://www.pas.rochester.edu/~erica/rhoweight100_rho.gif

Pressure: ​http://www.pas.rochester.edu/~erica/rhoweight100_pressure.gif

Test of Collapse:

Radial velocity lineout: ​http://www.pas.rochester.edu/~erica/rhoweight1vrad.gif

Log Density Lineout: ​http://www.pas.rochester.edu/~erica/rhoweight1rhoLineout.gif

Log Density: ​http://www.pas.rochester.edu/~erica/rhoweight1.gif

CRL618

14 feb '12

Clump: 400 km s^-1, dens_clump/dens_amb=50, larger grid

Ambient =r^-2 + rings, time~176yr, the slit is r_jet/2 displaced from the symmetry axis

Inclination	Clump, log(rho2)	Clump, pv		Jet, log(rho2)	Jet, pv
90o
60o

Ambient =r^-2, time~176yr for the clump and time~88yr for the jet (still running), the slit is r_jet/2 displaced from the symmetry axis

Inclination	Clump, log(rho2)	Clump, pv		Jet, log(rho2)	Jet, pv
90o
60o

9 feb '12

Clump: 400 km s^-1, dens_clump/dens_amb=50

Model C^*: ambient =r^-2 + rings, time~120yr

Clump	log(rho2)	pv
90o
60o

Model C^*: ambient =r^-2 + rings, time~120yr

Clump+ambient	log(rho2)	pv
90o
60o

^*​https://clover.pas.rochester.edu/trac/astrobear/blog/martinhe01202012

So i checked in what should be a newer version of the code to clover:/data/scrambler_devel that should be tested and (if successful) pushed onto clover:/data/scrambler.

Before testing on clover, someone should run hg update to update the source files to match the newly pushed head. I also was looking through all of the scripts and found…

• buildtests.s - Brandon's old suite of three tests intended to be run on bluehive.
• testbuilds.s - My script that attempts to compile every problem module
• test_scripts/test_Xproc.pbs - pbs submission scripts for Brandon's old suite of tests.
• modules/BuildModuleControl.pl - empty… I believe at one point Brandon or Matt intended to use this as a type of configure script to turn on/off the compilations of various modules
• tests/go.s - script that does post processing on an individual test problem passed as an argument
• postprocess.s - script that cycles through a list of tests and runs go.s (currently only setup to test FieldLoopAdvection)

But i currently don't see Matt's script or any script that would cycle through the problems and actually run them - or update the output on the wiki in astrobear in the scrambler repository

I did find an svn repository cloverdata/repositories/tests/ that has the images and the logs from the last test run that show up if you use matt's wiki plugin

​https://clover.pas.rochester.edu/trac/astrobear/wiki/u/noyesma#no1

And there is a copy of the repository at cloverdata/tests in which there is a process_chombos script which I imagine is supposed to update the test logs etc… and then check in the new images to the tests repository so they can be displayed on the wiki… but there should be another script somewhere that actually runs each test problem before calling this processchombos script…

Update

So the buildproblem script can be found in /home/noyesma/bin and Matt's screencast is available at

​http://dl.dropbox.com/u/2018623/Test%20Suite%20Tutorial.mov

Here's a movie of the "hydrostatic" disk. The hydrostatic solution is modified so that it never drops below some background density (when in practice it would need to to be in hydrostatic equilibrium)

movie​

Also redid the rotating true love collapse problem

movie​

Presentations
The LLE presentation is at:
​http://www.pas.rochester.edu/~shuleli/magclumps2.pdf

Clumps

The beta map of a bx case, vs a by case.

We can see there is a wake of low beta, low density region hiding behind the clump in the bx case. However, clump material is totally exposed to the erosion of the shock.
In the by case, there is a cylindrical cavity filled with strong magnetic field (low beta) surrounding the clump's "tail". Clump material are partly protected by this field. The field covering the boundary flow region is amplified by stretching. There is also an amplification region at the head of the clump initially, due to compression. Stretching is more effective in amplifying the field in this scenario. This is partly why the contained field case does not do much in terms of evolution. The only working field amplification mechanism is the compression by the transmitted shock, which is slowed down by a factor of 10.

The following shows the comparison of by cases between resistive and nonresistive runs at approximately 1 crushing time. The resistivity is isotropic and constant. We have microphysical resistivity which depends on the number density and the temperature, and nonlinear. We may want to try that out later.

We can see some blueish region disappeared, and the region behind transmitted shock has a higher beta. The clump tail is less obvious in the resistive case. The field diffusion speed is 2 orders of magnitudes smaller comparing to the shock speed. We can expect things getting more different when the resistivity is higher. Then again the problem is "is the realistic microphysical resistivity high enough to trigger noticeable effect".

The high res runs (50 zones per radius) have been sitting on the bluegene queue for five days.
It's currently near the top of the queue. The three queued jobs are bx strog, by weak, by strong respectively. Hopefully we can finish these plus the "single mode" contained field case within this month. Since these runs need to have a ending time twice as the AAS runs (we need at least 4 crushing time), I expect at least 48 hours on 512 processes for each run.

Multiphysics

Flux limited thermal conduction. Saw strange behavior while doing the test on the Parrish&Stone paper. Switched to the simple "ring test", found that the inner boundaries are doing some strange things when run in parallel. See the following picture.

Future Meetings, Abstracts?

My 1D Jet simulation now runs literally in 1D with Jonathan's recent 1D implementation. This makes the simulations much faster ~1 minute on alfalfa. I have run the simulations with the new cooling routines both with and without magnetic fields. I'm now working on getting the emission maps for H-alpha and SII. However, the emission routines in bear2fix are proving to be more troublesome than we initially thought.

I get results for H-alpha that look ok. I'm not sure if these are exactly right, but at a glance they make some sense. I am not getting anything for SII. So far I have traced the problem all the way back to how it calculates the temperature. It seems that some of the chombo information is not being read in properly. This is likely due to the chombo being in 2D, but the simulation and many of the arrays being allocated for 1D. Or there is some other bug that I can't locate. All I know for sure at this point is that TempScale is wrong. I have set TempScale to be 1000, but bear2fix thinks it is some crazy number ~10^-4. This makes input temperatures way too low to get any emission from anything.

Here are the movies of my 2D BE run:

​http://www.pas.rochester.edu/~erica/rhoFeb62012.gif

​http://www.pas.rochester.edu/~erica/rhoLineoutFeb62012.gif

​http://www.pas.rochester.edu/~erica/vRadLineoutFeb62012.gif

In the process of getting 3D runs posted. Currently a job is in Bluehive queue.

Jonathan and I are piloting a shared file system called Box.net. Currently it has 2GB single-file-limit. So it's probably better for documents rather than data files. It can sync your documents on your computer while keeping historical versions. I find it's convenient when two or more people cooperate on same documents. If you are interested or have better ideas of using it, please let me know.

​http://www.pas.rochester.edu/~bliu/Box.net/Box.net.png

Quotations for New Machine (​https://clover.pas.rochester.edu/trac/astrobear/blog/johannjc01182012)

ASA: 2 Xeon 2.4GHz Quadro Core Processors, 24GB Memory, 16TB Harddisk ​https://www.pas.rochester.edu/~bliu/ComputingResources/ASA_Computers.pdf

AberDeen: 2 Xeon 2.4GHz Quadro Core Processors, 24GB Memory, 16TB Harddisk ​https://www.pas.rochester.edu/~bliu/ComputingResources/Aberdeen.pdf

Pogo: 1 Xeon 1.6GHz Quadro Core Processor, 3GB Memory, 14.5TB Harddisk ​https://www.pas.rochester.edu/~bliu/ComputingResources/Pogo_linux.pdf

Current Load of my Teragrid allocation ​https://www.pas.rochester.edu/~bliu/ComputingResources/Teragrid_Load.png

• Cleaned up sweep scheme - and self gravity (soon to be checked into main repo)
• Implementing Bondi accretion
• Improved protection options. See ProtectionOptions
• Cleaned up cruft in PhysicsData namelist…

  !================================================================================
  ! AstroBEAR physics.data file.
  ! Input data describing built-in physics-related parameters of the code.
  ! Created for problem:  Template
  ! Default values for variables in [brackets]
  !================================================================================
  
  &PhysicsData
  
  !================================================================================
  ! Equation of State stuff
  !================================================================================
  gamma =  1.66666666666667  ! Adiabatic index for ideal gas [1.66666666667]
  Xmu   =  1d0               ! mean molecular weight [1d0]
  iEOS  =  0                 ! [0]-ideal gas, 4-Isothermal  (Equation of state)
  
  !================================================================================
  ! Field based refinement control
  !================================================================================
  InterpOpts = 0,0,0,0,0,0,0,0,0,0,0  ! [0]-constant, 1-minmod, 2-superbee, 3-vanLeer, 4-mc, 5-Parabolic (not conservative), 6-Linear (not limited).  It is recommended to use constant interpolation for any fields in q that correspond to aux fields (iBx, iBy, iBz, etc...)
  refineVariableFactor = 1d0,1d0,1d0,1d0,1d0,1d0,1d0,1d0,1d0,1d0,1d0 ! weight factors for each field used in triggering refinement.  Combined with qtolerance for that level. [1d0]
  
  !================================================================================
  ! Global source term switches
  !================================================================================
  lSelfGravity  = .false. ! Turns on self-gravity if true [.false.]
  UniformGravity = 0d0    ! Gravitational acceleration in the -y direction [0d0]
  iCylindrical  = 0       ! [0]-No cylindrical geometry, 1-Cylindrical with no angular momentum, 2-Cylindrical with angular momentum
  
  !================================================================================
  ! Density Protection Options
  !================================================================================
  lRestartOnDensityProtections = .false.  ! Do density protections trigger restarts?  [.false.]
  iDensityProtect  = 2                    ! 0-Set to MinDensity, 1-Set to minimum nearby density, [2]-Average nearby densities
  iMomentumProtect = 2                    ! 0-Conserve momentum, 1-Set to zero, [2]-Average nearby velocities
  MinDensity       = 1d-10                ! Minimum computational density before protection is triggered [1d-10]
  
  !================================================================================
  ! Pressure Protection Options
  !================================================================================
  lRestartOnPressureProtections = .false. ! Do pressure protections trigger restarts? [.false.]
  iPressureProtect = 4                    ! 0-Set to MinTemp, 1-Set to minimim nearby pressure, [2]-Set to average nearby pressure, 3-Set to minimum nearby temperature, 4-Set to average nearby temperature
  MinTemp          = 1d-10                ! [1d-10] minimum allowed temperature for the system in Kelvin before protection is triggered
  
  !================================================================================
  ! Description of various scaling parameters
  ! Define one of each of the following: [nScale/rScale], [TempScale/pScale], and set the other to 0d0.
  !================================================================================
  nScale          =       1d0,              ! number density scale parameter [particles/cc]
  rScale          =       0d0,              ! density scale [g/cc]
  TempScale       =       1d0,              ! temperature scale parameter [Kelvin]
  pScale          =       0d0,              ! pressure scale [dynes/cm^2]
  lScale          =       1.4959800E+015,   ! length scale parameter [cm] (AU=1.49598e13, pc=3.08568025e18, R_sun=6.955e10
  
  !================================================================================
  ! MHD related section
  !================================================================================
  lMHD             =  .false.      ! Magnetic Fields present? [.false.]
  lCheckDivergence =  .false.      ! Turn on divergence checking [.false.]
  /
  

• ideal gas solver, gamma=1.1
• temp_disk=1K; temp_amb=1000K
• dens contrast=10³
• r_soft=4= .5 disk_height
• r_disk=2 cu (=20AU)
• 64³ +2amr
• periodic BC

The gas evolution in this simulation is as expected and consistent with my previous tests (see tables 1 and 2 in ​https://clover.pas.rochester.edu/trac/astrobear/blog/martinhe12092011).

So… For folks interested in visualizing data remotely, I recommend installing visit on their local machine and using a visit data server on the host… I also recommend using the same version of visit on your local machine as on the remote host machine.

Here is a sample configuration for alfalfa which you can add under Options→Host Profiles

Check out updates in ​https://clover.pas.rochester.edu/trac/astrobear/blog/martinhe01202012

• Cold, temp=10K
• r_soft=0
• flared
• dens contrast=10³ (less than in tests 2, 2.1 and 3^*),
• High time resolution
• 64³ +2amr
• the disk has not been rotated
• r_disk=2 cu (=20AU; r_disk=1cu in test 2 and 2.1^*)
• periodic BC (note that the domain goes from -6 to 6cu (1cu=10AU), so the disk is till far from the boundary)

The central tilt happen again and quite quickly. There's no disk gas located at r < r_soft in this case. The knee shown by the dens distribution, i.e. the tilt, is also followed by the velocity distributions.

^*​https://clover.pas.rochester.edu/trac/astrobear/blog/martinhe01292012

This is the same as test 3 (​https://clover.pas.rochester.edu/trac/astrobear/blog/martinhe02032012) but the disk is initially inclined by 10^o with respect to the xy plane.

Central disk tilt is still evident, so seems that the grid coordinates, and split solving, are not introducing numerical artifacts which may cause the tilt.

Stratified ambient density:

LINEOUT KEY
Blue y=1
Red y=0.5
Black y=0
Green y=-0.5
Magenta y=-1

Clump

Rho color	50 yr	100yr	146yr
​http://www.pas.rochester.edu/~blin/feb2012/clump-strat.gif
Vx lineout	50 yr	100yr	146yr
​http://www.pas.rochester.edu/~blin/feb2012/clump-strat-vx.gif
Num dens -cm3	50 yr	100yr	146yr
​http://www.pas.rochester.edu/~blin/feb2012/clump-strat-numd.gif

Jet

Rho color	50 yr	100yr	146yr
​http://www.pas.rochester.edu/~blin/feb2012/jet-strat.gif
Vx lineout	50 yr	100yr	146yr
​http://www.pas.rochester.edu/~blin/feb2012/jet-strat-vx.gif
Num dens -cm3	50 yr	100yr	146yr
​http://www.pas.rochester.edu/~blin/feb2012/jet-strat-numd.gif

I've been trying to run the binary wind capture simulation using the reservation that we have in bluegene. I have however found the bug in ticket 121 (​https://clover.pas.rochester.edu/trac/astrobear/ticket/121). i.e., the code runs for a short while, aborts and reports:

hyperbolic/sweep/i_dependencies.f90", line 548: 1525-108 Error encountered while attempting to allocate a data object. The program will stop.

I do not see the error in bluehive (~ 32 afrank procs). I cannot go any further in bluegene with the following parameters:

• built at gbene, /home/mhuartee/27nov/
• revision 696
• ran at /scratch/mhuartee/astrobear2/feb-reserva/c, where all the relevant files (Makefile.inc, submit.cmd, *.data, binary ad chombo) are located and have open permissions
• job ID 36691
• 256procs
• produced one 25MB chombo00000 only (because I had set it to do so)
• the amr was off via: desieredFillRations=0 0 …; qTolerances=1e30 1e30…; refinemetVariableFactors=0 0 …
• 5 levels of sink particle refinement.

$> cat core.0 | addr2line -e ./astrobear

shows:

…

source/pointgravity.f90:352

source/pointgravity.f90:349

hyperbolic/sweep/sweep_scheme.cpp.f90:3571

hyperbolic/sweep/sweep_scheme.cpp.f90:3571

modules/objects/disks.f90:290

modules/objects/disks.f90:312

modules/objects/disks.f90:82

modules/objects/disks.f90:73

modules/objects/ambients.f90:128

modules/objects/outflows.f90:593

modules/objects/outflows.f90:279

modules/objects/outflows.f90:510

particle/particle_declarations.f90:568

modules/objects/outflows.f90:289

data/data_info_ops.f90:2513

data/data_info_ops.f90:2418

data/data_info_ops.f90:1961

data/data_info_ops.f90:1961

data/data_info_ops.f90:2024

data/data_info_ops.f90:2106

data/data_info_ops.f90:2091

…

I do not trust this error message since I've seen the error with a setup which had no sink particles, amr and much larger grids and chombos.

• More time resolution than tests 2^* and 2.1^*
• particle refinement, which is faster than the amr
• again^*, dens contrast=10⁶, temp=1000K
• the disk has been rotated by theta=pi/2.
• r_disk=2 cu (=20AU; r_disk=1cu in test 2 and 2.1^*)
• disk height=2r_soft (it was .75r_soft in test2 & 2.1^*). This increase of disk height seems to reduced the length of the central warped region in comparison with tests 2 & 2.1^*.
• periodic BC (note that the domain goes from -6 to 6cu (1cu=10AU), so the disk is till far from the boundary)

The vertical expansion has a velocity ~11.4km/s, which seems reasonable.

Density movie: ​http://www.pas.rochester.edu/~martinhe/2011/binary/3feb12-955.gif

Poloidal velocity movie [mach] ​http://www.pas.rochester.edu/~martinhe/2011/binary/3feb12-1050.gif

The poloidal velocity movie shows the gas located at r ~ < rsoft is quickly pulled towards the particle. Note this is not the ambient gas. This forms a thin -2dx wide-, small central region with colliding Mach 14 flows.
Shortly after, asymmetries, likely related to grid resolution and/or numerical diffusion, develop and produce shear. The effect grows in time. This seems to cause the early central warping. It doesn't make sense to me that this central region doesn't precess; I'd expect it to rotate about the disk ang mom axis, even if at a slower speed than the Keplerian one at r >~ rsoft.		.

The fact that we see central non-precessing warping also for a disk that rotates about the x-axis (and not about the z-axis as in previous tests^*), suggests there are no grid geometry bugs. Jonathan has found the same in his own, different, exploratory central disk tests. Still I'm now running a test with an inclination angle of 10^o, as suggested by Jonathan.

It's premature to say whether this instability will affect disk material located at r >> r_soft.

^*​https://clover.pas.rochester.edu/trac/astrobear/blog/martinhe01292012

The goal is to set up a hydrostatic isothermal rotating disk. Need to solve

\frac{\partial \rho}{\partial z}c_s^2 = \rho f_{z}

as well as

\frac{\partial \rho}{\partial s} c_s^2 = \rho f_{s} + \rho s \omega^2

for

\rho(s,z) and \omega(s,z) where f is the gravitational acceleration due to a point mass at the origin

First we can simplify the equations

\frac{\partial \ln \rho}{\partial z}c_s^2 = f_{z}

as well as

\frac{\partial \ln \rho}{\partial s} c_s^2 = f_{s} + s \omega^2

Since omega is a free parameter, we can choose a density field that satisfies the first equation and then integrate the second equation to find \omega as long as \omega^2 > 0

Let's try \rho(s,z) = \rho(s,0)\exp{\left( \displaystyle{\int_{z'=0}^z{\frac{f_z}{c_s^2}dz'}}\right )}

where \rho(s,0) describes the density along the midplane. We can set the height of the disk to be where the density drops to the ambient value - or we can set the background density distribution and replace it with the disk to the height where the pressures balance.

Also, if we don't use softening, then [[latex($f_{z}=-\frac{GMz}{r^3}

and \rho(s,z) = \rho(s,0)\exp{\left(\left [\frac{GM}{rc_s^2} - \frac{GM}{sc_s^2} \right ] \right )}

We can then solve the second equation for \omega^2=\frac{\partial \ln \rho}{\partial s} c_s^2 - f_{s}

If we now consider the density along the midplane \rho(s,0)

Since we expect f_{s} = -\frac{GM}{s^2} and since \omega^2 > 0 we have

\frac{\partial \ln \rho}{\partial s} c_s^2 > f_{s}

or

\frac{\partial \ln \rho}{\partial s} > -\frac{GM}{c_s^2s^2}

or \rho(s,0) = \rho_0\exp{\left(\frac{GM}{s c_s^2}\right)}

Of course this blows up at the origin, so we can either avoid the origin or use some form of softening.

Centrifugal acceleration can only support things from going inward. To keep things from going outward we have to use pressure or gravitational force - Since we want the density to decrease as we go outward, we have to use gravitational force - but this limits how quickly the density can drop. Lets choose

\rho(s,0)=\rho_0\exp{\left(\frac{(1-\alpha) GM}{s c_s^2}\right)}

where \alpha >= 0. If \alpha = 0 then we expect the equivalent of a BE sphere (but with a point particle instead of self-gravity).. Solving for \omega along the midplane…

s\omega^2 = -\frac{(1-\alpha) GM}{ s^2}+\frac{GM}{s^2} = \frac{\alpha GM}{s^2}

Combining all of this we get a solution for the density field

\rho(s,z)=\rho_0\exp{\left(\frac{GM}{rc_s^2} - \frac{\alpha GM}{s c_s^2}\right)}

and rotational velocity goes like

\omega^2=\frac{\alpha GM}{s^3}

Now for softening

First lets consider plummer softening. We can repeat the above procedure but instead of f_{grav}=-\frac{GM}{r^2}\hat{r} we have f_{grav}=-\frac{GM}{r^2+a^2}\hat{r} where a is a softening radius

Let's again try \rho(s,z) = \rho(s,0)\exp{\left( \displaystyle{\int_{z'=0}^z{\frac{f_z}{c_s^2}dz'}}\right )}

Now the integrand is -\frac{GM}{cs^2}\frac{z}{s^2+a^2+z^2}^{3/2} which is very similar to before except that r \rightarrow r_{eff} = \sqrt{r^2+a^2} and s \rightarrow s_{eff} = \sqrt{s^2+a^2}

so \rho(s,z) = \rho(s,0)\exp{\left(\left [\frac{GM}{r_{eff}c_s^2} - \frac{GM}{s_{eff}c_s^2} \right ] \right )}

And the integral along the midplane is also modified…

\rho(s,0) = \rho_0\exp{\left(\frac{(1-\alpha)GM}{s_{eff}c_s^2}\right)}

Combining these gives the solution for density

\rho(s,z)=\rho_0\exp{\left(\frac{GM}{r_{eff}c_s^2} - \frac{\alpha GM}{s_{eff} c_s^2}\right)}

and angular velocity

\omega=\frac{\alpha GM}{s_{eff}^3}

Spline Softening

If the softening function is compact (only modifies gravity inside of a certain radius) then

\rho(s,0)=\rho_0\exp{\left(\frac{(1-\alpha) GM}{r_{soft} c_s^2}+\displaystyle \int_{s'=r_{soft}}^s{f_s(s')ds'}\right)}

So we need the integral of the softened function…