%\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} .