I now have 4 different 2.5D pulsed jet runs, each at 2 different resolutions. The high resolution is at the same resolution that ​de Colle & Raga used. The low resolution is half of that. The low resolution runs seem less reliable in terms of cooling. The higher resolution runs show more structure, and the jet head evolution is different. This must be due to resolving the cooling length. In the image below (hydro runs), the low resolution run is on the left, and the high resolution run is on the right:

Since this resolution dependency also affects the emission, I have chosen to only do emission maps from the high resolution runs.

Emission Maps

Using the new 2.5D projections object, I was able to create emission maps. In theory, I could make maps for any emission line that I want, but I am only going to show Halpha and SII (6716 + 6731) for now. There are a couple of different ways to visualize the emission data, but I find the truecolor plots to be pretty neat. The scaling for the truecolor plots is done in a way that Jonathan described in a previous post (johannjc02202013).

Halpha is red, SII is green, and where they mix is yellow.

Beta	Inf (hydro)	5	1	0.4
Emission

Clump Emission Analysis

Figure 4 of de Colle & Raga shows a more quantitative analysis of how the magnetic field affects the Halpha emission for each of the 5 clumps within the jet. They plotted the ratio (MHD / Hydro) of the total Halpha emission for each clump. The result is pretty conclusive; as B increases, Halpha increases. I have tried constructing a similar plot with the low resolution runs, but the results don't follow a conclusive trend, so again I am only going to show the plot for the high resolution runs.

The total Halpha emission of a clump is found by summing the Halpha emission inside a circle of radius 1.4 centered at each clump. The clump positions are taken to be z = 6.2, 12.5, 19.37, 25.7, 32.53.

The numbers for the high resolution runs didn't make much sense either. Some clumps showed increase in Halpha with increase in B, while others showed decrease or little change in Halpha. I need to rethink how the Halpha and other emission lines are being calculated, or maybe come up with a better diagnostic.

Jonathan helped me alter the projections object to handle 2.5D data. The data is revolved around the symmetry axis (cylindrical z-axis = cartesian y-axis), and then summed along the line of sight (cartesian z-axis). This can be done for any field defined in fields.f90, but it is especially useful for creating emission maps.

How it Works in Astrobear

The projections object was also altered to handle an array of fields, instead of just one. Now you can have several projections in one .bov file, whereas before we would have a separate .bov file for each projection. Here's how you would initialize it:

CALL CreateProjection(Projection)
Projection%dim = CYLINDRICAL_PROJECTION
Projection%field(1)%id=Halpha_Field
Projection%field(1)%name="Halpha"
Projection%field(2)%id=SII_6716_Field
Projection%field(2)%name="SII_6716"
...

The CYLINDRICAL_PROJECTION parameter tells astrobear that this is 2.5D data, and it needs to revolve the data before integrating. Projections are always output to .bov files, but this is handled a little differently depending on how many fields your projection object has. If there is only one field, then the .bov file will contain the field name, and the variable in visit will appear as the field name. If there is more than one field, then the .bov file will contain the name 'projections', and the variable in visit will appear as a vector or array named 'projections'.

Visit

Visit can be a little quirky when handling a multi-component variable from a .bov file. I have seen that in visit 2.2.2, you cannot see each component separately. Instead, you have to write expressions to access each component individually. If the projection object has 2 or 3 fields, visit will read the variable projections as a vector. However, visit does not know how to use a 3-component vector within a 2D data set. So let's say for example that the 3 fields within your projection object are density, pressure, and temperature. In visit 2.2.2, you would have to write these expressions:

density = projections[0]     OR      density = array_decompose(projections,0)
pressure = projections[1]    OR      pressure = array_decompose(projections,1)
temperature = array_decompose(projections,2)  ! notice no choice for the 3rd component

If your projections object contains more than 3 fields, then visit will read the projections variable as an array, and then you have to use the array_decompose function for all elements of the projections array.

Jonathan found that visit 2.5 is a little better. This version of visit gives you the components of the projections vector or array individually. It names them consecutively as comp00, comp01, etc. There is still no way of naming the components by their appropriate field names before visit reads the .bov files. However, at least in this version of visit (and I would assume newer versions as well), you can rename the components with expressions without having to worry about that vector and array_decompose nonsense.

Thinking about the adiabatic BE sphere, I went through some calculations to check if we should change anything in the BE problem module. To begin, I refreshed my memory on deriving the Lane Emden equation, and then considered how the algebra would differ in an adiabatic case.

Isothermal case -

If we start with an ideal EOS,

and assume isothermality so we can use the isothermal sound speed,

we can rewrite the ideal gas equation to get,

Plugging this into the equation for HSE,

and integrating, gives

which plugging into

gives

Letting

allows one to scale the ODE into the following form called the Lane Emden equation:

Adiabatic case -

So, how does this change if we wish to use an adiabatic sound speed Cs, where

The ideal gas eqn when combined with the adiabatic sound speed is:

Plugging into HSE gives

implying,

This gives ~ the same ODE as above, but in order to scale it this time so that we get the same form as the Lane Emden equation above, one must make following modified variable substitutions:

Note the additional factor of gamma coming into xi.

Now the question is, how does this affect the BE module?

Well, the BE module reads in rhoc, r, xi, from which it computes an isothermal sound speed, aT, to get the temperature to set the pressure inside the clump.

The module uses the definition of xi above for the isothermal case to get aT:

NOTE that for the adiabatic case, the sound speed should have a factor of gamma coming in the numerator under the square root.

Now, considering the adiabatic sound speed,

gives

(The gamma factor cancels out anyway, leaving T unaffected as we should expect!!).

Results

1. T should be the same in an adiabatic BE sphere cloud (at least initially).
2. Xi, r, and rho_c are given in problem.data and together they fix the sound speed. Thus, the sound speed should be effectively different between the isothermal and adiabatic cases. However, as shown above, T is independent of gamma. Since T is used to set P, gamma doesn't effect the problem set up. Therefore, no changes should be made to the problem module.
3. Further, rho= rho(xi). Since xi has gamma in it for the adiabatic case, the stability criteria for the adiabatic case is going to be different. It would no longer be the case that for an adiabatic BE sphere, xi_crit = 6.451. Considering dp/drho shows that a factor of gamma would come into the P given by the adiabatic ideal gas law.

May have gotten VisIt plotting to work (seems to be a Mac - Linux problem that may disappear when script is run by server).

Working on writing a server - if that DOES work with the VisIt script, good. If not, then I will try alternate plotting.

Wrote Joint PhD proposal and met with Wendi Heinzelman.

TSF test run with Boss' BE sphere: When the incoming shock has Mach 1, the ram pressure from the shock matches the thermal pressure on the cloud edge since pressure remains the same across the isothermal shock. We get a slow blow of materials onto the sphere. The cloud density evolution in simulation:
​http://www.pas.rochester.edu/~shuleli/tsf/tsf02.gif

When the incoming shock has Mach 3.16, the ram pressure is about 100 times the thermal pressure on the cloud edge. We get a shock speed of about 5.78km/s. This shock should theoretically trigger the cloud collapse as indicated in Boss2010 but there should be very little injection. The cloud density evolution:
​http://www.pas.rochester.edu/~shuleli/tsf/tsf01.gif

When the incoming shock has Mach 20, the ram pressure is about 1.6e+5 times the thermal pressure on the cloud edge. In this case, we have a shock that is about 37km/s, falls into the triggered-injection window in Boss2010.
movie coming soon

I'll do the high res version of the above runs on stampede this week.

Gritschneder's setup: it seems the BE module only treats the isothermal case, so I need some input from Jonathan or Erica.

Papers: finished writing the MFU-IV draft.

The past couple of weeks have been completely devoted to fixing non-equilibrium cooling, 1D radiative shocks, and 2.5D pulsed jets.

Jonathan has done an awesome job of breaking down the DM cooling curve into its components. Now we can more accurately use equilibrium values to reproduce the DM curve, or we can use non-equilibrium values to get a more dynamic cooling curve due to ionization/recombination effects. In comparing with models done by Pat and his student Anna, it has become evident that our cooling routines probably have different physics. Or at the very least we are not putting the same importance or weight on various cooling processes. A closer look at the code that Pat is using will tell us more in the near future.

I was able to reproduce the pulsed jet simulations done by de Colle & Raga. The similarities in the plots look very good ehansen02212013. The emission maps need a little more work.

• astrobear alias - currently points to clover and we back up to botwin
  • Any reason not to move wiki server and repositories to botwin and backup elsewhere?
  • Need to find 'clover' references and rename to 'astrobear'
• Updated bov2jpeg program and documented it on the wiki
• Updated ImageObject and ProjectionObjects pages
• Implemented non equilibirum DM and improved cooling/chemistry/EOS module structure.
• Having trouble with Christina's run on BlueStreak Out of memory trying to allocate 482088528 bytes
  • Tried restarting with 2x processors but same error
  • Will try re-restarting to make sure nothing changed on system end
  • Will try restarting with 2x nodes but same number of cores (8 cores/node)
• This week:
  • Test new cooling routines
  • Implement HelmHoltz EOS
  • work on FLD

• New revision 1239:7aab6defde61 in Scrambler:
  1. working bov2jpeg from Jonathan.
  2. Currently has a problem with my bamboo account. So didn't run tests on bamboo.

• Tickets
  1. New: #276 (Update the testing suite scripts for HydroWaves IsoHydroWaves IsoMHDWaves MHDWaves on blue hive and blue streak)

• Weekly testing runs
  1. Results: Debugging

• New Hardware? (wireless mouse, Adapter Cable: Apple Mini DisplayPort)

• Worked on
  1. Clean up the users on wiki and astrolocals google group: #272 (Clean up forum/discussion board)
  2. Setup Stampede of Teragrid for AstroBEAR: #273(Install AstroBEAR on Stampede of TACC)
  3. #276

I was out of town last week so don't have too much to share, but since we last met here are some updates:

Science Outreach - I met with Prof. Luehmann at East High where she introduced me to Bertha, a sophomore in the Science Stars. I got to tutor her for her AP Physics class that she is struggling with. I scheduled to meet with her again in 2 weeks. Will also be planning a trip to the Planetarium for the group and will keep details posted.

AstroBear Movie - We have been given permission to use the astrobear comic image as our logo Would we like to crop this image down at all to exclude "AstroBEAR and Kid" (​http://nedroidcomics.livejournal.com/86373.html) ? In order to get the ball on the roll with making the PR video, can people send me or Baowei their schedules… What days/times are you unavailable during the week? I can forward these to our videographer to begin to schedule shooting times. Also, attached is the outline for the movies — we should decide on who will be interviewed for each section!

This is my last week working 100% with you guys. Starting March the 4th, I'll only be about 4 hours a week. Sometimes in the department but reachable via email.

Continuing projects: teragrid allocation, poloidal collimation, AGN with Eric Blackman and Alex Hubbard, probably the pulsars (see below).

QUESTIONS

• what office can I use?

• should I register for the APN6 conference?, ​http://www.astroscu.unam.mx/IA/index.php?option=com_content&view=article&id=664:apnvi&catid=44:undia-alias&lang=es

-Future projects?

Disk paper referee's answers. In Adam's hands now.

Poloidal collimation. Latest run almost finished. I'll then make a synthetic map of it. Waiting for Julien to give us an update on the paper writing.

AGN jet truncation, debugging the mhd module which showed problems with the latest code revisions.

pulsar winds+GRB, potential new collaboration. Enrico Ramirez-Ruiz (from Sta. Cruz), is an expert on GRB. He and Fabio de Colle want to try to do something like ​http://adsabs.harvard.edu/cgi-bin/bib_query?arXiv:astro-ph/0610454, and then to model GRBs after glows. I'll be in touch with them and Adam as this project progresses.

I'm rerunning my pulsed jet simulations, because the previous runs did not appear to be affected very much by the magnetic field. I"m using the exact set up used by de Colle and Raga, 2006. In the interest of time, my runs are at ½ the resolution of de Colle and Raga.

Beta	inf (hydro)	1 (weak field)	0.4 (strong field)
log(density)
emission

For the lineouts…solid is hydro, dashed is beta = 1 (weak field), dotted is beta = 0.4 (strong field). The head of the jet, which is to the left, is only included in the 1/beta plots.

Quantity [cu]	Image
density
ionization fraction
velocity (z-dir)
temperature
1/beta (r = 0)
1/beta (r = 0.6 Rjet)

I was thinking about Eddie's emission plots and how where there is h-alpha and sulfur 2 emission, the plot is brown instead of bright yellow - and this has to do with the way visit adds transparent pseudocolors… Basically it adds them as though they were pigments - not sources of light. Red+Green dye = brown, Red+Green light = yellow..

I thought about doing the colors in reverse - basically using anti-red (cyan) and anti-green (magenta) to get anti-yellow (blue), and then inverting the final image to get red, green and yellow - however, visit doesn't add transparent plots in a commutative fashion (the order of the semi-transparent images matters). Then I discovered visit's new true color plot which allows you to specify a rgb vector field and then plot the resulting image.

So if you create an expression for H-alpha emission and S-II emission that is scaled between 0 and 255,

H_sc=max(0, min(255, (H-Hmin)/(Hmax-Hmin)*256))) S_sc=max(0, min(255, (S-Smin)/(Smax-Smin)*256)))

then you can create a vector

{Halpha_sc, 0, SII_sc}

that you can then plot with a truecolor plot.

(I'm not sure why green is 3rd and not 2nd like 'rgb' - and I couldn't get blue to plot at all)

The following was produced by having the red value go from 0 to 255 across the left half and then stay at 255 on the right half, and the opposite for green. That way the center is at 255, 0, 255 = bright yellow….

TSF I set up a set of simulations with Boss cloud purturbed by an incoming wind:

• Pressure matched: no density jump, Mach = 1 (this is the case where the ram pressure from the wind equals the thermal pressure at the edge of the cloud)
  
• Density matched: the density of the wind matches that of the edge of the cloud, Mach = 1
  
• Boss condition: same as the density matched case but Mach = 10
  
• Magnetic Boss condition: same as the Boss condition but with a vertical magnetic field embedded in the domain (beta ~ 4)
  
  

Some of them should finish this week. I have also been going through some of the literature, the Gritschneder 2012 paper is interesting where they studied the triggered collapse of a marginally stable BE cloud with 10 solar masses and radius of 0.21, central mass about 10⁴ with DM cooling. The shock velocity is considerably higher than that of Boss 2010: ranging in 100km/s to 200km/s, which is far out of the velocity window proposed by Boss 2010 for successful injection and collapse(20~70 km/s). This raises the question of how does this "successful window" depend on the cloud radius and mass question, which I have not seen a complete set of simulation/modeling done yet. They proposed the observed abundance of Al26 due to shock injection, which is similar to Boss' work.
​http://arxiv.org/abs/1111.0012

The drawbacks are that these simulations are in 2D, no magnetic field or rotation, which they mentioned in the text:
" As this is an idealized case and the main focus of this study is the enrichment with SRLs, we don’t take any rotational or turbulent motion (e.g. Walch et al. 2010) inside the cloud core into account. While this initial rotation might be important in the later phases of the collapse - especially in determining the final disk size - its influence will be very small in the initial phase here, as the rotational velocities are orders of magnitude smaller than the shock speeds."
"For example, the core could be initially rotating. We assume here that the high shock velocities dominate over the rotation, but this should be further investigated. Furthermore, we neglect magnetic fields. They are of course present and should be involved in various processes, e.g. the details of the shock front and especially later on in the disk formation and the removement of angular momentum."

It could be easy to set up simulations using their setting with weak rotation (beta = 0.01~0.1) and weak~medium field strength (beta = 4~100) in a 3D simulation. One challenge is to replicate their resolution in 3D, which is effectively 512 zones per cloud radii, need a effective resolution of 1200*1200*6000-ish (self gravity in reflect boundary condition?).

Also been going through some of the rotating/magnetic BE cloud papers:
​http://arxiv.org/abs/0901.2127
​http://arxiv.org/abs/astro-ph/0408277
​http://adsabs.harvard.edu/abs/1993ApJ...417..220G
​http://adsabs.harvard.edu/abs/1993ApJ...417..243G

Papers I have been writing up the draft for the MFU-IV proceeding.

I finished the first round of pulsed jet simulations and generated some images. The pulsations are generated by a sinusoidal function with the average jet velocity being 60 km/s. The ambient density is at 100 cm^-3 and the jet is injected with a density of 1000 cm^-3. The three types of runs that I am showing differ in magnetic field strength. For the pseudocolor plots they are ordered hydro, beta=1, and beta=0.4 from left to right. The lineouts are straight down the axis of the jet, and for those we have the solid lines (hydro), dashed (beta=1), and dotted (beta=0.4). I should also clarify that the head of the jet is to the left for the lineouts.

Quantity	Image
log(density)
log(density)
Ionization Fraction
Velocity (z-direction) [km/s]
Temperature [K]
Emission (multiplied by cell area)
Emission (scaled to maximum)

• Weekly testing suite runs: Debugging
  1. Ran two time the 30 testing modules on blue hive. Each time 5 (different) modules problems failed due to the job queue system — the jobs died silently before running. They should all pass as it's the same revision as last week. Working on modifying the script to handle the case.

• Users:
  1. Created Wiki accounts for Andrew.
  2. Rui asked questions about hypre in AstroBEAR (Jonathan)

• Worked on ticket #273 (Install AstroBEAR on Stampede of TACC). Currently get error with the new installed hdf5 and qprec type.

Splitting DM into Hydrogen-Helium-Metal-Molecular processes

• So I spent a fair amount of time breaking down the curve in DM into 12 various parts that we use in non-equilibrium cooling. Here is a plot showing all of the different processes:

Now defining:

• x = HII/nH
• xHe = HeII/nHe
• yHe = HeIII/nHe
• y = ne/nH

each process can be identified with proper weight factors (before being multiplied by nH²)

And here is the reconstructed cooling function using equilibrium values from the code plotted alongside the total equilibrium curves extracted from the DM paper:

Other than the rise at 10⁴ due to helium excitation, the majority of the cooling comes from metals excitation by electrons (especially Oxygen at 10⁵). Of course the DM curve assumes some metal fractions (which are different from Pat's) so merging these two might be difficult. Though it would be possible to modify the DM curve to use Pat's metal fractions…

Breaking up the metal cooling

The form of the cooling function is dependent on the abundances of each element, and each of the sources used to construct the various tables use slightly different abundances:

DM Curve

Low Temperatures

The low temperature DM curve is dominated by carbon, silicon, and iron

The left plot shows the functional form for each of the lines considered from DM, while the right plot shows the total of all of the lines for each element along with data points extracted from a similar curve in the DM paper.

High Temperature

I also looked at how the abundance affects the cooling at high temperature. The DM curves are only for low temperature, although the completed curve uses cooling functions from Cox and Tucker. So, first I Combined the DM curves at low temperature, with the Cox and Tucker metal cooling curves (though using Gould and Thakur's hydrogen and helium excitation curves).

and then compared the resulting cooling function with total DM curve

The largest errors are around T=1e6 where the dm curve is weaker by a factor of ~2. Tucker and Koren give additional cooling functions for iron at high temperatures that I have not included…

Pat's Data

The data from Pat is a 3D-table PH_{cool}=n_H^2 PH(ne, \chi, T) while the DM metal curve is only a function of T DM_{cool}=n_H n_e DM(T)

So I divided Pat's table by \chi to get

PH'(n_e, \xi, T)=PH(n_e, \chi, T)/\chi

so that n_H n_e PH'=n_H^2 PH since Helium is not ionized in the range of Pat's tables…

I then wanted to construct curves at constant nH so I took slices through the 3D data set at n_e=\chi n_{Hi} where n_Hi = [1e2,1e3,1e4]. I then took lineouts at xi = [.01, .1, and 1] to construct 9 curves which I then plotted against the DM metal curve.

They all agree fairly well at high temperatures and deviate at low temperatures when the ionization fraction is high (which is significantly out of equilibrium). But I took an average of these 9 curves to get a single metal cooling curve from Pat's table.

I also looked at the DM cooling function at low temperatures but using Pat's abundances.

The cooling strength from Pat's tables is still much less than the DM curve. But the DM curve appears to have several iron transitions - and if all of them are included, it does not agree with Pat's tables. If I play around with the various iron and carbon contributions I can get a curve that agrees at low temperatures with Pat's tables.

I then combined that curve, with Pat's 1D average curve, and the original DM curve to construct a new metal cooling curve that is more consistent with Pat's cooling tables.

Comparing the cooling curves

This plot shows the DM cooling curves built out of Pat's abundances, the DM abundances, as well as the total metal curve extracted from the DM graph. Also are shown metal cooling functions from Pat's tables at ionizations of .01, .1, and 1 and densities of 1e2, 1e3, and 1e4, along with their average value, and a complete cooling curve that blends the DM graph with Pat's data.

And here is a close up of the range covering Pat's tables

And here is a fun image showing the shape of Pat's metal cooling functions for each species.

As well as the total cooling function

and a breakdown of which species dominates the cooling

I also did a breakdown of Pat's tables where I averaged in X and ne to get functions of T for each species

As well as taking the maximum value in X and ne

Since the iron cooling seems to stick out - I looked in more detail at the DM iron cooling processes:

Oxygen and Nitrogen cooling contributions

Since Oxygen and Nitrogen ionization states are tied to hydrogen via charge exchange, I looked at dividing out 'xi' dependence to make all of the cooling graphs 2D instead of 3D.

Here is a comparison of the original cooling functions for OI and OII with the reconstructed version built from a 2D table.

Also talked with Pat about extending tables to higher temperatures and what new species would need to be tracked.

About the attachments

See the attachments for excitation cooling of hydrogen and helium as well as metals. The DMCooling spreadsheet has all of the calculations, while the DMTables has the various temperature dependent functions for the original DM. And the DMZ Blend has the modified metal cooling function that is more consistent with ZCooling (Pat's tables). And the PHCooling.bov/dat files can be loaded into visit to see the function PH'(ne, x, T)

Hey all, my result for my Godunov code is on this page —

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

It's not immediately obvious why the solution is falling off beyond x~1.2. The initial condition is plotted in red… Up until that point, the solution matches Toro's. Not sure if I should take the time to diagnosis this. Probably should move on, got the big picture here of a) the Godunov method, b) marching through the time evolution, c) how to test/write/debug/etc fortran code..

Jonathan and I have met a few times now to discuss nonequilibrium cooling since he is implementing the molecular hydrogen routines and variable gamma. While trying to understand how NEQCooling modifies DMCooling, we realized that it was not being modified correctly. First, let me define some of the main parts in the different types of cooling.

Hex_eq = equilibrium hydrogen excitation
Hex_neq = nonequilibrium hydrogen excitation
DM_em = metal excitation via electron/metal collisions (from DM)
DM_Hm = metal excitation via hydrogen/metal collisions (from DM)
Z_em = metal excitation via electron/metal cololisions (from Pat)

Now, here is how the different types of cooling are related to each other:

DM = Hex_eq + DM_em + DM_Hm
NEQ = DM - Hex_eq + Hex_neq
Z = NEQ - DM_em + Z_em

NEQ

For the NEQ case, I think there was some confusion in the past and this was not handled correctly. Hydrogen excitation was being double counted in a weird way. There was a routine for an equilibrium ionization fraction table, which should have been used to get Hex_eq, but it was not used and I could not find a table. To fix this, I implemented a function to calculate the equilibrium ionization fraction:

! solves for equilibrium ionization fraction of hydrogen (derived from Saha equation)
FUNCTION EqIonFrac(Temp, nH)
   REAL(KIND=qPREC) :: EqIonFrac, Temp, nH, b

   b = (2d0*Pi*emass*Boltzmann*Temp/planck**2d0)**(1.5d0)*EXP(-13.6d0*ev/(Boltzmann*Temp))/nH
   IF (b>10d10) THEN
      EqIonFrac = 1d0
   ELSE
      EqIonFrac = (-b + SQRT(b**2d0 + 4d0*b))/2d0
   END IF

END FUNCTION

The equation for EqIonFrac is quadratic, and b is used with the quadratic formula to get EqIonFrac. The IF statement is in there because I learned a valuable lesson about limits and large numbers. You can see as T becomes large, the exponential in b goes to 1 and b scales as T^1.5. So b also becomes large as T becomes large. Due to rounding, the computer takes the square root part to be approximately b as b becomes large, and thus EqIonFrac to be zero. However, if you take the limit of EqIonFrac as b goes to infinity you find that it is 1 (which makes more sense physically). I found that 10¹⁰ is a good limit to use for b in that IF statement.

Z

The implementation of Z is slightly more complicated. Pat's cooling tables have a relatively limited range. Outside of this range, Z should revert back to NEQ. This transition needs to be smooth, so DM_em is not completely subtracted from NEQ and Z_em is not completely added to NEQ. These two components are weighted depending on the temperature.

Another issue is that in the code, we do not have a good way of separating DM_em and DM_Hm from the DM curve. So for Z cooling, Z_em is actually replacing DM_em + DM_Hm. This is okay as long as DM_Hm is small in the range of Pat's table. According to DM, DM_Hm is only important when the ionization fraction is low (approximately < 0.001). In Pat's table, the ionization fraction never goes below 0.01, so DM_Hm is not important and can safely be ignored.

Below is an example of how the fixes affected the temperature profile behind a particular shock. DM is black, NEQ is blue, Z is red…the old version is dashed, and the new version is solid. vs = 30 km/s, namb = 100 cm^-3, B = 0 (hydro run)

• Meeting with East High girls for outreach program on Friday morning

• Can't get a copy of Toro from printing center on campus

• Finished Godunov solver - this was for the 1D inviscid burger's equation. Just finishing up some tests of the code and will update my page on the results of this as soon as that's ready. (www.pas.rochester.edu/~erica/God.f90)

• Worked on writing results section — ~½ through (www.pas.rochester.edu/~erica/BEsphere.pdf)

• Next chapter in Toro is for Godunov method for systems of non-linear PDE's (Euler equations)

• Going out of town next week for some family things back in Philadelphia

Pulsed Jets

I altered my 2.5D jet module a bit and added forced refinement of the injection region (it's at 16 cells per jet radius). I completed 4 runs: beta = inf (hydro), 5, 1, 0.4…need to visualize the data now.

Rad Shocks

Spent almost all weekend practically rewriting my 1D rad shock module because I found that at longer run times the shocks were not steady. I rewrote the jump conditions in a different way, and was much more explicit about tracking the electron number density. I also started accounting for electron mass. For a long time I couldn't get it to work, but I just figured out this morning that it was due to the velocity scale. My density was being scaled much more accurately now, but in physics_control.f90 VelScale ~ SQRT(1/mu). In all other scale factors (except rScale), mu cancels out so I actually don't need to multiply the temperature by mu like I originally thought in cooling.f90. I still need to run some tests, but the shock looks rock steady now.

• Wiki
  1. Set yearly backup
  2. Updated Doxygen

• New Revision
  1. 1237:82b26a9a1a33 and 1238:604fb418ad9a checked in: just some updates with running scripts for blue streak.
  2. Weekly tests on blue hive and blue streak all passed: Added a naive weekly update for the testing to Debugging. So everyone can see the weekly testing result — The last week testing status won't be accurate for now.

• Users
  1. Contacted IO (Jonathan has the email) and Rui of LLE (the student Ariji is trying AstroBEAR now). Let them know the latest revision in developing branch.

I added a new EoS called EOS_MULTISPECIES_GAS which uses an array of species mass, free electron number, and gamma, to calculate pressure, temperature, etc…

And I've also started to look through the neqcooling routines to see how the code handles some chemistry stuff… Apparently the last version of AstroBEAR did not support any molecular hydrogen - just hydrogen and helium ionizations…

I tried looking back in the repository for molecular hydrogen routines and found the 1st revision checked into the repository had

 OIcool=OI_cool(nvec(iH2),nvec(iH),Temp)
 H2cool=H2_cool(nvec(iH2),nvec(iH),Temp)
 dmcool = DMCoolingRate(Temp)*ne*(nvec(iH)+nvec(iHII))

 H2diss = H2_diss(nvec(iH2), nvec(iH), nvec(iHe), ne, Temp)
 H2term = -H2diss + H2_dust_Recomb_table(nvec(iH)+nvec(iHII),nvec(iH),Temp)
 hion = H_ioniz_table(nvec(iH),ne,Temp) 
 hrec = H_recomb_table(nvec(iHII),ne,Temp)
 heion = He_ioniz_table(nvec(iHe),ne,Temp)
 herec = He_recomb_table(nvec(iHeII),ne,Temp)

 dqdt(iE)=-(OICool+H2Cool+dmcool+BindH2*H2diss+IonH*hion+IonHe*heion)

 dqdt(iH2)    = H2term
 dqdt(iHII)   = (hion-hrec)
 dqdt(iH)     = (-(hion-hrec) - 2.d0*H2term)
 dqdt(iHe)    = (herec-heion)
 dqdt(iHeII)  = -dqdt(iHe)

And there were two sections commented out. The first one had to do with a critical density for molecular H2 cooling

 CDH2=CriticalDenH2(nvec(iH),Temp)
 ! function used to smothly conect high HI and low HI density cooling rates
 CD=one/(one+CDH2)
 IF(CD<0.01) THEN
   H2cool=H2_critical_cool_table(nvec(iH2),Temp)*(one-CD)
 ELSE IF(CD<0.99) THEN
   H2cool=H2_cool(nvec(iH2),nvec(iH),Temp)*CD+H2_critical_cool_table(nvec(iH2),Temp)*(one-CD)
 ELSE
   H2cool=H2_cool(nvec(iH2),nvec(iH),Temp)
 END IF

and the other had to do with dust cooling

 dustcool=dust_cool_table(nHneuc-nvec(iH2),Temp)  ! dust cooling due to H sticking
 dqdt(iE)=-(OICool+H2Cool+dmcool+BindH2*H2diss+IonH*hion+IonHe*heion+dustcool)

Then in 2009 Bobby made a few changes

He added Helium double ionization though forgot to modify the source term for HeII to account for ionization / recombination

heIIion = HeII_ioniz_table(nvec(iHeII),ne,Temp) 
heIIIrec = HeIII_recomb_table(nvec(iHeIII),ne,Temp)
dqdt(iHeIII) = (heIIion-heIIIrec)
dqdt(iHeII)  = -dqdt(iHe)

He also (accidentally?) removed molecular hydrogen reactions

He then modified the dm cooling contribution from

DMCoolingRate(Temp)*ne*(nvec(iH)+nvec(iHII))

to

DmCoolingRate(Temp)*nvec(iHII) * (nvec(iH)+nvec(iHii))

which amounts to assuming that the number of free electrons is equal to the number of ionized hydrogen atoms… Which isn't true if there is helium around to ionize…

He also added hydrogen excitation cooling

Hex_cool=ne*nvec(iH)*HexCoolingRate(Temp)

and recombination energy losses

Rec_cool=(Boltzmann*Temp)*hrec + (Boltzmann*Temp)*herec

While the average energy of the recombining electron is 3/2 k_B T, I vaguely remember Pat saying something about integrating over the cross sections and the velocity distribution etc…. and that the energy loss is just k_B T

And finally he added BBC Cooling and changed the name of the dmcooling table file to DMcooling.tab instead of cooling.tab

Then 1 month later Jacob turned off HeII ionization and HeIII recombination completely which seems reasonable as it was only half-implemented…

And then Bobby implemented a Van der Waals equation of state…

Cooling Processes

1D Radiative Shocks

I have been trying to come up with a better way to figure out the cooling length and thus the appropriate length scale for these runs. While I was reading through the literature, I found a better way to calculate the immediate post-shock conditions, and I made some neat plots.

The points represent some of the actual runs that I am doing based on the grid of parameters that Pat sent me. Each line represents a different shock velocity. You can see that stronger shocks get closer to that maximum compression ratio of 4. I should note that this is the initial compression ratio across the shock jump; with cooling, the density compresses much more behind the shock. Also, a decrease in beta (stronger field) results in a smaller compression ratio. Furthermore, the field has no real noticeable affect until around beta = 10, and then a much greater affect as it gets close to and less than 1.

Similarly, the immediate post-shock temperature decreases with increasing field strength.

Figuring out the cooling length is much less straightforward. This is because it depends on the cooling function lambda. Even if I had a simple function for lambda (which I don't), the cooling length that you might calculate at the shock front is not what the actual cooling length will be. This is because lambda, as well as all other variables, are changing as you move farther behind the shock. The only solution I see is to figure out the cooling length numerically on a case by case basis with trial and error. However, that's essentially what my simulation does and it defeats the purpose of trying to predetermine the cooling length.

2.5D Pulsed Jets

I found a very useful paper that covers a lot of what I am trying to do (​de Colle & Raga 2006). There are several main differences that my research will explore:

• I will be using Pat's cooling tables which, as a reminder, include forbidden line emission.
• The differences between helical and toroidal field geometries.
• I have implemented the capability for pseudo-random jet velocities instead of sinusoidal pulses.
• Results will show other emission lines besides H-alpha.
• Eventually, moving to full 3D simulations might change the evolution.

I have a few runs of these jets with different beta. These runs are not pulsed yet, they have a constant jet velocity of 60 km/s, and they are characterized by their minimum beta.

A note on pressure profiles and beta

A lot of times, you'll see these jets defined by a plasma beta as: \frac{8 \pi P_{amb}}{B_{max}^2} However, if you look at the pressure profiles, you realize that this is not the minimum beta inside the jet. The minimum beta is actually: \frac{8 \pi \alpha P_{amb}}{B_{max}^2} This is the beta that I use when I characterize my jets in the above table. Alpha is a parameter that is defined by the radius at which B reaches its maximum (Rm = 0.6 for my runs). Alpha, Rm, and Bmax are not independent. You have to be careful with how you define these. For a given Rm, there is a maximum B you can use that will keep alpha and thus the thermal pressure positive. The equation for MHD equilibrium, the resulting pressure profiles, and other useful equations are given in several papers, including the one I already linked to above. Here is what the profiles look like when plotted:

Polloidal collimation.

Running the high res sim. Field lines seem to be fixed to the boundaries. However, I see spurious features which seem to develop from protections, see attachment (gray scale rim). I'm working no getting rid of this.
Found that the top part, and not the bottom —as we've been thinking—, of the lines outside the spherical wind BC (circular-ish gray region) move towards the axis, despite the fact that I'm using normal mag field BC. This needs to be fixed.

The above way of plotting the field lines is misleading. I've been playing around with the boundary conditions and with visit too.

1. killing the radial flow velocity in the bottom domain boundary introduces spurious flows.
2. the stream lines that we use to see mag fields in visit are misleading: I tried using different seed distributions in the grid to plot these lines. Some cases shows that the upper part of the lines moved towards the axis, hence giving the impression that the wind pushed the lines radially. Other seed configurations, on the other hand, showed that the bottom part of the lines were the ones moving.

So, I went back to use simple BC, where bz is forced to be parallel to both the bottom and top boundaries of the domain, and I do not touch the radial velocity at the bottom boundary. Here's a movie showing the results, descent resolutions, 2d cut (but the sim in 3d), and with log(number density) and magnetic field arrows. The fields seem to bend the way we want: stay frozen to the bottom and go around the wind. They only show a radial B component once the wind is well collimated (close to the end of the movie).

​http://www.pas.rochester.edu/~martinhe/andrea/6feb2013.gif

Done:

• Script written that
  1. Modifies module.
  2. Compiles.
  3. Runs.
  4. Plots with VisIt.
  5. Creates video of data.

Current Roadblock:

• VisIt doesn't work with ssh -Y connection - I don't need Xwindows to work, but VisIt assumes it does and crashes.

Next:

• Get VisIt to work via ssh.
• Begin implementing server interface for working with this framework.

• Rerunning 1D radiative shocks with a fixed length scale and adjusted plot ranges to send to Pat. Also, working on an analytic form for cooling length dependent on magnetic field strength. I'll use that to add the length scale as another parameter that changes for each run.

• For the 2.5D pulsed jets, I fixed the magnetic field and pressure profile following the form used by Lind in his 1989 paper. So now I'm doing some runs with varying maximum magnetic field strength. Below are images and movies of density, H-alpha, SII, and beta. This particular jet is not pulsed (v = 60 km/s). I changed some of the other parameters as well: namb = 100, njet = 1000, Tamb = 1000 K. Tjet is no longer set to a constant, it now follows the aforementioned pressure profile.

movie​

movie​

• Golden Version AstroBEAR & Blue streak
  1. New Revision in /clover/data/scrambler: 1235:f61e035a8ee7 including important bug fixes in cooling.f90 (#275) and clumps.f90
  2. Blue Streak: mercurial installed on blue streak. (Ticket #245)

• New publication with AstroBEAR: [AstroBearPublication]

• Tickets:
  1. New: #275: Segmentation fault found with RadShock testing module on Blue hive
  2. Closed: #234 and #235(both for Golden Version Test Modules), #245(Preparing Environment for AstroBEAR on Blue Streak (Q)), #256(GNU General Public License Licence for AstroBEAR), #275

• AstroBEAR Vedio Meeting with Erica, Brendan and Will

• Worked on:
  1. Working on building problem with MHD, following documentation page using magnetic fields: Field Advection, clumps with rotation and magnetic field: ​https://clover.pas.rochester.edu/trac/astrobear/blog/bliu01222013
  2. Updated documentation page: [ModulesOnAstroBear], [MultiPhysics] and [ClumpObjects]
  3. Debug and check in code.Ticket: #275

• Met with April Luehmann — see blog here for details: ​https://clover.pas.rochester.edu/trac/astrobear/blog/erica01302013

• Talked with guys at CRC about making astrobear movie. There are 2 short movies we are thinking of doing: 1. Public Release of astrobear2.0, and 2. Short intro for new users. The first will require actors. The second a narrator.

• Finished intro and methods for paper. Just have a few misc equations to add: ​https://clover.pas.rochester.edu/trac/astrobear/wiki/u/u/erica/BEPaper

• Magnetized Colliding Flow Update

movie​	movie​

• Looked at riemann solvers for general EOS johannjc01312013
• Van der Waals EOS implemented and tested johannjc02012013
• New equations of state possible - just need to provide functions (or tables) for:
  • Pressure given internal energy and density + tracers
  • Sound Speed given internal energy and density + tracers
  • Internal energy given pressure and density + tracers
  • Temperature given internal energy and density + tracers (technically not needed for euler equations)
• Does this automatically go into 'Golden Version' since collaborators are interested in using it? or is there a centralized development branch we should be pushing/pulling into?
• Organized existing documentation into Developer's Guide and various team pages.
• Updated the code to be doxygen compatible and enabled the doxygen tab
• Plans for this week
  • Write and implement Helmholtz EOS
  • Sesame tables? - maybe someone else wants to take this on?
  • FLD solver

So for any new equation of state, we need to be able at the very least to calculate:

If we want to know the actual temperature - we also will need a thermal EOS

and if we want to do any roe averaging we will need

As an example, I looked at a Van der Waals gas

So for a Van der Waals gas the pressure takes the form P(\rho,T) = \frac{\rho RT}{1-\rho b}-a\rho^2 However to solve the euler equations we need a caloric equation of state

So we need the internal energy density (which can't be derived from the pressure function alone) e=\frac{f}{2} \rho RT-a\rho^2 where f is the degrees of freedom.

This allows us to write our caloric EOS e(p,\rho)=\frac{f(p+a\rho^2)(1-\rho b)}{2}-a\rho^2

as well as solve for the pressure given the internal energy density p(e,\rho)=\frac{2(e+a\rho^2)}{f(1-\rho b)}-a\rho^2

so we can calculate the specific enthalpy

h(p,\rho)=(e+p)/\rho=\frac{f(p+a\rho^2)(1-\rho b)}{2\rho}-a\rho+\frac{p}{\rho}

and then use the specific enthalpy to calculate the sound speed using various thermodynamic identities…

c^2=-\frac{\left (\frac{\partial h}{\partial \rho}\right )_p}{\left ( \frac{\partial h}{\partial p} \right) _\rho - \frac{1}{\rho}}

c^2=-\frac{f\left [ 2a\rho^2(1-\rho b)-b\rho(p+a\rho^2) -(p+a\rho^2)(1-\rho b) \right ] - 2\left(p+a\rho^2 \right)}{f\rho(1-\rho b)}

c^2=\frac{f+2}{f(1-\rho b)}\left(\frac{p}{\rho}+a\rho\right) - 2 a \rho

which when a = 0 reverts to

c^2=\frac{f+2}{f(1-\rho b)}\frac{p}{\rho}=\gamma\frac{p}{\rho}

and when b = 0 reverts to

c^2=\frac{f+2}{f}\frac{p}{\rho} - \frac{f-2}{f}a \rho

which gives damping sound waves when

p < \frac{f-2}{f+2}a\rho^2

and when a=b=0 reverts to

c^2=\frac{f+2}{f}\frac{p}{\rho}=\gamma \frac{p}{\rho}

and finally we need to invert the equation for the specific enthalpy

p(h,\rho)=\frac{2\rho h + 2a\rho^2-f(a\rho^2)(1-\rho b)}{(2+f(1-\rho b))}

So we ran the basic shock tube tests comparing ideal gas and VanDerWaals gas.

Note that as \rho \rightarrow \frac{1}{b} then pressure should begin to dramatically increase with density due to limited volume.

And if \rho \rightarrow \frac{RT}{a} then attractive Van der Waals forces begin to significantly modify the pressure.

The parameters in cgs were

Note since \rho_L=10 \mbox{ so } \rho_L b = .3258 \approx 1

And \rho_L << \frac{RT}{a} so molecular attraction is negligible.

In what follows, the red line is Van der Waals, and the blue line is ideal gas.

If we match the temperatures in the left and right states we get

and if we match the pressures

and if we match the internal energy

The energy match looks identical to the temperature match because the " a " constant is so small (although the bounds are different on the plots).

There are 2 movies we talked about making today - 1) Public release video, and 2) instructional video.

I took down notes and can explain the general gist of things on Monday morning.