Posts for the month of June 2012

Martin's update, 7/2 '12

Ranger resources almost finished. Kraken ans Steele queues are insanely long.

Binary

  • 15 AU, forms a disk since .3 orbits. http://www.pas.rochester.edu/~martinhe/2012/binary/15au.png

clarifications on resolution for my 3D cooling jet runs

Here is an updated table from my previous blog post. Again, run A was not submitted to bluegene since it requires a lot of resolution. The required x res is the fixed resolution that I would need to guarantee 10 cells per cooling length.

run n jet eta chi required x res base grid levels actual cells per cooling length cells per jet radius
A 0.1 0.01 0.0015 78,335
B 1 0.1 0.034 3495 56 x 28 x 28 6 10.15 298.67
C 10 1 0.358 335 42 x 21 x 21 3 10.02 28
D 100 10 1.367 88 44 x 22 x 22 1 10.02 7.33
E 1000 100 2.428 50 50 x 25 x 25 0 10.12 4.17
C' 10 1 0.358 335 48 x 24 x 24 3 11.46 32
D' 100 10 1.367 88 48 x 24 x 24 3 43.74 32
E' 1000 100 2.428 50 48 x 24 x 24 3 77.70 32

I also figured out that if I want at least 32 cells per jet radius, then I need a minimum fixed x resolution of 348 which corresponds to a base grid of 48 x 24 x 24 with 3 levels of AMR. Unfortunately, when I set up these runs I was only considering cooling length to set the resolution, so runs C, D, and E are below this minimum resolution based on jet radius although C is very close.

The primed runs are new ones that meet the minimum resolution requirement. This essentially puts "extra" resolution on the cooling length. Fortunately, the old runs were still waiting in the queue, so I updated their global.data files to have the higher resolution. So C, D, and E were successfully replaced by C', D', and E'. I'll have to take a look at the results when I get back from vacation.

Meeting Update 0626

I have ran the low res version of the poloidal contained cases as discussed to determine the beta to be used in the "strong" cases. From the 4 simulations (lower res only the first 0.5 crushing times) of beta = 0.25, 0.5 aligned and perpendicular setups, it seems the beta = 0.25 is viable although you can observe obvious distortions happening at the beginning. Since we are more interested in the behavior after 1 crushing time, the initial distortion of beta = 0.25 should be no big deal. I then went on to run the cases needed for the paper on bluegene, with beta = 0.25 toroidal perpendicular, beta = 1 toroidal perpendicular, beta = 0.25 poloidal aligned and perpendicular. The first one is done (see result below), the second one is at 20%, the last two should be finished by Friday.
http://www.pas.rochester.edu/~shuleli/paperpics_clump1/Artical_File.pdf

I restructured the discussion part of the paper to have subsections with different beta. It seems too much with more than 10 3D figures and the coding is a bit confusing:

I still have queued runs of 3D clumps with more progressive beta values:
beta = 0.4, beta = 0.5, beta = 0.6.
Although we may not use these runs in the current paper, they maybe useful if we want to investigate the role of beta more thoroughly in the future (another paper probably): for the toroidal case, there seems to be a "threshold" beta around 0.6 that determines whether the clump collapse or open up, even in the perpendicular cases.

I've picked up the resistive clump runs from earlier this year. One problem we saw is that there are artificial effect at the head of the clump where the field is piling up in the B_y case. Currently looking into whether the subcycling is working correctly or not by comparing the result with a non-subcycled version.

Meeting update

I have updated my BE page with movies of the lineouts for rho and vrad of the matched case, where there seemed to be a bounce (see earlier blog for this reference as well), and also have added comparison plots to the earlier B&P collapse study. Both of these additions can be found here: https://clover.pas.rochester.edu/trac/astrobear/wiki/u/BonnorEbertMatched#Discussion.

I am also playing with getting the sphere to be in a bigger box by simulating just an octant of the original run and working through some calculations of Bondi accretion.

Baowei's Meeting Update 06/26/12

  • Worked on:
    1. Weak scaling test on Ranger: #193
    2. Hybrid scaling test on Kraken: #202
    3. #226
  • Will work on:
    1. Proposal for Teragrid allocation
    2. More scaling test on Kraken

Meeting Update 06/26/2012 - Eddie

I ran some very low res sims for the 3D cooling jets to make sure the problem worked and I had the dimensions right. The left image is what it looked like at frame 70 (100 frames total). I made the domain a bit bigger in all directions so that frame 100 now looks like the right image.

movie

Since these test sims looked okay, I submitted 4 jobs to bluegene this morning (Monday, 6/25). The wait in bluegene looks like it's about 3 days, so I guess I'll have to run the analysis when I get back from vacation.

With the ambient density held constant at 10 cm-3, I change the jet density to change the cooling strength chi. When chi is really small, this implies very strong cooling, and it becomes difficult to resolve the cooling region. All runs have a jet velocity of 250 km/s and a jet radius of 2e16 cm. Here are the varying parameters for 5 different runs:

n jet eta chi required x res
0.1 0.01 0.0015 78,335
1 0.1 0.034 3495
10 1 0.358 335
100 10 1.367 88
1000 100 2.428 50

The required resolution in the x-direction is based on a desired resolution of 10 cells per cooling length. You can see that the low chi runs require a lot of resolution, but there are ways around this. I can probably come up with a clever way to only track the jet head and region immediately behind it. For now, I did not submit a job for the first run in the table (njet = 0.1), and I threw some AMR levels at the other runs.

So here is my plan for when I get back:

  • analyze these runs that I have submitted
  • code something to force refinement only where I need it (this will help with low chi sims)
  • run some sims with Francisco's parameters (Adam, can you forward this info to me please?)

Movies of my lineouts

In the group meeting we noticed the density of the sphere of the matched case drop below the previous times when the density wave was moving into the sphere by 3.7 myrs (t=0.176 in computational units). We wondered if this was due to compression wave bouncing off of the dense core forming in the center. It turns out that a sink particle had not formed by 3.76 myrs… Also there does not seem to be a concurrent outward wave of velocity appearing that would seem to occur with a bounce…

Movies of the lineout in log space of density and radial velocity for the BE sphere in the matched condition:

https://clover.pas.rochester.edu/trac/astrobear/wiki/u/BonnorEbertMatched/LineoutMoviesMatched

Martin's update, 6/26 '12

Disks

15au and 30au cases running.

Note that the color scale on the right is smaller than that on the left to stress gradients.

Paper figures page: https://clover.pas.rochester.edu/trac/astrobear/blog/martinhe06202012

New cool movie from the 10 au case (the stupid arrow should be ~1 inch to the left): http://www.pas.rochester.edu/~martinhe/2012/binary/10au-3dDisk.gif

Co-rotating binaries disk formation, PAPER FIGURES

To do list form the meeting on June 21.

  1. Confirm we have a disk in the 20au case by comparing the disk's radial velocity profile vs a keplerian one.
  1. State in the paper that we are not far (far) from the WRLOF limit in the 10au (20au) cases, and refer to
  1. Assess the numerical viscosity in our sims by comparing current 10 and 20 au sims with different effective resolutions

10AU case. Not significant difference once steady. http://www.pas.rochester.edu/~martinhe/2012/binary/10fixedVSamr.png

  1. Analyze the accretion that we find with
    • profiles of radial mass flux on a longitudinal plane, vs polar angle and time
    • 2d maps of radial mass flux on the orbital plane vs time
    • compare profiles of Mdot vs time for 10au sims with
      • orbital motion (exiting data)
      • no orbital motion (new quick run)
  1. Runs
    • keep running current 20au
    • start running a 15 au
    • start running a 30 au

IN NO PARTICULAR ORDER

1. 10au accretion rate evo. Still higher than expected. val-borro et al. table3, report after 2 orbits for a 70AU case. The {\bf gradients} in the corresponding accretion rate evo profile (their fig. 12) are similar to ours; a very fast and brief initial increase followed by rather mild variations. http://www.pas.rochester.edu/~martinhe/2012/binary/10au-64x64x32-2amr-ranger.png
2. 20au accretion rate evo. http://www.pas.rochester.edu/~martinhe/2012/binary/20au-64x64x32-4amr-ranger.png
3. 10au Disk mass evo. http://www.pas.rochester.edu/~martinhe/2012/binary/10auDiskMass.png
4. 20au Disk mass evo. http://www.pas.rochester.edu/~martinhe/2012/binary/20au-diskMass.png
5. 10 au, disk density structure at 4 times, orbital plane view. The disk is consistently asymmetric. By time=.1 orb it has a width of about 1.5 AU. By time=1orb, the disk radius increases 4-fold and its morphology changes mildly. Density gradients are steeper in the disk part that faces the incoming wind. http://www.pas.rochester.edu/~martinhe/2012/binary/10au01.png
6. 10 au, disk density structure at 4 times, longitudinal plane view. The disk has a flared structure. http://www.pas.rochester.edu/~martinhe/2012/binary/10au02.png
7. 20 au, disk density structure at 4 times, orbital plane view. http://www.pas.rochester.edu/~martinhe/2012/binary/20au01.png
8. 20 au, disk density structure at 4 times, longitudinal plane view. http://www.pas.rochester.edu/~martinhe/2012/binary/20au02.png
9. 10 au (black) vs. 20 au (red), disk orbit streamlines comparison at 3 times: 1orb=thin; 2orb=thicker; 3orb=thickest, orbital plane view. Grid squares are 1 AU2. http://www.pas.rochester.edu/~martinhe/2012/binary/10vs20lines.png

updated documentation for my projects

My project pages were getting a little crowded, so I separated them and organized them in a logical manner on my page.

The work that I've done on 1D radiative shocks is now on a separate page. It has relevant equations and plots for the different types of cooling. The last plot on that page is an interesting one that I haven't shown before.

This is where I explain non-equilibrium cooling. I go into detail about what was changed from AstroBEAR 1.0 to 2.0. A lot of this information corresponds to ticket #147.

This Z cooling page is now mainly just about the implementation into astrobear. I also created ticket #224 for this, and it is part of the milestone called Implementing Multi-Physics.

This page contains remnants from my 1D pulsed jet simulations that I attempted last fall. I'm not sure if I'll be picking this project back up, or just moving on to something else for now.

This is the new page where I will post things for my 3D cooling jet simulations. I will also put the stuff about cooling length and cooling strength here, but the page is blank at the moment.

Meeting Update 0619

I'll talk about the clump paper related stuff.

http://www.pas.rochester.edu/~shuleli/clump0619/kinetic_selected_differentbeta.png
http://www.pas.rochester.edu/~shuleli/clump0619/thermal_selected_differentbeta.png
http://www.pas.rochester.edu/~shuleli/clump0619/BE_selected_differentbeta.png
http://www.pas.rochester.edu/~shuleli/clump0619/fieldcut.png
http://www.pas.rochester.edu/~shuleli/clump0619/beta25.png
http://www.pas.rochester.edu/~shuleli/clump0619/beta50.png
http://www.pas.rochester.edu/~shuleli/clump0619/beta75.png
http://www.pas.rochester.edu/~shuleli/clump0619/beta100.png
http://www.pas.rochester.edu/~shuleli/clump0619/totalpressureb1.png
http://www.pas.rochester.edu/~shuleli/clump0619/totalpressureb2.png

beta_avg=0.5
beta_avg=0.75
beta_avg=1.0

Meeting update

I've been studying the effects of a heavier medium on the collapse of a BE sphere. The sphere was found to be stable to collapse in a light medium (1/100 the density of the sphere's outer edge), but collapsed fairly quickly in the case where the density matched the sphere's outer edge. The infalling material increased the density, and hence pressure (sim. is isothermal so they are directly proportional) at the sphere's outer edge. Since the sphere was a critical sphere, any increase in Pext(Rbe), sends the sphere out of equilibrium. Thus it seems a rise in thermal P alone (not strictly Ram presure) was sufficient in inducing collapse.

My new project page can be found here:

https://clover.pas.rochester.edu/trac/astrobear/wiki/u/BonnorEbertMatched

Baowei's Meeting Update 06/19/12



Meeting Update 06/19/2012 - Eddie

  • some issues with Z cooling resolved, everything works now ehansen06082012.
  • figured out cooling length and cooling strength for 3D jet simulations which are based on the Blondin paper ehansen06122012.
  • testing Martin's jet module locally, will move to bluegene for initial set of runs
  • also started reviewing for the evil prelim

Martin's update, 6/19 '12

Binaries

This run produced sensible particle accretion rates: http://www.pas.rochester.edu/~martinhe/2011/binary/20au-64x64x32-4amr-ranger.png

Ranger is very friendly, e.g. consistent access to 64-512procs for 1 day with a queue wait time of <~ 30min. I suggest we apply for more time there soon.

-10au. New run with the same resolution as the above one, is running in ranger. Aim to see whether this 10au run produces lower particle accretion rates; I found higher values than the expected ones, for factors of ~30, for the 10au run in the table below.

-paper. The model/setup/IC sections are almost ready. I've been writing some of the 10au results. Missing: intro, more about the 10au results, the 20au results, some discussion and conclusions.

pre-run analysis for 3D cooling jets

Cooling Length

I've been following this Blondin paper to prepare for the 3D cooling jet simulations. http://adsabs.harvard.edu/abs/1990ApJ...360..370B. Section III is especially important because this is where the cooling length is defined. The analytic expression is there, but for analysis Blondin uses the numerical approximation from 1D simulations. I am therefore using my 1D radiative shock module to come up with my own formula for the cooling length.

Blondin's cooling length formula is going to be different because he did not use the DM curve. His cooling curve is based on a nonequilibrium ionization calculation by Kafatos (1973). Hartigan, Raymond, and Hartmann (1987) also have a different cooling curve and therefore a different cooling length formula. Blondin's is proportional to vs4 and Hartigan's is proportional to vs4.67.

The goal is to get a formula such that where dcool is the cooling length, na is the ambient number density, and vs is the shock velocity. b and a are the free parameters to be solved for. The formulas from Blondin and Hartigan do not match my data very well.

Source b (1016) a Sum of Error2
Blondin 4.5 4 125.96
Hartigan 1.80 4.67 121.39
Mine 3.51 3.2 10.42

Here, vs is in units of 100 km/s to give dcool in cm. Obviously, my formula is going to do the best because it's based on my data.

It is also important to note that the cooling lengths that I am calculating use a floor temperature of 8000 K just like Blondin. Hartigan's cooling length allows for cooling down to 1000 K.

Cooling Strength

Now, the cooling strength can be defined as: where rj is the jet radius, vj is the jet velocity, and is the density contrast (jet/ambient). The idea is to keep rj, vj, and na fixed. Then, you change the jet density and therefore changing in order to test different cooling strengths.

For some reason, I cannot reproduce the chi's that are reported in the Blondin paper. I've looked over my work several times, and I cannot find any errors. Maybe the paper has a calculational error or a typo somewhere. When I calculate chi for their "standard run" using their formula I get 0.27, but they claim that it should be 0.55. Not sure whether I should ignore the discrepancy and just do my own calculations, or if this needs to be resolved before moving on.

UPDATE

After doing more calculations for different sets of parameters, I realized that all of my cooling strengths were exactly ½ of the cooling strengths in the Blondin paper. It seems that Blondin either forgot a factor of 2 in their formula for chi, or all of their cooling strengths should really be ½ of what they reported.

Z cooling plots

Here are the other two important plots for comparing my shock module with Pat's model…

The overall shape/trend is good. Not sure why Pat's ionization fraction is a bit higher.

Density varies greatly towards the end of the cooling region. This could be because Pat's model keeps cooling past 2000 K where mine does not. Z cooling stops at 2000 K and NEQ takes over. The cooling rate for NEQ at those temperatures is very low.

Also, here is temperature plotted on a log scale. It's easier to see the differences at the lower temperatures this way.

New Revision 948:9cb35866cda0 in the main repository

I just pushed Revision 948:9cb35866cda0 in the main repository. It contains a modification to restarts where the master sends a message to each worker when it is ready to receive and process the data. WIthout this - the master can get overrun with messages on some platforms like Kraken.

A update list can be found at https://clover.pas.rochester.edu/trac/astrobear/log/Scrambler/?action=stop_on_copy&mode=stop_on_copy&rev=9cb35866cda06e8f69cce5ef787af702df0f1db9&stop_rev=946%3Aff6bdbea174a&limit=200

Details of modification to the code can be found at: https://clover.pas.rochester.edu/trac/astrobear/changeset?old_path=%2FScrambler&old=9cb35866cda06e8f69cce5ef787af702df0f1db9&new_path=%2FScrambler&new=946%3Aff6bdbea174a

Test results can be found at:

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

Gource rendering of scrambler...

http://youtu.be/grhEYDd0HnA

(The video shows all of the files in the scrambler repo - and who was working on which file when… - though drastically sped up :)

After a while the module directory begins to take on a life of it's own - and later on you can spot bear2fix - and the test modules… Also .f90 files are in green and data files in purple :) See if you can spot the module-objects directory

Don't forget to blast your favorite heavy metal tune in the background while watching…

And for those who reminisce about AstroBear 1.0

http://youtu.be/Q4l5B32glbI

ZCooling Works!

Thanks to Jonathan, ZCooling is finally working. He found the last tiny error in the cooling routines, and it made the shock steady. My shock is now even closer to resembling Pat's, and I believe I understand why it's different. The issues I was having with implementing a variable mu are fixed now too. So this is the full result…Zcooling, varying mu, MHD should work but these runs are not MHD. Basically, the difference between NEQ and Z cooling is that Z will give you much more cooling below about 10,000 K or so.

Pat's is in black, mine is red

There are two discrepancies here:

  • The big one is the difference between our immediate post-shock temperatures. I believe this can be accounted for by mu. Pat's model uses a bunch of heavier species, and I calculated his mu to be about 1.26 where mine is approximately 1.0. If my simulations had an initial mu of 1.26 my post-shock temperature would be much closer to Pat's. I should be able to resolve most of this by implementing He into my module.
  • My model has a bit of a jump near 2000 K. There is also a less noticeable one at 16,500 K. This is due to the way I am using Z cooling. NEQ cooling is a combination of a few things: hydrogen excitation, ionization, recombination, and metal excitation. Z cooling is supposed to replace the metal excitation component. Right now, Z cooling and NEQ metal excitation switch on or off at the upper and lower limits of the Z cooling table. This transition just needs to be "smoothed out".



UPDATE

I fixed the two discrepancies. I could not get He ionization/recombination to work, but just adding the right amount of neutral He increased mu and brought my solution much closer to Pat's. I also "smoothed out" my solution with what I call a parabolic weight function. I don't know if that's what the math world calls it, but it worked. I will play around a little more with the He ionization/recombination, but this is already looking good enough. I'll also post the density and ionization fraction plots soon. Again, my run is red and Pat's is black.

There is probably a good reason why my run seems to flatten out at 2000 K and Pat's keeps cooling a little more. The Z cooling table that I have stops at 2000 K, and so NEQ cooling takes over below that point. I suspect that whatever Pat uses is better at cooling at such extremely low temperatures.

Self-gravity thoughts...

Thoughts about the potential having negative curvature

3D

So in 3D - with radial symmetry - you can have a potential that radially has a negative curvature (ie ) and still satisfy poisson's equation since but rather

So the potential can go like r2 inside a uniform sphere and then switch to -1/r outside.

1D

For an infinite slab of uniform density - the potential inside goes like x2 and switches to just x1 outside… - so it has zero curvature outside as expected - and never has negative curvature anywhere.

Thoughts about the 1st derivative of the potential - or the Gravitational Field Strength |g|

Uniform density

Imagine a uniform density sphere in a very low density environment. As you move towards the sphere from infinity - the gravitational field increases because you are getting closer to all of the mass. Now as you enter the mass two competing effects happen. On the one hand you are getting closer to most of the mass, but at the same time you effectively stop feeling the rest of the surrounding mass.

Or in equation form…

Here the first term represents the fact that you are getting closer to the remaining interior mass - while the second term reflects the fact that the amount of matter inside is decreasing.

For a uniform density mass distribution, so we have

So the strengh of the gravitational field peaks at the outer edge of the cloud. This means that the outer material feels a stronger acceleration and falls in faster… This makes sense since if everything accelerated at the same rate - you would get material piling up at the center - and the uniform collapse is uniform. So the uniform collapse is essentially an outward in collapse

1/r mass distribution

If the mass distribution goes like 1/r something very interesting happens…

and we have

So every point feels the same inward acceleration. Matter will accumulate fastest at the center and you have an inside out collapse..

Point mass

The other extreme is to have a very centrally concentrated mass distribution (ie a point mass). In that case the second term is zero since as you get closer to a point mass - you don't see any less mass - so the strength of the field increases like - and so the acceleration is fastest near the point particle - and you would get an inside out collapse.

BE Sphere

A BE sphere is something in between a uniform density and a 1/r potential… And the inflection point will probably be where the density switches from a constant to 1/r - or it looks like at a non-dimensional radius of ~ 3… Of course we have completely neglected pressure gradients which are extremely important for a BE sphere…

BE sphere in a uniform density background

So now if we imagine a BE sphere surrounded by a uniform density ambient….

First the ambient effectively sees a point mass (provided that the BE sphere is much heavier than the ambient) - so it will want to collapse from the inside out. So a rarefaction should develop in the ambient. Except at the edge of the BE sphere - where the infalling material must pile up (or heat up and re-expand)… In any event this would add additional mass to the outer edge of the BE sphere and send a pressure wave through the cloud. The center of the BE sphere - must be completely ignorant of any additional mass piling up outside - so it will only see pressure waves. If the pressure wave is slow enough - or low enough amplitude - the central region should remain stable while the outside is free to undergo collapse…

Meeting Update 0606

Mostly working on NLUF proposal stuff and Hedla proceedings.
Current status:
http://www.pas.rochester.edu/~shuleli/NLUF/NLUF2.png

Just received Eric's comments on the proceedings, working on making changes to the manuscript. Will send the revised version once Adam get back his comments.
Also got some time to work on the conduction module. Now it looks like this on parallel fixed grid (works): AMR gives NaN fluxes, haven't had time to look at what's wrong.
http://www.pas.rochester.edu/~shuleli/NLUF/propagationtest.png

Meeting Update 06/06/2012 - Eddie

First, see my recent blog post ehansen06052012.

NEQ cooling with a varying mu is steady now, but the shock for Z cooling is still moving. I haven't been able to find the error yet, but I'm still working on it. Here's a temperature movie of both NEQ and Z. NEQ is dashed black line, Z is solid red line. NEQ_Z_temp.gif

I also started reading up on the jet simulations with cooling, so that I know what to do for these upcoming 3D runs. I have Martin's set up that he used for CRL618 and the Blondin paper http://adsabs.harvard.edu/abs/1990ApJ...360..370B.. Basically, the idea here is to run 3D jets with different cooling strengths. The cooling strength is defined in the aforementioned Blondin paper. For their different runs, they have some cooling curve (similar to DM) and they change the jet density in order to change the cooling strength. My question is…do I want to replicate the stuff in this paper by using their parameters, or should I start with the default parameters that Martin has used in his CRL618 runs?

Meeting Summary

  • Have good initial scaling results on Kure up to 256 cores #202… Moving to Kraken and 10,000 cores…
  • Worked on threading development for scaling results and various tickets…

Baowei's Meeting Update 06/06/12

  • Working on Teragrid allocation application information
  • We need thoughts about how to build a version of Astrobear working stably for most users and on most machines —our goal.
    Currently our new official revisions most come from bug-fixing for tickets from Martin and Jonathan while other group members may have their own revision working well for themselves. What would be a good way to merge these revisions together to get more features and better performance but not bugs? Also how to build a stable official version of astrobear? I'm preparing a document to collect the information about revisions/modules and machines each of our group member working with.

https://clover.pas.rochester.edu/trac/astrobear/wiki/ProjectStatistics

Bluegene/Q is coming in late June. I know there will be a gallery with big/huge TV monitor showing research photos/videos to the visitors. The university communication team is working on it. We don't have the details yet. But I think we should start thinking how to prepare for this gallery show.

Issues with NEQ and Z Cooling

I haven't been able to get steady shocks for these types of cooling yet. I originally thought that it had to do with mu (mean molecular weight). I need to calculate mu for every cell, because it is always changing. Heavier species will increase mu, and ionization will decrease mu. However, I just realized this morning that the shocks are not steady for a different reason. I went back to the beginning and just set mu = 1, and the shocks move. I'm guessing that there's a scaling issue somewhere since that has often been the case when I am close but not quite steady. There could also be a typo in the way temperature is being calculated. I'll look into this further and update this post later today.


Why mu is important

The following equation can be derived from the Rankine-Hugoniot conditions in the strong shock limit:

This shows that post-shock temperature is proportional to mu. Temperature in general is actually always proportional to mu because of the ideal gas law. We are just so used to having mu = 1.0 for an ideal hydrogen gas. For my pure hydrogen runs, mu is just 1.0, but Pat's runs have other species. The aforementioned post-shock temperature is immediately after the shock interface and before ionization occurs. Therefore, since the gas has no pre-ionization, the mu for Pat's runs is higher.

Let's look at the temperature for both runs (Pat's is in black, mine is in red):

You can get a rough estimate of what mu should be in Pat's runs by just looking at the difference in temperature. I suspect that Pat has an initial mu of approximately 1.22.

Here is how mu is calculated in the cooling routines:

 mu = (nvec(iH2)*muH2 + (nvec(iH)+nvec(iHII))*muH + (nvec(iHe)+nvec(iHeII)+nvec(iHeIII))*muHe)/npart

nvec is the number density vector, and it contains the number densities for all the different species. npart is the total number density of particles which includes all the electrons as well. This is why mu significantly decreases with ionization. For example, a fully ionized gas of HII would have mu = 0.5 instead of 1.0 for a neutral gas of all H. In other words, you double the number of particles without adding any mass (electrons are treated as mass-less). Recombination has the opposite effect as it reduces the number of electrons.

So in general, mu is important in calculating temperature. Both for the immediate post-shock temperature, and the decreasing temperature within the cooling region.

Here is what mu looks like throughout the problem domain:

Notice that mu starts off a little less than 1.0. That is because there is a trace amount of pre-shock ionization. Mu then drops rapidly (to about 0.65) after the shock due to ionization,and gradually increases due to recombination.

Martin's update, 6/6 '12

Disks

-20au preliminary early evo test with more resolution than before, http://www.pas.rochester.edu/~martinhe/2011/binary/20au-1000k-64x64x32-4amr-correctWINDdens.gif. This sim has a resolution of 64x64x32+3amr. AGB wind inflow enters the -x, +y and +-z boundaries. Improved versions of this sim should run in kraken soon.

  • 10au. Evo plots:
    • NEW http://www.pas.rochester.edu/~martinhe/2011/binary/6jun12.png
    • OLD: http://www.pas.rochester.edu/~martinhe/plot.pdf, these are the disk mass and the particle accretion rate, red and green respectively. Mdot (green) is higher than expected for an AGB mass-loss rate of 10-5Moyr-1, however, I found a mismatch between the wind's density at the position of the grid's inflow boundaries, relative to the primary's position. Thus the injected mass-loss is actually 3x10-4Moyr-1. With this in mind, a Mdotparticle ~ 4x10-6Moyr-1 (green profile) is sensible. I'm running a 10au and a 20au (with more res) versions, with a corrected mass-loss of 10-5Moyr-1. I can already see a lower Mdotparticle for the 10au case.
  • Here's the very preliminary paper's draft. A LOT is needed still, http://www.pas.rochester.edu/~martinhe/2011/binary/draft1.pdf
  • Progress in kraken (managed to run for ~8hrs on 4092p), but issues too, see tickets 209 and 217. Will try the latest rev.
  • The code is behaving well in steele, but queues are very, very busy.

Towers. Working on the ApJ referee's comments on the paper.

AGN jet truncation by red giant stellar winds. 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.