Version 58 (modified by 12 years ago) ( diff ) | ,
---|
Z Cooling
The Physics of Z Cooling
Z cooling is an energy loss due to forbidden line emission. The routines use Pat Hartigan's cooling table Phcooling.tab. The table does not include recombination or permitted line emission.
Implementation into AstroBEAR 2.0
PHcooling.tab
Understanding the format of this table is crucial to interpolating correctly. There are a couple of properties that this table must have:
1.) It must be ordered as log(ne), T, log(X), (9 species cooling fractions), cooling rate. Where ne is the electron density, T is temperature, X is ionization fraction, and the 9 species cooling fractions are the fractions of how much each species contributes to the cooling rate. Note that here log refers to the base 10 logarithm. The 9 species are ordered as follows: OI, OII, NI, NII, SII, FeII, SIII, MgII, CII.
2.) The step sizes of log(ne), T, and log(X) must remain constant. In other words, each column is evenly spaced.
Also, it is important to note that the cooling rate has units of erg * cm-3 * nH-2. So before the rate can be added to the change in energy, dqdt(iE), it must be multiplied by nH2 where nH is the number density of hydrogen.
Zcooling.tab
PHcooling.tab contains a lot of data that is not needed. Reading in this data one line at a time is inefficient, especially when astrobear is being run on multiple processors. This is why a separate script called ph2zformat.s* is used to reformat PHcooling.tab into Zcooling.tab before running astrobear. Remember that PHcooling.tab still requires the aforementioned properties in order to be reformatted correctly. Zcooling.tab has the following format:
nDensities nTemps nXs lognemin lognemax tempmin tempmax logxmin logxmax coolingrate(1,1,1) coolingrate(2,1,1) ... coolingrate(nDensities,nTemps,nxs)
Where nDensities, nTemps, and nXs are the number of unique values for log(ne), T, and log(X) respectively. So if these parameters are equal to 40, 30, and 21, then there are a total of 25,200 cooling rates (40*30*21 = 25200). lognemin is the minimum value of the log(ne) column in PHcooling.tab, and the other minimum and maximum values are defined in the same way. Now the data does not need to be read in one line at a time for 25,200 or so values. The entire array of cooling rates can be read in with a single READ statement. No iteration or nested DO loops required.
*The only thing that ph2zformat.s assumes is that it is being run from the TABLES directory, and that the file you want to reformat is named PHcooling.tab.
Preparing for the Interpolation
There are several steps that need to occur before the actual interpolation. When Zcooling.tab is read in, the cooling rates are assigned to a properly allocated 3D array called Zcoolingtab(i,j,k). The indices i, j, k of Zcoolingtab are integers ranging from 1:nDensities, 1:nTemps, and 1:nXs respectively. These indices can also be thought of as coordinates. Before the cooling rate can be calculated, values for ne, T, and X are calculated. The electron density ne and ionization fraction X are converted to log(ne) and log(X). Then these values for log(ne), T, and log(X) can be scaled to the integer ranges of the coordinates i, j, and k. This is where the minimum and maximum values in Zcooling.tab are used.
In order to do a 3D interpolation, we need to find the 8 closest points to our newly scaled point of interest. You can think of these 8 points as the vertices of a cube or rectangular prism, and our point of interest is somewhere inside this box. First, we find the 8 closest coordinates. Since the coordinates i, j, k are all integers, the 8 closest coordinates are easily found by using the FLOOR and CEILING functions. Using these functions, you get a "low" and "high" value for each coordinate. Since there are 3 coordinates, 2*2*2 = 8 unique sets of coordinates.
The next step is to match these 8 closest coordinates with their corresponding cooling rates. This is done easily with 8 simple assignments. For example, rate1 = Zcoolingtab(lognelow, templow, logxlow) would be the rate that corresponds to the lower bound coordinates. Now we have the 8 closest points.
The final step before interpolating is to do another rescaling of the coordinates. We take the 8 closest coordinates to be the coordinates of a unit cube aligned with the origin. The coordinates of our point of interest are then scaled accordingly. This makes the interpolation very simple. Remember that in all these rescaling steps, nothing happens to the cooling rate because this is the value that we are interested in.
Trilinear Interpolation
Now we have 8 points that are positioned on the vertices of a unit cube aligned with the origin. Each point has a unique cooling rate. Our point of interest is inside this cube, but we do not yet know what its cooling rate should be. Trilinear interpolation is the method used to find this rate. The formula used is from Paul Bourke's website: http://paulbourke.net/miscellaneous/interpolation/. In cooling.f90, this is done on 8 lines to avoid having one lengthy, confusing line…
ZCoolingRate = rate1*(1d0-lognein)*(1d0-tempin)*(1d0-logxin) ZCoolingRate = ZCoolingRate + rate2*lognein*(1d0-tempin)*(1d0-logxin) ZCoolingRate = ZCoolingRate + rate3*(1d0-lognein)*tempin*(1d0-logxin) ZCoolingRate = ZCoolingRate + rate4*(1d0-lognein)*(1d0-tempin)*logxin ZCoolingRate = ZCoolingRate + rate5*lognein*(1d0-tempin)*logxin ZCoolingRate = ZCoolingRate + rate6*(1d0-lognein)*tempin*logxin ZCoolingRate = ZCoolingRate + rate7*lognein*tempin*(1d0-logxin) ZCoolingRate = ZCoolingRate + rate8*lognein*tempin*logxin
Where rate1,rate2,…,rate8 are assigned as aforementioned, and the "in" values represent the coordinates of our point of interest inside the unit cube.
Connection to NEQ Cooling
One component of the cooling rate calculated by NEQ cooling is called metal_loss. This is the energy loss due to metal excitation, and it is based on the ionization fraction and DM cooling curve. Z cooling passes the parameter iCooling to the subroutine Cool_Derivatives. Now, if iCooling==ZCool, metal_loss will be subtracted from the cooling rate, and it is going to be replaced by a different rate.
Back in the main cooling file cooling.f90, ZCoolingRate calculates the rate due to forbidden line emission. It is based on temperature, electron density, and ionization fraction. In the ZCoolingStrength subroutine, the metal_loss and ZCoolingRate are combined to give a weighted cooling rate. This is done by a parabolic weight function which I have called Zweight. The purpose of this weighted rate is to make the transition from NEQ cooling to Z cooling smooth.
The table for ZCoolingRate is currently from 2000 to 16,500 K. Regular NEQ cooling will only cool effectively to about 8000 K or so. So you need to use Z cooling if you need to look at cooling below about 1e4 K.
1D Jet Simulations
The set up for these simulations is very simple. The grid is initialized with an ambient object and a wind object. The goal is to study the effect of magnetic fields on the emission behind shocks in a pulsed jet.
Resolution Study
Resolution Study
Simulations involving cooling tend to require a lot of resolution. For this reason, I have been running a simplified version of the 1D Jet simulation to see how resolution affects the dynamics. Here are the initial parameters:
Ambient: density = 1000 particles/cc velocity = 0 km/s temperature = 1000 K By = 0 uG Wind: density = 10000 particles/cc velocity = 50 km/s temperature = 1000 K By = 0 uG Global Data: Domain length = 30 AU Final time = 1.7 computational units (approx 7 cooling times)
Note that these runs feature a stationary ambient, and they do not include MHD or jet perturbations. To get an idea of the form of the shock waves that are produced, here is a lineout of velocity:
This is from a high resolution run (100 base cells + 7 levels of AMR = 12800 effective). You can see that even with this resolution, the front shock is somewhat smooth instead of a sharp jump. Now in order to make some predictions based on the initial parameters, we need some equations…
Important Equations
A strong shock implies that the mach number M >> 1. We will also be assuming an ideal gas, so
. The following equations can be derived from the usual shock jump conditions:
Where s stands for shock, ps for post-shock, and a for ambient. T is temperature, v is velocity, n is number density, mH is the mass of a hydrogen atom, and k is Boltzmann's constant.
In order to figure out these post-shock values, we need to somehow relate them to the jet. This is where these equations come in:
Where j stands for jet, and
is basically the density contrast of the jet material to the ambient material. Now we can combine these equations and constants to come up with an equation for the post-shock temperature in terms of jet parameters:
Where vj is more conveniently in units of km/s. Now we need equations for the characteristic cooling time and cooling length:
Where
is the cooling rate in units of cm3erg/s. Using equations (3) and (4), these cooling equations can be rewritten in terms of our known initial parameters, and the post-shock temperature which we will want to have calculated already:
Some Predictions
Based on all these equations, we can make some predictions about the post-shock temperature and what resolution we might want. Using the aforementioned initial parameters, you get the following predictions:
Tps (103 K) | tcool (107 s) | lcool (AU) |
---|---|---|
58.3123 | 1.26155 | 3.20345 |
A rate of 4.7863x10-22 was used, which is fairly consistent with Dalgarno McCray. So if we have a domain that is 30 AU and we want 10 cells per cooling length, then we would need a resolution of about 94 cells. That doesn't seem too bad, so I will round up to 100 and start from there.
Results
Here we look at what post-shock temperature the simulation gives as a function of the effective resolution. The effective resolution = (# cells) x 2(# AMR levels). Each run used a base grid of 100 cells, so the only thing that changed from run to run was the number of AMR levels. For some reason, I could not get it to run with more than 7 levels, so the 25600 run is actually a base grid of 200 cells with 7 levels. For Tps I will use the highest temperature that I get from the simulation. The highest temperature always occurs early on in the simulation and then the post-shock temperature slowly decreases over time, but stays relatively constant at some lower temperature.
Z
Effective Resolution | cells/lcool | Tps (103 K) | Relative Error (%) |
---|---|---|---|
100 | 10.68 | 34.37 | 41.06 |
200 | 21.36 | 37.11 | 36.36 |
400 | 42.71 | 40.21 | 31.04 |
800 | 85.43 | 43.68 | 25.09 |
1600 | 170.85 | 46.66 | 19.98 |
3200 | 341.70 | 49.25 | 15.54 |
6400 | 683.40 | 51.18 | 12.23 |
12800 | 1366.81 | 53.20 | 8.77 |
25600 | 2733.61 | 53.66 | 7.98 |
And in graphical form, the table looks like this where I have plotted the expected Tps as the green horizontal line:
So the big question is, why does it appear that we need way more resolution than is predicted from the cooling length? Perhaps it has to do with how the refinement is being triggered. So I will try a resolution of 400, but with no AMR and see if it does any better than the 400 effective resolution which used 2 levels of AMR. With this run I got Tps = 40.59 x 103 K. So there really was not a significant difference from AMR to fixed grid.
My next idea was to revert back to just DM cooling instead of the new Z cooling. Perhaps, the Z cooling routine is not implemented quite right, and there is some "double counting" for cooling strength or something like this. So here is a data table where I used DM instead of Z. The improvement column is the decrease in relative error from the comparable Z cooling run.
DM
Effective Resolution | cells/lcool | Tps (103 K) | Relative Error (%) |
---|---|---|---|
100 | 10.68 | 25.71 | 55.91 |
200 | 21.36 | 31.63 | 45.76 |
400 | 42.71 | 41.07 | 29.57 |
800 | 85.43 | 47.65 | 18.28 |
1600 | 170.85 | 51.37 | 11.91 |
3200 | 341.70 | 53.30 | 8.60 |
6400 | 683.40 | 53.98 | 7.43 |
12800 | 1366.81 | 54.14 | 7.16 |
25600 | 2733.61 | 58.16 | 0.26 |
The DM cooling gets closer to what I expect but the Tps is still a little low. In both cases, Tps seems to be converging to some other value different from my expected value of 58.3123. However, when I do a run with no cooling, I get post-shock temperatures reaching 58.3123, and some even higher…closer to 60.97 with just 100 cells and no AMR. So I decided to take a closer look at the adiabatic case to see if I could figure out the discrepancy. Again, here is the data for no cooling:
None
Effective Resolution | cells/lcool | Tps (103 K) | Relative Error (%) |
---|---|---|---|
100 | 10.68 | 60.97 | 4.56 |
200 | 21.36 | 61.09 | 4.76 |
400 | 42.71 | 61.00 | 4.61 |
800 | 85.43 | 61.06 | 4.71 |
1600 | 170.85 | 61.76 | 5.91 |
3200 | 341.70 | 63.45 | 8.81 |
6400 | 683.40 | 64.17 | 10.04 |
12800 | 1366.81 | 65.20 | 11.81 |
25600 | 2733.61 | 64.16 | 10.02 |
You can see that these simulations actually give higher than expected post-shock temperatures. If I used these values as the "correct" values then the post-shock temperatures from the cooling simulations are even further off. Here is a plot that contains all the data:
Shock Structure
I've made some lineouts of density, velocity, pressure, and temperature to get a better understanding of how the simulation is interpreting the shock structure. Every image is taken at the middle of the simulation for the highest resolution run. The solid black line is density, dashed blue is velociity, dotted green is pressure, and dash/dot red is temperature. Everything is in computational units, and it's plotted on a log scale. The structure is more difficult to see in the cooling cases, so I included some magnified images for those two cases.
Attachments (9)
- vel_sample.png (23.8 KB ) - added by 13 years ago.
- temp_res.png (13.3 KB ) - added by 13 years ago.
- temp_res1.png (24.8 KB ) - added by 13 years ago.
- dm_lineouts.png (28.9 KB ) - added by 13 years ago.
- dm_lineouts_zoom.png (30.3 KB ) - added by 13 years ago.
- no_lineouts.png (28.4 KB ) - added by 13 years ago.
- z_lineouts.png (28.7 KB ) - added by 13 years ago.
- z_lineouts_zoom.png (30.0 KB ) - added by 13 years ago.
- nocool_noB_vel.gif (1.7 MB ) - added by 12 years ago.
Download all attachments as: .zip