\subsection{Hydrodynamical Test Problems} \subsubsection{Noh's shock tube} In this shock tube problem a piston of constant speed $v = -1$ pushes a cold gas with $\gamma = {5\over3}$ into the origin. Courant \& Friedrichs~[1] have discussed this problem qualitatively and showed that an instantaneous shock of infinite Mach number is generated. The shock speed turns out to be ${1\over2}(\gamma-1)|v|={1\over3}$. Noh~[9] defined this problem as a challenging benchmark for hydrodynamics codes that can be run in cartesian, cylindrical, and spherical geometry. In addition, there are analytical formulae against which the numerical results can be tested. We follow Noh's examinations at time $t=0.6$ before the shock collides with the piston. We represent the computational domain with 100 grid points and employ the density, pressure and energy in logarithmic scaling as our grid parameters together with $\alpha = 1.5$ and $\tau = 10^{-8}$. We obtain solutions to the Noh problem in all three geometries. Let us first consider the planar case which we calculate with a fixed and the adaptive grid. The superiority of the adaptive grid is well established: since the shock is the only nonlinear feature in the flow it is resolved with some 60 zones. The time $t=0.6$ is reached after 160 time steps while constraining the time step to $\le 10^{-2}$. The adaptive grid is locked into the shock frame and thus the code evolved almost statically. In contrast to that the shock has to move from one grid point to the next in the Eulerian computation thus allowing the time step to only vary between $10^{-5}$ and $10^{-4}$ (note that this still orders of magnitudes larger than the practical CFL number!). This computation exhibits all the problems induced by the use of the artificial viscosity as discussed in~[9]. The adaptive grid, however, overcomes these difficulties very early on ($t<10^{-2}$). % after about 90 time steps. \psfig{file=noh.ps,width=\textwidth} \vspace{-2.0in} {\tenrm Fig. 1 (Noh's shock tube): density ($+$), energy ($\Diamond$), and pressure ($\triangle$) profiles in planar, cylindrical,} \hspace{10mm}{\tenrm and spherical geometry.} \vspace{5mm} In figure~1 we plot profiles of the density, pressure, and energy at time $t=0.6$ for planar, cylindrical, and spherical geometry. The adaptive grid behaves very similar in all cases producing sharp steps at the shock front. The step in energy is always ${1\over2}$ because of the isothermal nature of the shock. The step in density equals to $3\cdot4^\mu$. The density profiles in the unshocked regions satisfy the analytical relations: $\rho(r,t) = (1+t/r)^\mu$ which are valid for $t<{3\over4}$ and $r>{1\over3}t$. In the shocked regions the density plateau at $\rho(r,t) = 4^{\mu+1}$ for $0\le r\le {1\over3}t$ is well maintained in planar geometry and closely established in cylindrical and spherical geometry. We have evaluated the relative errors and found them to be between one tenth of a and a few percent in all cases. Since the theoretical plateau rises rapidly with $\mu$ the error becomes most visible in spherical geometry. We would like to remind that fixed grid computations usually generate errors of order unity for this problem! Therefore their solutions do not resemble the correct answer. The adaptive grid allows, for the first time, to generate good numerical results with little computational effort. Next we simulate two different kinds of blast waves, one in planar and the other in spherical geometry. \subsubsection{Woodward's interacting blast waves} First we repeat a planar problem of two interacting blast waves that was invented by Woodward and is described in~[14]. The computational domain is divided into three sections with three different but constant pressure regimes in them. In the left part we have $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$ with an ideal gas resting at unit density. The adiabatic index is $\gamma = 1.4$. The density is initially constant, $\rho = 1$, across the domain $0\le r\le 1$. We adopt again the same grid parameters as before but this time run the code with 200 grid points. The artificial viscosity is kept linear at $\ell_o = 10^{-4}$. \vspace{2mm} {\hbox{% {\vbox{\hsize=3.0in \strut The large pressure steps at both ends of the cylinder create two shocks which travel towards the middle. The graph to the right (figure 2) 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{\hsize=10mm}} {\vbox{% \psfig{figure=wood0.ps,width=3.0in} {\tenrm Fig. 2 (Woodward's blast waves): density profile} \hspace{10mm}{\tenrm at time t = 10$^{-4}$.} }} }} \vspace{5mm} 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.\\ \centerline{\psfig{figure=wood1.ps,width=5.5in}} {\tenrm Fig. 3 (Woodward's blast waves): shock profiles at four different times ($+$'s indicate grid points).} \vspace{5mm} 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 4. 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.\\ \centerline{\psfig{figure=wood2.ps,width=5.5in}} {\tenrm Fig. 4 (Woodward's blast waves): space-time diagram of the logarithmic density; 25 contour lines have} \hspace{10mm}{\tenrm been selected equidistantly.} \vspace{5mm} The five density profiles (figures 2 \& 3) 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} 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 our code 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.\\ \centerline{\psfig{figure=wood3.ps,width=5.5in}} {\tenrm Fig. 5 (Woodward's blast waves): space-time diagram of the grid motion; the trajectory of every fourth} \hspace{10mm}{\tenrm zones is depicted.} \vspace{5mm} Figure~5 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. In conclusion of this exercise we can highlight the generic behavior of the adaptive grid: Once the nonlinear features of the flow are resolved the adaptive grid locks into them and the code can evolve quasi-statically at time steps that are considerably larger than those of a fixed grid. When the grid needs to redistribute, as during the reflection of the shock from the wall or its interaction with the other shock or contact discontinuity, the time step drops by several orders of magnitude to accomplish this in a very short time interval. \subsubsection{Sedov-Taylor self-similar blast wave} We conclude our hydrodynamics benchmarking 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=3.0in \strut This problem is set up so that the shock wave reaches the surface of the ball after $t = 2\times 10^5~sec$. The code 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{\hsize=10mm}} {\vbox{% \psfig{figure=sedov1a.ps,width=3.0in} \psfig{figure=sedov1b.ps,width=3.0in} {\tenrm Fig. 6 (Sedov-Taylor blast wave): density and} \hspace{10mm}{\tenrm velocity profiles at five consecutive times.}}}\\ }} \vspace{5mm} 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{(31)}$$ $$ v_s(t) \propto t^{-3/5} \eqno{(32)}$$ $$ T_s(t) \propto P_s(t) \propto t^{-6/5} \eqno{(33)}$$ are satisfied. These are the analytical solution as discussed in the textbook of Landau \& Lifshitz.\\ \centerline{\psfig{figure=sedov2.ps}} {\tenrm Fig. 7 (Sedov-Taylor blast wave): position R$_s$~(km), density $\rho_s$~(10$^{-8}$ g/cm$^3$), temperature T$_s$~(K),} \hspace{10mm}{\tenrm pressure P$_s$~(dyn/cm$^2$), and speed v$_s$~(km/sec) of the shock as a function of time; notice the} \hspace{10mm}{\tenrm power law dependencies.} \vspace{5mm} The self-similarity of the computed flow is demonstrated in figure~8. 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{(34)}$$ $$ \rho/\rho_s\propto (R/R_s)^{3/(\gamma-1)} \eqno{(35)}$$ 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!\\ \centerline{\psfig{figure=sedov3.ps}} {\tenrm Fig. 8 (Sedov-Taylor blast wave): self-similarity in the profiles of the density $\rho/\rho_s$, pressure P/P$_s$, and} \hspace{10mm}{\tenrm velocity v/v$_s$ generated through overlays from six different times t $\in$ (140, 532, 576, 1160, 2290,} \hspace{10mm}{\tenrm 4540) sec.} \subsection{Radiation Test Problems} In this section we follow L. Ensman's suggestions for radiative tests~[4]. We consider spheres of about $0.5~M_\odot$ and radial extension of $R=2.84\times 10^3~R_\odot$ that are filled with pure hydrogen at a constant density $\rho = 2.9677\times 10^{-11}~g/cm^3$. Further we assume the hydrogen to be fully ionized so that the opacity source are scattering electrons; hence we set the mean flux opacity proportional to the Thomson scattering opacity: $\rho\kF = \chi_e = 1.19\times 10^{-11}~cm^{-1}$. In addition we always equate the energy and Planck mean opacities, $\kE = \kP$, and allow them to differ from the mean flux opacity by a constant factor, $\kE = \fr \kF$, in order to illustrate the effects of a scattering versus absorption dominated gas. \subsubsection{Radiative heating tests} We set up this sphere to be in diffusion equilibrium at a luminosity ${\cal L} = 826~{\cal L}_\odot$. From the available data one computes the radiative 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{(36)} $$ \noindent and from that the temperature profile via $T(r) = {^4}\!\sqrt{\er(r)/\ar}$. At time $t=0$ we turn on a strong light source at the center with a $10^{1.5}$ times higher luminosity, ${\cal L}_c=26,130~{\cal L}_\odot$, which illuminates the gas sphere. The sphere now has to evolve into a new radiative equilibrium at this higher luminosity. \vfill \centerline{\psfig{figure=heat0.ps,width=5.5in}} {\tenrm Fig. 9 (Heating towards radiative equilibrium): profiles of the radiant energy ($+$) at different times} \hspace{10mm}{\tenrm during the relaxation of a scattering dominated atmosphere with $\rho\kP$ = 10$^{-5}$ $\chi_e$.} \vspace{5mm} 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 (36) 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} \centerline{\psfig{figure=heat1.ps,width=5.5in}} {\tenrm Fig. 10 (Heating towards radiative equilibrium): profiles of the gas ($+$) and radiation (---) temperature} \hspace{12mm}{\tenrm at different times during the relaxation of a scattering dominated atmosphere with $\rho\kP$ = } \hspace{12mm}{\tenrm 10$^{-5}$ $\chi_e$.} \vspace{5mm} The heating of the gas ($+$) and radiative (---) temperature by the central light source is given in figure~10. 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 (36). \vspace{2mm} This experiment is repeated with different factors $\fr$. The sphere relaxes to the same equilibrium state every time. Figure~11 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 (---). \break The relaxation of the temperature profiles from the previous figure is represented here by the asterisk. Note that the time the radiation wave takes to reach the surface and hence be- {\hbox{% {\vbox{\hsize=3.0in \noindent \strut gins to raise the luminosity is, with- in 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 almost pure absorptive opacities attain it later.\strut}} {\vbox{\hsize=10mm}} {\vbox{% \psfig{figure=heat2.ps,width=3.0in} {\tenrm Fig. 11: Relaxation of the surface luminosity for} \hspace{13mm}{\tenrm six different opacities}}} }} \vspace{5mm} 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``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} 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. We calculate the relaxation with the full transport scheme under the assumption that the absorptive opacity is considerably smaller than the total opacity: $\rho\kE = 10^{-5}~\rho\kF = 10^{-5}~\chi_e$. This defines the case of a scattering dominated gas. The computation is carried out with our adaptive grid. This time we take the radiation energy and the luminosity as the grid parameters in logarithmic scaling and set $\alpha=3$, $\tau=10^3$. \subsubsection{Radiative cooling tests} We now take the equilibrium model described above as the initial condition. At $t=0$ we turn the light source at the center off. This time we choose a fixed (Eulerian) grid instead of the adaptive grid and perform the runs with the full radiative transport scheme. \vfill \centerline{\psfig{figure=cool0.ps,width=5.5in}} {\tenrm Fig. 12 (Cooling away from radiative equilibrium): profiles of the radiant energy ($+$) at different times} \hspace{12mm}{\tenrm during the relaxation of a scattering dominated atmosphere with $\rho\kP$ = 10$^{-7}$ $\chi_e$.} \vspace{5mm} 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$. Just like before one can consider different strengths of the gas-radiation coupling. {\hbox{% {\vbox{\hsize=3.0in \noindent \strut In figure~13 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 $\fr > 0$. Their analyses tells that the long time behavior \strut}} {\vbox{\hsize=10mm}} {\vbox{% \psfig{figure=cool2.ps,width=3.0in} {\tenrm Fig. 13: Relaxation of the surface luminosity for} \hspace{13mm}{\tenrm six different opacties}}} }} \noindent 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} In figure~14 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 (36) and therefore shows the same behaviour as $\er$ in figure 12. \centerline{\psfig{figure=cool1.ps,width=5.5in}} {\tenrm Fig. 14 (Cooling away from radiative equilibrium): profiles of the gas ($+$) and radiation (---) temper-} \hspace{12mm}{\tenrm ature at different times for a scattering dominated atmosphere with $\rho\kP$ = 10$^{-7}$ $\chi_e$.} \subsection{Radiation Hydrodynamics Test Problems} The fundamental understanding of radiating shocks was developed by Zel'dovich \& Raizer and additional discussion of the phenomena involved is given by Mihalas \& Mihalas. This section describes three radiation hydrodynamics benchmarks. All tests are carried out in spherical geometry. The full transport scheme for the radiative field is employed. We follow again L. Ensman's test setups.\vspace{2mm} We first 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$. Further, we assume a shallow temperature profile $$ T(r) = 10 + 75~{{r-R_o}\over{R_i-R_o}}~K . \eqno{(37)} $$ 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} In the following calculations we represent the sphere 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$. \subsubsection{Subcritical shock} The piston speed is $6~km/sec$. Figure~15 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. \vfill \centerline{\psfig{figure=sub1.ps,width=5.5in}} {\tenrm Fig. 15 (Subcritical shock): evolution of the gas and radiative temperature, the radiative flux, and the } \hspace{12mm}{\tenrm velocity profiles while a piston of speed 6 km/sec pushes the left boundary towards right.} \vspace{5mm} {\hbox{% {\vbox{\hsize=3.0in \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 \strut}} {\vbox{% \psfig{figure=sub1a.ps,width=3.0in}}} }} {\hbox{% {\vbox{\hsize=3.0in \noindent \strut 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{\hsize=10mm}} {\vbox{% \psfig{figure=sub1b.ps,width=3.0in} \psfig{figure=sub1c.ps,width=3.0in} {\tenrm Fig. 16 (Subcritical shock): three different repres-} \hspace{12mm}{\tenrm entations of the evolution of the tempera-} \hspace{12mm}{\tenrm ture profile}}} }} \vspace{5mm} 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~17 the density, gas and radiative temperature, flux and Eddington factor (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} \centerline{\psfig{figure=sub2.ps,width=5.5in}} {\tenrm Fig. 17 (Subcritical shock): snap shots of the density, gas ($+$) and radiative ($\Diamond$) temperature, radiative} \hspace{12mm}{\tenrm flux, and gas velocity as functions of the optical depth $\tau$ relative to the optical depth of the} \hspace{12mm}{\tenrm shock $\tau_s$ at half piston crossing time; overlayed over the flux plot is three times the Eddington} \hspace{12mm}{\tenrm factor ($\cdots$).} \vspace{5mm} 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{(38)} $$ 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. \subsubsection{Supercritical shock} The piston speed is $20~km/sec$. Figure~18 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} \centerline{\psfig{figure=sup1.ps,width=5.5in}} {\tenrm Fig. 18 (Supercritical shock): evolution of the gas and radiative temperature, the radiative flux, and} \hspace{12mm}{\tenrm the velocity profiles while a piston of speed 20 km/sec pushes the left boundary towards right.} \vspace{5mm} 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 (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} 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{(39)} $$ 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}$. \vfill \centerline{\psfig{figure=sup2.ps,width=5.5in}} {\tenrm Fig. 19 (Supercritical shock): snap shots of the density, gas ($+$) and radiative ($\Diamond$) temperature, radia-} \hspace{12mm}{\tenrm tive flux, and gas velocity as functions of the optical depth $\tau$ relative to the optical depth of} \hspace{12mm}{\tenrm the shock $\tau_s$ at quater piston crossing time; overlayed over the flux plot is three times the} \hspace{12mm}{\tenrm Eddington factor ($\cdots$).} \vspace{5mm} 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 {\hbox{% {\vbox{\hsize=3.0in\noindent\strut thin 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{\hsize=10mm}} {\vbox{% \psfig{figure=sup3.ps,width=3.0in} {\tenrm Fig. 20: Temperature spike of supercritical shock}}} }} \subsubsection{Radiative blast} We 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 we assume a power law atmosphere: $$ \rho(r) = 10^{-6}~\Bigl({R_i\over r}\Bigr)^7~g/cm^3 \eqno{(40)} $$ The gas is purely absorptive and the mean opacity is taken to be constant and equal to the Thomson free electron scattering opacity (hence we imply the gas to be 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$. We obtain a self-consistent temperature profile for the initial model by relaxating the sphere's luminosity to this value using the full transport scheme, \cf~the discission of radiative heating.\vspace{2mm} At time $t=0$ we deposit a large amount of energy inside the sphere (corresponding to a central temperature of $10^9~K$). 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} \centerline{\psfig{figure=rblast1.ps,width=5.5in}} {\tenrm Fig. 21 (Radiative blast): snap shots of the density, gas temperature, radiative flux, and gas velocity} \hspace{12mm}{\tenrm as functions of the exterior mass log M$_{ext}$ at six different times t $\in$ (437, 2162, 7134, 10329,} \hspace{12mm}{\tenrm 15562, 29099) sec; the corresponding optical depths of the shock front are $\tau\in$ (3.2$\times$10$^4$,} \hspace{12mm}{\tenrm 6.25$\times$10$^3$, 184, 61.5, 15.9, 0.17).} \vspace{5mm} In figure~21 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 radiation preheats an increasing region ahead of {\hbox{% {\vbox{\hsize=3.0in\noindent\strut 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~22. The gas and radiative temperatures are out of equilibrium at the head of the precursor. The heating induces a quasi-\strut}} %{\vbox{\hsize=10mm}} {\vbox{% \psfig{figure=rblast2.ps,width=3.0in} {\tenrm Fig. 22: surface luminosity as a function of time;} \hspace{12mm}{\tenrm the $+$ signs mark individual time steps.}}} }} \noindent 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} \centerline{\psfig{figure=rblast3.ps,width=5.5in}} {\tenrm Fig. 23 (Radiative blast): optical depth-time diagram of the logarithmic density in g/cm$^3$.} \vspace{5mm} The evolution of the radiative blast simulation can be illustrated via ``optical depth-time diagrams'' of the density, velocity, and temperature (figures~23 -- 25). 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. \centerline{\psfig{figure=rblast4.ps,width=5.0in}} {\tenrm Fig. 24 (Radiative blast): optical depth-time diagram of the velocity in 1000 km/sec.} \vspace{5mm} \centerline{\psfig{figure=rblast5.ps,width=5.0in}} {\tenrm Fig. 25 (Radiative blast): optical depth-time diagram of the logarithmic gas temperature in K.} \vspace{5mm} The grid motion is depicted in figure~25 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.\\ \centerline{\psfig{figure=rblast6.ps,width=5.5in}} {\tenrm Fig. 26 (Radiative blast): optical depth-time diagram of the grid point motion; every fourth zone is} \hspace{12mm}{\tenrm shown.} .