%\documentstyle[11pt]{article} 

%\input{rep12}
%\input{symbollga}
%\input{psfig}

%\begin{document}
%\setcounter{page}{42}
\newpage
 
\begin{center}
\Large{\bf V. TEST PROBLEMS}
\end{center}   
\addcontentsline{toc}{section}{\protect\numberline{V.}{\bf ~Test Problems}}
 
\vspace*{5mm}
{\flushleft\large{\bf A. Brief Description}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsection}{\protect\numberline{V.A.}{ ~Brief Description}}
 
\begin{center}
\it References
\end{center}
 
{\flushleft\rm
\hspace{10mm}G. A. Sod, \it J. Comput. Phys., \bf 27, \rm 1, 1978\\
\hspace{10mm}P. R. Woodward and P. Collela, \it J. Comput. Phys., \bf 54, \rm 115, 1983\\
\hspace{10mm}L. D. Landau and E. M. Lifshitz, ``Fluid Mechanics'', Pergamon\\
\hspace{15mm}Press, London, 1959\\
\hspace{10mm}L. Ensman, Ph. D. Thesis, University of California at Santa\\
\hspace{15mm}Cruz, 1991\\
\hspace{10mm}Ya. B. Zel'dovich \& Yu. P. Raizer, ``Physics of Shock Waves and\\
\hspace{15mm}High-Temperature Hydrodynamic Phenomena'', Academic Press,\\
\hspace{15mm}New York, 1966\\
\hspace{10mm}D. Mihalas and B. W. Mihalas, ``Foundations of Radiation\\ 
\hspace{15mm}Hydrodynamics'', Oxford University Press, Oxford, 1984}\\
\vspace*{2mm}
 
\noindent\rm
TITAN comes with a suite of test problems that were designed to illustrate
how versatile and powerful the code is. The test problems divide naturally
into three categories: hydrodynamical tests, radiation transfer tests in
static media, and radiation hydrodynamical tests. \vspace{2mm}

\noindent In the first class fall
Sod's shock tube in which a dense hot gas expands into a dilute cold one,
Woodward's problem, in which pressure jumps generate two interacting shocks,
and the Sedov-Taylor problem, in which an inital energy deposition creates
a self-similar spherical blast wave. \vspace{2mm}

\noindent In the second class fall radiative
heating problems, in which a light source heats a gas until it settles
into an equilibrium, and radiative cooling problems, in which the light
source of the previous situation is being turned off. The radiative problems
allow the solutions to be obtained through equilibrium and non-equilibrium
diffusion as well as through full transport with variable Eddington factors.\vspace{2mm}

\noindent In the third class fall the sub- and supercritical shocks, in which a moving
piston generates a non-adiabatic shock with radiative features, and a
radiative blast wave, in which an energy deposition leads to an explosion
with shock fronts and radiative waves. \vspace{2mm}
 
\noindent\rm
In TITAN the flag {\tt lprob} identifies a particular test problem and 
selects the appropriate start routine as well as the correct formats of
the output. This flag is assigned in the input file as described in 
section IV.B.\\
 
\vspace*{2mm}
{\flushleft\large{\sl B. Hydrodynamical Test Problems}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsection}{\protect\numberline{V.B.}{ ~Hydrodynamical Test Problems}}

\noindent For the hydrodynamical tests an ideal gas $\Pg = \left(\gamma-1\right)\rho\eg$
is considered to be initally at rest within a cylinder of unit length. The adiabatic
index is chosen to be $\gamma = 1.4$.
 
\vspace*{2mm}
{\flushleft\large{\sl B.1. Sod's Shock Tube}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsubsection}{\protect\numberline{V.B.1.}{\sl ~Sod's Shock Tube}}

%\noindent\rm (1) Sod's problem ({\tt lprob = 3}):\vspace{2mm}

\noindent
The first test problem is the classical problem in which a dense hot ideal gas 
occupying the left half of a cylinder penetrates into a colder dilute gas residing 
in its right half. It has already been shown by Riemann that a shock wave, a contact 
discontinuity, and a rarefaction wave are formed. A well-known analyses of 
corresponding numerical simulations has been carried out by Sod. In his setup the 
gas is initially separated into two parts by an eightfold density and a tenfold pressure 
jump. \\
\sodinit\\

\noindent
The computation is carried out with time centering parameter $\theta = 0.55$ on
100 grid points. An artificial mass and energy diffusion of $\sr = \se = 10^{-5}$ is used 
in order to treat the heating at the right wall from the shock reflection.\vspace*{2mm}
 
\noindent\rm The fixed grid (Eulerian) case is run with a linear artificial viscosity 
$\ell_0 = 10^{-3}$, \ie~one tenth of the grid spacing. Using typically four iterations per 
time step the time $t = 1$ was reached after 670 time steps. For the adaptive grid case 
a $\ell_0 = 10^{-4}$ was chosen. As the grid parameters $f_k^n~$'s were employed the density, 
pressure and energy in logarithmic scaling together with $\alpha = 1.5,~\tau = 10^{-15}$. 
TITAN needed about eight iterations per time step and reached the time $t = 1$ after only 
230 time steps. Since one additional equation had to be solved the actual CPU time for 
completion of the run was approximately the same in both the adaptive and fixed grid 
version.\vspace{2mm}

\noindent The first figure shows snap shots of the density profile on the adaptive grid at four 
different times $t\in (0.034, 0.283, 0.316, 0.996)$ from top left to bottom right.\\

\psfig{figure=sod1.ps}
\centerline{\small Fig. 4: Density profiles at different times ($+$'s indicate grid points)}
\vspace{2mm}

\noindent
The top left plot shows the formation of Riemann's solution: the right-most step being the
shock front, the next step indicating the contact discontinuity, followed by the rarefaction
fan. In the top right figure the shock rams into the wall while being followed by the 
discontinuity at a slower speed. The shock has reflected (bottom left plot) and moves 
leftward towards to discontinuity. In the last drawing the shock has moved beyond the middle
of the shock tube and drags the rarefaction wave behind it which connects it to the
contact discontinuity.\\

\noindent
One can obtain a complete overview of the actions of all nonlinear waves by drawing a 
space-time diagram of the density. This is done by stacking many of the its profiles and
plotting the contour lines of the result:\\

{\vbox{%
\psfig{figure=sod2.ps,width=4.7in}
\centerline{\small Fig. 5: Density contours}
}}
\vspace{2mm}

\noindent
For times $t<0.3$ there are three nonlinear waves: the right-traveling shock front and contact 
discontinuity and the left-traveling rarefaction waves. The first two can be identified
by their respective speeds which can be read off the slope of the right-going lines. The
shock speed is approximately +1.77 while the speed of the discontinuity is about +1.0. At
$t=0.283$ the shock reflects from the wall and travels backwards with a slower speed of
about -1.2. It collides with the contact discontinuity around $t=0.41$. Both waves interact and
penetrate each other. The shock and the contact discontinuity are slowed down considerably and a
second right-traveling shock emerges as a new wave. It reflects from the wall around $t=0.5$
and shortly thereafter interacts with the contact discontinuity. As a result the latter one
is slowed down even further to an almost halt but the computation does not tell what happens
further. The contact discontinuity and the first, left-propagating shock are connected by a 
rarefaction fan after their interaction. The primary left traveling rarefaction waves reflect
also about $t=0.5$, and when they reach the first shock, it starts to speed up.\vspace{2mm}

\noindent
The primary shock front and the contact discontinuity appear as sharp lines in the space-time
diagram above. This means that they are well resolved by many grid points which is possible
through the action of the adaptive grid. The following space-time diagram shows the motion
of every fourth grid point during the calculation:\\

{\vbox{%
\psfig{figure=sod3.ps,width=4.7in}
\centerline{\small Fig. 6: Action of the adaptive grid}
}}
\vspace{2mm}

\noindent
As soon as the nonlinear waves form the adaptive grid detects the two right- propagating
waves. This is possible through the choice of grid parameters: density, pressure, and
energy. The grid points zoom into the discontinuities on a very short time scale until
they are numerically resolved and track them. When the shock reflects at the right wall 
the adaptive grid tries to relax (because momentarily there is no shock) and begins to
move the grid points towards a more uniform (\ie~Eulerian) distribution. This is, however,
not done instantanteously because the mesh remembers the previous configuration via the 
grid time scale $\tau$. When the shock moves away from the wall the adaptive grid recognizes
it again and re-distributes the grid points so as to resolve it. The shock and the contact
discontinuities remain resolved and tracked until the run terminates. By contrast, the
rarefaction waves and the secondary shock are not detected because they do not show any 
strong variations in the grid parameters. This situation could be improved if they were
chosen differently.
 
\vspace*{2mm}
{\flushleft\large{\sl B.2. Woodward's Interacting Blast Waves}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsubsection}{\protect\numberline{V.B.2.}{\sl ~Woodward's Interacting Blast Waves}}

%\noindent\rm (2) Woodward's problem ({\tt lprob = 4}):\vspace{2mm}

\noindent
Woodward introduced a planar problem of two interacting blast waves. The cylinder is 
divided into three sections with three different pressure regimes.\\
\woodinit\\

\noindent
In the left part the gas pressure is $P(0\le r\le {1\over10}) = 1000$, in the middle part 
$P({1\over10}< r< {9\over10}) = 0$, and in the right part $P({9\over10}\le r\le 1) = 100$. 
The density is initially constant $\rho = 1$ across the domain $0\le r\le 1$. As grid 
parameters we adopted the density, pressure and energy in logarithmic scaling, together with 
$\alpha = 1.5,~\tau = 10^{-6}$. This time the code is run with 200 grid points. The artificial 
viscosity is kept linear at $\ell_0 = 10^{-4}$.\vspace{2mm}

{\hbox{%
{\vbox{\hsize=2.5in \noindent \strut
The large pressure steps at both ends of the cylinder create two shocks which travel 
towards the middle. The graph to the right (fig. 7) shows their density profiles at an
early stage ($+$'s marking individual grid points). One can see that the left front is 
much stronger than the right one. Further, rarefaction waves begin to appear in form 
of the density dips, which connect to the unit density near the walls.\strut}}
{\vbox{
\psfig{figure=wood0.ps,width=2.5in}
\hspace{12mm}{\small Fig. 7: Density profile}}} 
}}
\vspace{5mm}

\noindent
In the following figure, density profiles at four different times $t\in($1.57, 2.94, 3.85, 
5.00$)\times 10^{-2}$ are plotted. The top left graph shows the approach of the two shocks. 
Note in particular that the left one is much faster than the right one. The fast shock has 
a resolution of $10^5$ and the slow shock one of 
several $10^4$. The density between the shocks is still at the $\rho = 1$ plateau thus 
indicating that they move independently. The shocks are followed by contact discontinuities, 
both being resolved by $10^3$. They connect to the rarefaction fans which stretch to the walls.
The high resolution of the discontinuities leaves only a crude representation of the 
rarefaction fans with resolution of about $20$. It is by choice that the adaptive grid is 
allowed to zone only into the shocks and has to neglect the rarefaction waves. At time
$t=0.028$ the shocks collide near position 0.69 and penetrate each other. When they separate
a new third contact discontinuity forms in between them (lower right plot). At $t=0.05$
the fast shock has reflected from the right wall but has not yet run into the third 
contact discontinuity.\\

{\vbox{%
\psfig{figure=wood1.ps,width=4.7in}
\centerline{\small Fig. 8: Shock profiles at different times ($+$'s indicate grid points)}
}}
\vspace{5mm}

\noindent
A complete overview of the actions of all nonlinear waves can be obtained by investigating
the space-time diagram of the density. The respective contours in logarithmic scale are
drawn in figure 9. The shock fronts (due to their high resolutions) appear as thick lines, 
the contact discontinuities as broad lines, and the rarefaction waves as an ensemble of 
single contour lines.\\

{\vbox{%
\psfig{figure=wood2.ps,width=4.7in}
\centerline{\small Fig. 9: Density contours (logarithmic scale)}
}}
\vspace{2mm}

\noindent
The five density profiles (figures 7 \& 8) are horizontal cross sections at the respective
times. For times $t<0.3$ there are the fast left traveling and slow right traveling shocks
(speed given by the slope of the thick contour lines) each followed by a contact
discontinuity which connects to a rarefaction fan. The fans are reflecting soon from the
walls indicating that the region near the walls are being depleted quickly. At $t=0.028$
the two shocks collide and slow down after the interaction while creating a secondary contact
discontinuity. Shortly after the fast shock runs into the primary discontinuity from the
right and accelerates again to near its original speed while the discontinuity is reflected.
At $t=0.045$ this shock reflects at the right wall and interacts with the discontinuity at
about $t=0.05$. Both nonlinear waves are slowed down considerably. The situation is comparable
Sod's shock tube. New, however, is that the shock interacts with the secondary discontinuity 
and speeds up again. Shortly before $t=0.045$ the slower right propagating shock collides
with the other primary discontinuity and is greatly accelerated while the discontinuity is
reflected. This shock now propagates uninterruptedly to the left wall. There it reflects 
at time $t=0.063$ and moves with similar speed towards the right again until it runs again
into the discontinuity. This time, however, the other shock front meets also with them leading
to a multiple interaction of nonlinear waves.\vspace{2mm}

\noindent\rm
A comparison of TITAN's computation with the original one by Woodward and Collela (using a 
special version of the PPMLR scheme on a domain of 3096 zones) tells that the shocks are 
well represented but the contact discontinuities lack sharpness. This can be accounted for 
in part to a fairly large artificial mass and energy diffusion which allowed TITAN to 
complete the run in a little over 1000 time steps (about a factor 5 less than the number of 
time steps PPMLR required). For comparison reasons the space-time diagram of the grid motion 
is presented in the following figure.\\ 

\psfig{figure=wood3.ps}
\centerline{\small Fig. 10: Action of the adaptive grid}
\vspace{2mm}

\noindent
Figure~10 demonstrates well that the adaptive grid keeps track of all discontinuities. 
Note the strong correlation between the grid motion and the density contours.
Its typical behavior at the momentary absence of the shock(s) due to reflection from a wall
or shock-shock interaction is also evident.
 
\vspace*{2mm}
{\flushleft\large{\sl B.3. Sedov-Taylor Self-Similar Blast Waves}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsubsection}{\protect\numberline{V.B.3.}{\sl ~Sedov-Taylor Self-Similar Blast Waves}}

%\noindent\rm (3) Sedov-Taylor problem ({\tt lprob = 5}):\vspace{2mm}

\noindent
TITAN's hydrodynamics benchmarking is concluded with the Sedov-Taylor blast wave. An
energy deposition of $10^{50}~ergs$ at the center of a ball with radius $R_{max} = 10^9~km$
generates a spherical shock that propagates in a self-similar fashion. The gas is assumed
to have a polytropic index $\gamma = {5\over3}$ and is initially at rest with a constant
density of $10^{-8}~g/cm^3$ and temperature of $50~K$. This time, density and
pressure serve as the only grid parameters. Logarithmic scaling is choosen for all quantities
involved. The grid $\alpha$ is again kept at 1.5, the shock width is $\ell_1 = 10^{-4}$
and we take 100 zones to cover the ball. \vspace{2mm}

{\hbox{%
{\vbox{\hsize=2.5in \noindent \strut
This problem is set up so that the shock wave reaches the surface of the ball 
after $t = 2\times 10^5~sec$. TITAN completes the calculation in about 1000 time steps. The
resolution of the shock front is about $10^6$ during the evolution. In the figure to the
right the evolution of this blast wave is illustrated. Snapshots of the density and the 
velocity at five different times $t\in($140, 2900, 26590, 93023, 197670$)~sec$ are overlayed.
The upwind density is always $10^{-8}~g/cm^3$ with the material at rest. At the shock front 
the density jump remains constant $4\times 10^{-8}~g/cm^3$. The velocity jump decreases
from $117$ $km/sec$ to $1.5~km/sec$ during these five instances. Both the density and the 
velocity profiles fall off to zero behind the shock.\strut}}
{\vbox{
\psfig{figure=sedov1a.ps,width=2.4in}
\psfig{figure=sedov1b.ps,width=2.4in}
{\small Fig. 11: Density and velocity\\ profiles at five different times}}}
}}
\vspace{3mm}

\noindent
The position, density, pressure, temperature, and velocity of the shock as a function of 
time are given in next figure. All variables vary strictly monotonically except for the 
shock density which remains at a steady value $\rho_s = 4\times 10^{-8}~g/cm^3$ as it should. 
In particular the expected relations\\ 
$$       R_s(t) \propto t^{2/5}                                   \eqno{(54)}$$ 
$$       v_s(t) \propto t^{-3/5}                                  \eqno{(55)}$$
$$       T_s(t) \propto P_s(t) \propto t^{-6/5}                   \eqno{(56)}$$
are satisfied. These are the analytical solution as discussed in the textbook of Landau 
\& Lifshitz.\\

\psfig{figure=sedov2.ps,width=4.8in}
\centerline{\small Fig. 12: Power-law dependencies of the shock quantities}
\vspace{2mm}

\noindent
The self-similarity of the computed flow is demonstrated in figure~13. It shows the profiles of
density, pressure, and velocity in terms of the position, all measured with respect to the 
shock data. The graphics were generated by overlaying the actual grid values at six different 
times $t\in($140, 799, 12090, 42398, 110070, 174900$)~sec$. One can verify that the shocked
region obeys the analytical relations: 
$$          v/v_s\propto R/R_s                                    \eqno{(57)}$$
$$    \rho/\rho_s\propto (R/R_s)^{3/(\gamma-1)}                   \eqno{(58)}$$
and similarly for the pressure. Their functions jump to the preshock values and stay at a 
constant level. Compare this figure with the one in Landau \& Lifshitz!\\

\psfig{figure=sedov3.ps,width=4.8in}
\centerline{\small Fig. 13: Self-similar profiles at six different times}
 
\vspace*{2mm}
{\flushleft\large{\sl C. Radiative Transport Test Problems}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsection}{\protect\numberline{V.C.}{~Radiative Transport Test Problems}}

\noindent In the following section L. Ensman's suggestions for radiative tests are
implemented. Consider spheres with mass of about $0.5~M_\odot$ and radial extension of 
$R=2.84\times 10^3~R_\odot$ which are filled with pure hydrogen at a constant density 
$\rho = 2.9677\times 10^{-11}~g/cm^3$. Further assume that the hydrogen is fully ionized 
so that the sole opacity source is electron scattering. The mean flux opacity is then set
proportional to the Thomson scattering opacity: 
$\rho\kF = \chi_e = 1.19\times 10^{-11}~cm^{-1}$. In addition the energy and Planck mean 
opacities are equated, $\kE = \kP$, and allowed to differ from the mean flux opacity by a 
constant factor, $\kE = \fr \kF$. This setup allows us to illustrate the effects of a 
scattering versus absorption dominated gas.
 
\vspace*{2mm}
{\flushleft\large{\sl C.1. Radiative Heating}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsubsection}{\protect\numberline{V.C.1.}{\sl ~Radiative Heating}}

%\noindent\rm (1) radiative heating ({\tt lprob = 6}):\vspace{2mm}

\noindent
The sphere is set up to be in diffusion equilibrium at luminosity 
${\cal L} = 826~{\cal L}_\odot$. From the available data one computes the radiant energy 
according to\\
$$
  \er(r) = {{\cal L}\over{4\pi c R}}~
           \Bigl[\sqrt{3} + 3\chi_e\Bigl({R\over r}-1\Bigr)\Bigr]            \eqno{(59)}
$$
and from that the temperature profile via 
$$  T(r) = {^4}\!\sqrt{\er(r)/\ar}                                           \eqno{(60)}
$$ 
At time $t=0$ a strong light source at the center is turned on which has a $10^{1.5}$ times 
higher luminosity, ${\cal L}_c = 26,130~{\cal L}_\odot$, and illuminates the gas sphere.
The sphere now has to evolve to a new radiative equilibrium at this higher luminosity.\vspace{2mm}

\noindent One can calculate the relaxation with full transport under the assumption that the 
absorptive opacity is considerably smaller than the total opacity: 
$\kE = 10^{-5}~\kF = 10^{-5}~\chi_e$. This defines the case of a scattering dominated gas. 
The computation is carried out with the adaptive grid. This time the radiation energy and 
the luminosity are taken as the grid parameters in logarithmic scaling and $\alpha=3$, 
$\tau=10^3$ is set.\vspace{2mm}

{\vbox{%
\psfig{figure=heat0.ps,width=4.7in}
\centerline {\small Fig. 14: Relaxation of the radiative energy profile}
}}
\vspace{2mm}

\noindent
This graph illustrates how the higher central luminosity changes the profiles of the 
radiative energy in time. Initially the sphere is at equilibrium at the lower luminosity, 
${\cal L}_c = 826~{\cal L}_\odot$, which is represented by curve $\er(r)$ at the bottom 
($+$'s indicating individual grid points). This curve relaxes and the radiative energy 
increases until the new equilibrium at the higher luminosity, ${\cal L}_c = 
26,130~{\cal L}_\odot$, is reached. For both the initial and the final $\er(r)$ curve the
analytical solution from equation (59) were overlayed. The agreement between the analytical
(assuming diffusion) and numerical (assuming variable Eddington factors) result is remarkable.
This should be expected because the shell is optically thick except right underneath the
surface and hence the diffusion regime is established everywhere else.\vspace{2mm}

{\vbox{%
\psfig{figure=heat1.ps,width=4.7in}
\centerline{\small Fig. 15: Relaxation of gas and radiative temperatures towards equilibrium}
}}
\vspace{2mm}

\noindent
The heating of the gas ($+$) and radiative (---) temperature by the central light source is 
given in figure~15. It shows the relaxation of their profiles as a function of the optical 
depth at various times. Both bottom curves are the initial condition in which the two 
temperatures agree because of equilibrium diffusion. When the light source is turned on a 
radiation wave moves through the sphere thereby heating up the gas. Since a weak coupling 
between the gas and the radiation field is assumed, $\fr = 10^{-5}$, the gas temperature rises 
noticebly slower than the radiative temperature. As soon as the higher luminosity value has 
reached the surface both temperatures agree once again and the new radiative equilibrium is 
established. Because most of the gas is optically thick the diffusion regime prevails and the 
temperature follows again the profile from equation (59). \vspace{2mm}

\noindent
This experiment is repeated with different factors $\fr$. The sphere relaxes to the same 
equilibrium state every time. Figure~16 depicts the evolution of the surface luminosity for six
different choices $0\le \fr \le 1$. The correspondence between the symbols and the ratios $f_r$
for the gray opacity is as follows: $\fr = 0$ ($\Box$), $10^{-9}$ ($\triangle$), $10^{-7}$ 
($\Diamond$), $10^{-5}$ ($\ast$), $10^{-2}$ ($+$), 1 (---). The relaxation of the temperature 
profiles from the previous figure is represented here by
{\hbox{%
{\vbox{\hsize=2.5in \noindent \strut
the asterisk. Note that the time the radiation wave takes to reach the surface and hence begins 
to raise the luminosity is, within a factor of 2, the same in all cases, namely around 
$6\times 10^6~sec$. This number is in reasonable agreement  with the classical diffusion time 
scale $\chi_e R^2/c$ $\approx$ $1.4\times 10^7~sec$. The spheres with almost pure scattering 
opacities reach the higher luminosity plateau first. Those spheres with\strut}} 
{\vbox{
\psfig{figure=heat2.ps,width=2.4in}
{\small Fig. 16: Relaxation of the surface\\
        luminosity for six different opacities}}}
}}
\noindent
almost pure absorptive opacities attain it later.\\

\noindent The experimental setup allows us further to investigate the different radiative 
transport schemes in TITAN. There are three choices: 

\begin{description}
\item``equilibrium diffusion'':\\
The gas and radiative temperature are not distinguished, $\er = \ar T^4$, and the 
radiative flux is evaluated according to equation (7). The differentiation between the 
mean opacities becomes meaningless and only the total opacity is considered. 

\item``non-equilibrium diffusion'':\\
The radiative flux is still computed according to equation (7). But now the radiation 
energy $\er$ is computed from equation (6). Two temperatures are possible and the
absorptive character of the gas can be taken into account. Both diffusion cases demand a
constant Eddington factor $\fe = {1\over3}$. 

\item

\item``full transport'':\\
The radiation field is described by the two equations (4) \& (6). Since they necessitate 
an estimate for $\fe$ one obtains also information about the angular distribution of the 
intensity field. 
\end{description}

\noindent The consistency of the three approaches can be checked in the context of 
the hydrogen spheres. The equilibrium diffusion computation can also be simulated with
non-equilibrium diffusion by setting $\kE = \kF$. Since the gas is optically thick (aside
right underneath the surface) the corresponding full transport calculation should behave
in the same manner. Performing the actual runs demonstrates a tight agreement between the
three transport schemes.
 
\vspace*{2mm}
{\flushleft\large{\sl C.2. Radiative Cooling}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsubsection}{\protect\numberline{V.C.2.}{\sl ~Radiative Cooling}}

%\noindent\rm (2) radiative cooling ({\tt lprob = 7}):\vspace{2mm}

\noindent
The final equilibrium model from above is now taken as the initial condition. At $t=0$
the light source at the center is turned off, ${\cal L}_c = 0$. This time a fixed (Eulerian) 
grid instead of the adaptive grid is chosen. The runs are carried out with the full radiative 
transport scheme.\vspace{2mm}

{\vbox{%
\psfig{figure=cool0.ps,width=4.7in}
\centerline {\small Fig. 17: Cooling of the radiative energy profile}
}}
\vspace{2mm}

\noindent
This graph illustrates how the radiative energy decreases in time once the central luminosity 
is turned off. Initially the sphere is at equilibrium at ${\cal L}(0) = 26,130~{\cal L}_\odot$
which is represented by the curve $\er(r)$ ($+$'s again marking the grid points) at the top 
with the analytical solution overlayed. As time progresses, $\er$ decreases throughout the sphere 
as radiant energy leaks out, and the gradient of $\er$ flattens, consistent with ${\cal L}(t)$
deminishing to zero. At the surface a sharp gradient continues to drive ${\cal L}(t)$
outward. For $t\rightarrow\infty$ it goes to zero, $\er\rightarrow 0$, which represents the new 
equilibrium when ${\cal L}_c = 0$.

%{\vbox{\raggedbottom
{\hbox{%
{\vbox{\hsize=2.5in \noindent \strut
Just like before one can consider different strengths of the gas-radiation 
coupling. In figure~18 the decline of the surface luminosity for various gray opacities
defined by the ratios $0\le \fr\le 1$ are plotted. Again these ratios correspond to the
symbols in the following way: $\fr = 0$ ($\Box$), $10^{-9}$ ($\triangle$), $10^{-7}$ 
($\Diamond$), $10^{-5}$ ($\ast$), $10^{-2}$ ($+$), 1 (---). The curve of the diffusion 
regime ($\fr = 1$) is representative of all regimes with \strut}}
{\vbox{
\psfig{figure=cool2.ps,width=2.4in}
{\small Fig. 18: Relaxation of the surface\\ 
 luminosity for six different opacties}}}
}}
\noindent 
$\fr > 0$. Their analyses tells that the long time behavior of the luminosity decay follows a 
power law, ${\cal L}(t) \propto t^{-s}$, with the power being close to $s={1\over3}$. In 
contrast to that the regime of complete decoupling between the gas and radiative temperature, 
$\fr = 0$, is characterized by an exponential decay, ${\cal L}(t) \propto e^{-t/s}$, with the 
time scale $s\approx 1.7\times 10^7~sec$ being close to the classical  diffusion time. This 
regime behaves differently because the cooling wave alters only the radiative temperature but 
leaves the profile of the gas temperature unchanged.\vspace{5mm}
%}}

\noindent In figure~19 the profiles of the gas ($+$) and radiative (---) temperature for 
the case $\fr = 10^{-7}$ are plotted. Initially they are in equilibrium but the decoupling 
begins immediately. One observes that the profile of the radiative temperature is merely 
shifted downward in time. Simultaneously the profile of the gas temperature is modified to 
become a constant. Since the bulk of the total opacity comes from electron scattering the 
radiation cools off faster than the gas so that the gas energy contribution to the luminosity 
rises soon to domination. Note that the radiative temperature is related to the radiative
energy via equation (60) and therefore shows the same behaviour as $\er$ in figure 17.

\psfig{figure=cool1.ps,width=4.7in}
\centerline{\small Fig. 19: Relaxation of gas and radiative temperatures towards equilibrium}
 
\vspace*{2mm}
{\flushleft\large{\sl D. Radiation Hydrodynamics Test Problems}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsection}{\protect\numberline{V.D.}{~Radiation Hydrodynamical Test Problems}}

\noindent This section describes three radiation hydrodynamics benchmarks. The fundamental 
understanding of radiating shocks was developed by Zel'dovich \& Raizer and additional 
discussion of the phenomena involved is given by Mihalas \& Mihalas. All tests are carried out
in spherical geometry. The full transport scheme for the radiative field is employed. 
L. Ensman's test setups are again followed in detail.\vspace{2mm}

\noindent Consider a thin shell with an inner radius $R_i = 8\times 10^6~km$ and outer
radius $R_o = 8.7\times 10^6~km$ filled with a cold gas at constant density $7.78\times
10^{-10}~g/cm^3$. A shallow temperature profile 
$$           T(r) = 10 + 75~{{r-R_o}\over{R_i-R_o}}~K                 \eqno{(61)}
$$
is assumed . The Rosseland mean opacity is taken to be constant 
$\rho\kR=3.115\times 10^{-10}~cm^{-1}$, and all of it is absorptive, $\fr=1$. The gas is 
initially at rest. At time $t=0$ a piston at the inner radius is moved outwards at a 
constant high speed. A shock front is thus formed which couples strongly to the radiation 
field. If the maximum temperature immediately ahead of the shock is less than the downstream 
temperature one speaks of a ``subcritical'' shock, the case of equality one speaks of a
``supercritical'' shock.\vspace{2mm}

\noindent In the following calculations the sphere is represented by 100 grid points. 
Density and mass in logarithmic scaling are selected as the sole grid parameters along with
$\alpha=1.5$ and $\tau = 10^{-20}$. The shock width is taken to be proprotional to
the radius, $\ell_1 = 10^{-5}$, to demand a high resolution. All runs are executed fully 
implicit, \ie~with time centering parameter $\theta=1$.
 
\vspace*{2mm}
{\flushleft\large{\sl D.1. Subcritical Shock}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsubsection}{\protect\numberline{V.D.1.}{\sl ~Subcritical Shock}}
 
%\noindent\rm (1) subcritical shock ({\tt lprob = 8}):\vspace{2mm}

\noindent
The piston speed is $6~km/sec$. Figure~20 shows the radiative shock front propagating 
towards the surface. The four quantities gas $T$ and radiative temperature $T_r$, 
radiative flux $\Fr$, and density $\rho$ are plotted as functions of the radius at seven 
different times. The shock front is located where the flux is largest. The temperatures 
just ahead of the shock are always smaller than those in the shocked region thus defining 
this run as a subcritical shock. \vspace{2mm}

{\vbox{%
\psfig{figure=sub1.ps,width=4.7in}
\centerline{\small Fig. 20: Propagation of a subcritical shock}
}}
\vspace{2mm}

{\hbox{%
{\vbox{\hsize=2.5in \noindent \strut
It is often useful to look at the dependencies of the physical quantities on other coordinates.
The following figures illustrate this for the example of the gas temperature. The same seven
snapshots were taken but instead of showing them as a function of the radius they are now plotted
as a function of the exterior mass, of the optical depth, and of the grid point (zone 
index). The dependency $T=T(\Me)$ is the Lagrangean representation of the data which emphazises 
the outermost layers of the gas sphere. The graph $T=T(\tau)$ shows the temperature profile 
from a radiation transport perspective which stresses how transparent the features are. These 
functional dependencies help the physical interpretation of the results. On the other hand, 
the profile $T=$ \{$T_k$ ; $k=1,\dots,100$\} is interesting from a numerical point of view 
because one can learn from them how well the calculation was carried out. In the given example 
one can, for instance, see that the temperature jumps in the shock are well resolved by about 
10 grid points.\strut}}
{\vbox{%
\psfig{figure=sub1a.ps,width=2.4in}
\psfig{figure=sub1b.ps,width=2.4in}
\psfig{figure=sub1c.ps,width=2.4in}
{\small Fig. 21: Temperature profiles}}}
}}
\vspace{2mm}

\noindent
Its characteristic features, as explained in Zel'dovich \& Raizer, are an overshoot of the 
gas temperature, a nearly symmetric profile of the radiative flux, and an extended 
non-equilibrium region upstream from the shock, which is generated by the effects of radiative 
preheating from the flux. In figure~22 the density, gas and radiative temperature, flux 
and Eddington factor ($\cdots$, scaled so that 1.0 corresponds tp $\fe = {1\over3}$), and velocity 
are plotted in dependency on the optical depth away from the front, $(\tau-\tau_s)$, at the half 
piston crossing time ($t = 56,581~sec$ to be exact). \vspace{2mm}

\psfig{figure=sub2.ps,width=4.7in}
\centerline{\small Fig. 22: Profiles as functions of the relative optical depth}
\vspace{2mm}

\noindent
A non-equilibrium diffusion analysis, \cf~Zel`dovich \& Raizer, demonstrates that the density, gas
pressure, radiative flux and energy all decrease upstream depend exponentially on the optical
depth away from the shock:
$$  \rho(\tau) \propto \Pg(\tau) \propto \Fr(\tau) \propto \er(\tau) 
	       \propto e^{-\sqrt{3}|\tau-\tau_s|}                  \eqno{(62)}
$$ 
One consequence is that the radiative temperature must decrease slower than the gas temperature. 
This functional dependency on the relative optical depth holds also downstream for the flux and 
the gas temperature. The maximum flux is determined by the postshock temperature $T_2\;$:
$F_s \approx {1\over2\sqrt{3}}\ar c T_2^4$. The analytical values differ a factor less than
2 from the actual numerical data. The behavior of our simulation agrees quite well with the
analytical estimates, and the existing deviations can be explained by the use of variable
Eddington factors.
 
\vspace*{2mm}
{\flushleft\large{\sl D.2. Supercritical Shock}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsubsection}{\protect\numberline{V.D.2.}{\sl ~Supercritical Shock}}

%\noindent\rm (2) supercritical shock ({\tt lprob = 9}):\vspace{2mm}

\noindent
The piston speed is $20~km/sec$. Figure~23 shows the radiating shock while it propagates 
across the shell. The temperatures in the regions immediately ahead of and trailing the shock 
front are equal. Therefore this run is classified as a supercritical shock. \vspace{2mm}

{\hbox{%
\psfig{figure=sup1.ps,width=4.7in}}}
\centerline{\small Fig. 23: Propagation of a supercritical shock}
\vspace{2mm}

\noindent
Its characteristic features are again a strong overshoot of the gas temperature, this time a 
highly asymmetric flux profile, and an extended radiative precursor whose head is out of 
equilibrium with the gas. In contrast to the subcritical case the immediate surroundings of 
the shock satisfy the condition of radiative equilibrium. The temperature overshoot itself 
is out of equilibrium and its width is narrower than the unit mean free path of the photons. 
This implies that the radiation flow from this feature cannot be treated correctly by the
diffusion approximation. In figure~24 the density, gas and radiative temperature, flux and 
Eddington factor ($\cdots$ scaled so that 1 corresponds to $\fe = {1\over3}$), and velocity 
are plotted in dependency on the optical depth away from the front, $\tau-\tau_s$, 
at one quarter piston crossing time ($t = 8,541~sec$). \vspace{2mm}

\noindent
A non-equilibrium diffusion analysis according to Zel'dovich \& Raizer predicts again that
the radiative flux and radiative energy decline exponentially and so does the gas 
temperature in the upstream part of the precursor that is out of equilibrium. In the 
equilibrium region around the shock one finds for the gas temperature a dependency on 
the optical depth away from the front like
$$  T(\tau) \propto ~{^3}\!\sqrt{1 + {3\sqrt{3}\over4}|\tau-\tau_s|}       \eqno{(63)}
$$

\noindent 
There is also an estimate for the maximum flux which is not as accurate as for the 
subcritical case. Radiative equilibrium is well established within $(\tau-\tau_s)<-5$ aside 
of the shock front itself. Density, flux, and velocity undergo jumps from the up- to the 
downstream regions. Our numerical results reflect 
the analytical predictions. Because the analytical results are not as reliable as in the 
previous case a close comparison does not give good results despite the fact we actually 
have $\fe = {1\over3}$.\vspace{2mm}

{\vbox{%
\psfig{figure=sup2.ps,width=4.7in}
\centerline{\small Fig. 24: Profiles as function of the relative optical depth}
}}
\vspace{2mm}

\noindent 
The spike in the gas temperature is prominent and very narrow. Its width is given by 
$(\rho\chi_e)^{-1}$ $=$ $3,210~km$. This is much larger than the rezoning limitation of 
$\ell \approx 80~km$ imposed by the artificial shock width. In figure~25 on the right gas 
($+$) and radiative ($\Diamond$) temperature are plotted in the vicinity of the shock 
front. Their dependency in the relative optical depth is almost independent on time. The 
spike exists clearly only for the gas temperature while the radiative temperature varies 
smoothly. The spike is optically thin

{\hbox{%
{\vbox{\hsize=2.5in\noindent\strut
and hence practically transparent to the radiation (although the actual location is mostly 
at large optical depths). Note that there are many (around 25) grid points in this region, 
each zone with an optical depth of at most $5\times 10^{-3}$. The spike is therefore 
numerically well represented and represents the physical, in particular radiative, 
properties correctly.
\strut}}
{\vbox{
\psfig{figure=sup3.ps,width=2.4in}
{\small Fig. 25: Temperature spike}}}
}}
 
\vspace*{5mm}
{\flushleft\large{\sl D.3. Radiative Blast Wave}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsubsection}{\protect\numberline{V.D.3.}{\sl ~Radiative Blast Wave}}

%\noindent\rm (3) radiative blast ({\tt lprob = 10}):\vspace{2mm}

\noindent
Now consider a shell stretching from an inner radius $R_i = 10^7~km$ to an outer radius
of $R_o = 10^8~km$ containing a mass of $0.16~M_\odot$. Instead of a constant density
profile a power law atmosphere is assumed:
$$      \rho(r) = 10^{-6}~\Bigl({R_i\over r}\Bigr)^7~g/cm^3              \eqno{(64)}
$$
The gas is purely absorptive and the
mean opacity taken to be constant and equal to the Thomson free electron scattering opacity
(hence implying the gas is completely ionized hydrogen). The adiabatic index is
$\gamma={5\over3}$. The sphere is illuminated by a central light source with a luminosity
${\cal L} = 10^{34}~ergs/sec$. A self-consistent temperature profile for the initial model
is obtained by relaxating the sphere's luminosity to this value using the full transport 
scheme, \cf~the discission of radiative heating.\vspace{2mm}
 
\noindent At time $t=0$ a large amount of energy, corresponding to a central temperature of 
$10^9~K$, is deposited inside the sphere. This sets off a radiative blast wave (differing
from the Sedov-Taylor blast wave by the fact that now the effects of radiation are included).
Its evolution is followed employing the adaptive grid. As grid parameters we choose the mass, 
density, and radiative energy in logarithmic and the velocity $u$ ($1000~km/sec$) in linear 
scaling, together with $\alpha=1.5$, $\tau=10^{-20}$, and $\theta=1$. The shock width is 
$\ell_1=10^{-4}$. After initiating the blast wave it reaches the surface in 700 
time steps. \vspace{2mm}

{\vbox{%
\psfig{figure=rblast1.ps,width=4.7in}
\centerline{\small Fig. 26: Evolution of a radiative blast wave}
}}
\vspace{4mm}

\noindent In figure~26 the evolution of the gas sphere is depicted at six different times. 
Initially the shock moves down the density and temperature gradient. There is a very prominent 
flux peak and the shock profiles for density, temperature and velocity are sharp. The shock 
looks like an adiabatic one. In time the

{\hbox{%
{\vbox{\hsize=2.5in\noindent\strut
radiation preheats an increasing region ahead of it. When the shock has penetrated to an 
optical depth of about 200 (which occurs at time $t\approx 7000~sec$) this radiative precursor 
has reached the surface. The surface luminosity rises very fast and reaches its maximum at 
about $t=1.13\times 10^4~sec$, see figure~27. The gas and radiative temperatures are out of 
equilibrium at the head of the precursor. The heating induces a\strut}}
{\vbox{
\psfig{figure=rblast2.ps,width=2.4in}
{\small Fig. 27: Light curve of the blast}}}
}}
\noindent
quasi-isothermal wind in the outer shell and the gas equilibrates again. Before reaching an 
optical depth of 60 the shock has dissipated through the radiative processes. The shocked gas, 
however, maintains its momentum and begins a ``snow plow'' phase where it sweeps up wind material. 
The surface luminosity declines and an almost constant flux profile becomes established. The 
sphere begins to cool off. Behind the ``snow plow'' a very dense shell is formed which is well 
established by $t\approx 1.5\times 10^4~sec$ at optical depth $\tau = 16$. While the plow steepens 
up at its head to a strong supercritical shock (with the typical temperature spike) another,
reverse shock forms at its tail. It propagates slowly into the homologously expanding inner
regions. The dense shell expands rapidly and moves at nearly the shock speed outwards.
Through the radiation the gas keeps cooling down, and, since the almost
constant flux level decreases, the surface luminosity keeps dropping.\vspace{2mm}

\psfig{figure=rblast3.ps,width=4.7in}
\centerline{\small Fig. 28: Density contours (logarithmic scale)}
\vspace{2mm}

\noindent The evolution of the radiative blast simulation can be illustrated via ``optical 
depth-time diagrams'' of the density, velocity, and temperature (figures~28 -- 30). One can 
recognize the adiabatic shock during the initial phase up to several 1000 $sec$. While it 
dissipates through radiative losses the radiative precursor surfaces at about $10^4~sec$ and 
the luminosity there breaks out. The shocked regions plow into the isothermal wind while 
the gas keeps cooling off. A new supercritical shock is born and runs fast towards the surface. 
Simultaneously a reverse adiabatic shock forms. It moves into the inner regions which already 
undergo a homologous expansion and a high density shell at a constant speed is created.

{\vbox{\raggedbottom
\psfig{figure=rblast4.ps,width=4.7in}
\centerline{\small Fig. 29: Temperature contours (logarithmic scale)}
\vspace{5mm}
\psfig{figure=rblast5.ps,width=4.7in}
\centerline{\small Fig. 30: Velocity contours (linear scale)}
}}


\noindent
The grid motion is depicted in figure~31 in which every fourth zone is plotted 
as a function of optical depth and time. After a very short relaxation period the grid 
points begin to zoom into the forming shock, $t\approx 2000~sec$. They actually lock into 
the radiative precursor until it surfaces at $t\approx 7000~sec$. The grid then rezones 
quickly to resolve a new shock system that forms around $t\approx 10^4~sec$. It tracks a 
reverse shock and follows a forward shock all the way to the surface.\\

\psfig{figure=rblast6.ps,width=4.7in}
\centerline{\small Fig. 31: Action of the adaptive grid}
 
%\end{document}
.
