{'fraction_non_alphanumeric': 0.03681458048335613, 'fraction_numerical': 0.01037391700866393, 'mean_word_length': 4.450390276105782, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The RAGE radiation-hydrodynamic code

We describe RAGE, the ``Radiation Adaptive Grid Eulerian'' radiation-hydrodynamics code, including its data structures, its parallelization strategy and performance, its hydrodynamic algorithm(s), its (gray) radiation diffusion algorithm, and some of the considerable amount of verification and validation efforts. The hydrodynamics is a basic Godunov solver, to which we have made significant improvements to increase the advection algorithm's robustness and to converge stiffnesses in the equation of state. Similarly, the radiation transport is a basic gray diffusion, but our treatment of the radiation-material coupling, wherein we converge nonlinearities in a novel manner to allow larger timesteps and more robust behavior, can be applied to any multi-group transport algorithm.


Introduction
RAGE, an acronym for Radiation Adaptive Grid Eulerian, is a multidimensional (1D, 2D, 3D), multi-material, massively parallel, Eulerian hydrodynamics code for use in solving the Euler equations coupled with a radiation diffusion equation (gray or multigroup) in a variety of highdeformation flow problems. It uses computational cells of unit aspect ratio (squares in Cartesian 1D, 2D, 3D; tori in cylindrical 1D and 2D, and shells in spherical 1D). The mesh is refined (or unrefined) locally every computational cycle by an Adaptive Mesh Refinement (AMR) algorithm that looks at local spatial variation of physical quantities (pressure, density, material boundaries, etc.) and then bisects or recombines cells so as to maintain a desired degree of resolution. RAGE is intended for general applications without tuning of algorithms or parameters, and is portable enough to run on a wide variety of platforms, from desktop PC's (Windows, Linux, and Macintosh) to the latest MMP and SMP supercomputers.
After a brief historical review, the next section will discuss methods of generating meshes from various geometry descriptions and will introduce certain data-structure concepts, while section 2 will discuss our parallel-centric data structures in greater depth, because they allow us to scale to kilo-pe's without interfering with the physics. Section 3 will discuss the adaption logic briefly, because it uses a very conservative (i.e. mesh-profligate) algorithm, whose only merit is that all material boundaries and sharp gradients will only be found at the finest level of the mesh.
Because our Eulerian hydrodynamic algorithm has more of a Lagrangian flavor than most Eulerian codes, section 4 will describe the relation between the Riemann solution and the advection in detail. Our equation of state lookup techniques will be briefly discussed in section 5, and the logic that is required to modify the Riemann solver to handle self-gravity is covered in section 6.
RAGE's novel numerical treatment of the energy coupling between matter and radiation is covered in section 7. The exponential differencing of the heating term allows a smooth transition to equilibrium diffusion. The code uses an iterative solution to the nonlinear terms in this part of the equations, allowing larger timesteps and less computer time.
We conclude with a brief discussion of some of the work done at Los Alamos National Laboratory to verify the correctness of the algorithms and validate the appropriateness of those algorithms for a broad range of problems of interest to that community.
A turbulence model and a multifluid hydrodynamic treatment are under development at LANL, and work has also begun on improved AMR methods. We intend to discuss these and other work (e.g. multigroup radiation, ionization physics) in future articles.

Historical Background
The SAGE (SAIC Adaptive Grid Eulerian) code is a multidimensional (1D, 2D, 3D), multimaterial, massively parallel, Eulerian hydrodynamics code for use in solving the Euler equation (with optional thermal conduction) in a variety of high-deformation flow problems, begun by M. Gittings in 1990 to support DNA/DTRA work [1,2].
The RAGE code is built on top of SAGE by adding various radiation packages (gray diffusion and multigroup diffusion packages were written by T. Betlach and N. Byrne as an Internal Research and Development project in 1990-1991 [3]) and is used for problems in which radiative transport of energy is important.
RAGE is the result of a long prior evolution from earlier incarnations as a modeling code for ground-effects experiments at the Nevada Test Site and for underwater blast-effects. It was brought to Los Alamos National Laboratory in 1995 in order to develop an AMR Eulerian code suitable for problems of interest to that community. The original code was written for the supercomputers of the early 1990's: Cray vector processors. We found that these machines were not fast enough to complete 3D simulations (even rudimentary ones) in any reasonable amount of calendar time; based on 2D runs that were done at this time, we estimated that 3D runs would take more than 60 years to complete on the CRAY architecture.
With the advent of the Department of Energy's ASCI (Advanced Simulation and Computing Initiative) program 1 in 1996, LANL and SAIC began a joint effort to convert the original CRAY vector Fortran-77 code to a new paradigm: parallel and modular code via effective software quality engineering. We chose to employ "vanilla" Fortran-90/95 with C-language interfaces for I/O and the explicit use of message passing to effect communication among the processors. Domain decomposition was used to spread the computational cells of the problem among the processors.
The massively-parallel version of the codes was designed to run extremely large problems in a reasonable amount of calendar time. Our target is scalable performance to 10,000 processors on a 1 billion cell problem that requires hundreds of variables per cell, multiple physics packages (e.g. radiation and hydrodynamics), and implicit matrix solves for each cycle. Currently, the largest simulations we do are three-dimensional, using over 600 million computation cells and running for literally months of calendar time using 2000 processors.
SAGE itself provided an early testbed for development of massive parallelism as part of the ASCI program and has been shown to scale well to thousands of processors [4].

SETUP and the RAGE Computational Grid
The RAGE computational grid is set up by specifying a coarse uniform mesh in the appropriate number of dimensions, D, and then refining that mesh according to various criteria. Each processor is required to have a multiple of 2 D zones at setup, which allows us to treat our base grid as if it too were composed of adapted cells (i.e. no special run-time coding is required for this coarsest level of cells). For efficiency, that multiple should be on the order of thousands, or circa 50k zones/pe. The RAGE computational grid is an octree decomposition of the model space (i.e. adaption bisects each zone in each dimension). The tree depth is selected to resolve field gradients within a material to one user-specified resolution and material boundaries to a second user-specified resolution. Because RAGE allows a mixture of materials in any cell, it is not necessary to subdivide the cells to the size of the smallest geometric features -if a zone is not mixed on cycle zero, it soon will be. We refer to the depth of the subdivision required to attain the requested resolution as the "level" of the grid. The 2:1 edge length ratio of mother to daughter cells is enforced as the maximum for any adjacent cells, whether or not they have a common ancestor.
The only geometric query that RAGE grid generation requires 2 is a "point containment" test, which determines if a specified region contains a query point. The user assigns a priority value for each geometric region to be gridded (normally, just the order in which regions are defined). When the grid generation algorithm queries to determine which geometric region contains a query point, each region is tested in order from highest to lowest (but one) priority. The first region for which the point containment test returns inside is returned as the containing region. This process assures that all assembly gaps and voids within the model are assigned a material, since region#1 covers the entire mesh. The problem domain is subdivided into as many subdomains as the number of processors assigned to the problem. The grid for each subdomain is then generated independently, except for two adjustments that require occasional interprocessor communication: 1) load balancing, and 2) ensuring the edge length ratio never exceeds 2:1.
With the mesh structure fixed, we approximate the material properties within the finest cells ("leaves") of the octree by using multiple queries in the spirit of a Monte Carlo integral over the distributions of density and temperature (or pressure and temperature, or density and internal energy -the three alternate ways to specify the initial state of a region).

Grid Generation via Geometric Queries
Since RAGE only requires a point containment test, a wide variety of geometric modeling techniques are feasible. Geometric support is currently divided among SAGE primitives, an integrated geometry query library, Spica, and a Lagrangian link file query library, Merak.

SAGE Primitives
The intrinsic geometry specification support includes spherical, ellipsoidal, conical frustum, triangular prism and quadrilateral prism regions, as well as perturbed boundary regions (sine and cosine perturbations), and radially swept 2D regions (used, e.g., for overlaying a 1-d spherical dump into a 2-d cylindrical mesh). Use of the intrinsic geometric modeling involves describing a physical problem within the input file and then verifying through RAGE execution that the desired model has been achieved (various graphics options allow the user to examine the "cycle zero" mesh and properties before starting the first physics cycle). These geometric primitives combined with judicious application of the RAGE region-priority ordering scheme permit a surprising number of complex simulations.

The Spica Library
The Spica library allows RAGE to generate grids from OSO geometric regions, OSO being a Los Alamos supported geometric modeling program [5]. Geometry created in (or imported into) OSO and used for initial grid generation, includes CSG combinations of triangular facetted closed b-rep surfaces (which in turn are imported into OSO as stereo lithography or "STL" files), solid representations, and implicit surfaces. OSO supports importing data from several geometric modeling and commercial CAD packages.
Spica currently accepts STL files as input, simply because many packages write that format; many other equivalent formats can easily be translated to STL. The solid geometry description nodes contain a variety of primitive objects: box, cone, cylinder, sphere, parallelepiped, torus, and a surface of revolution generated by a piecewise linear curve. These primitives are composed into regions using OSO. While it can be challenging to create extremely detailed models using these primitives, there is a considerable query speed advantage to doing so. OSO also supports STL files as first class primitives that can be translated, rotated, scaled, and combined with the previously specified primitives to compose regions. In practice, combinations of all the geometry description node types are used in computational models.
Separating the region selection from the geometry description is useful because it moves the exact definition of inside or outside to the tree, rather than leaving it up to the geometry generation package, with the advantage that all these CAD/CAM descriptions can be used in RAGE at the minimal cost of the query routine. Figure 2 shows an internal combustion engine model that was created in Pro/ENGINEER, output as separate STL files for each region of the assembly, and collected for viewing, scaling, positioning, and rendering by OSO. The OSO dataset was then processed via the Spica routines while reading the RAGE input file.   Figure 3 presents a close-up of a cross section of a RAGE grid generated from the Pro/ENGINEER-OSO model. The boundaries of the cylinder block, head, valve springs, and retainers are modeled at level-5. Note the gradual transition from level-1 to level-5 caused by the enforcement of 2:1 edge length ratios between adjacent cells. This grid has been restricted to five levels for purposes of illustration; typically, a well resolved computational grid will descend to level-8.
One and two-dimensional meshes can be set up from a cross-section of a 3-d model, and 3-d models can be created by rotating cylindrically symmetric 2-d models.

Data Structures and Parallel Implementation
In a parallel environment, RAGE is a MIMD (Multiple Instruction Multiple Data) code with synchronization at points in the instruction streams where inter-processor communication occurs. Such communication is handled at the lowest level by a machine-specific implementation of the MPI (Message Passing Interface) library, and at the application code level by a set of RAGE routines (the "token library") that sits atop the machine-specific library and hides most of the details of inter-processor communication from the application code. Although the code is mainly designed with a distributed-memory multiprocessor architecture in mind, it also runs on sharedmemory machines and indeed on uniprocessor machines, including Mac's and PC's.

Computational Cells and Cell Faces
Cells are labeled by a single index, l, regardless of dimension. The highest level of refinement may change as the simulation proceeds, and in addition, cells may be coarsened, i.e., a previously refined block of two, four or eight cells may be recombined into a single larger cell. Spatially adjacent cells are not allowed to differ by more than one level of refinement. If changing conditions, e.g., the passage of a shock front, require local refinement by more than one level, the refinement process is made to propagate outward to satisfy this requirement. Also, the refinement algorithm "looks ahead" so that it can, for example, anticipate the arrival of a moving shock front; this tends to put the most intense action into a locally uniform subgrid at the locally highest level of refinement.
At problem setup, blocks of two, four or eight level-1 cells are assigned index values in a linear sequence along a compact (Hilbert) space-filling curve [6], [7] that covers the computational domain, or (by default) in Fortran column-major order (we were surprised to find that on some architectures, this naive partitioning scaled better than the fancier space-filling curve). Cells are then allocated to processors by means of a simple partitioning of this curve. This produces compact spatial processor-subdomains and at the same time a linear global ordering of the cells. It also approximately minimizes the processor surface-to-volume ratio (i.e., the number of cells that share a face with a cell that resides on another processor divided by the number of cells that do not-roughly equivalent to the communication-to-computation ratio). Subsequent cell refinement, either at problem setup or later during the simulation, proceeds via the creation of new blocks of daughter cells. The newly created cells are inserted into the existing index sequence following the last cell in the block containing the appropriate mother cell; this preserves a degree of spatial locality in the cell ordering. It also allows us to enforce the rule that a cell and its siblings (cells that share the same mother cell) are always assigned sequential cell indices and that they always reside on the same processor -during the cycle the daughters are created; there is no guarantee they stay together on subsequent cycles after further adaption. The price that is paid for this locality-preserving indexing is that the cell index-set must be rebuilt at the end of every computational cycle, as part of the cell refinement process. Aside from the specific cases that have just been enumerated, there is no necessary connection between a cell's index and its spatial location in the mesh, so that a relatively large amount of indirect addressing is required.
Active or "leaf" cells are cells at the finest local level of refinement, i.e., cells with no daughters. The physics algorithms usually involve only the active cells, although the self-gravity package's Poisson solver does use inactive cell information to construct its multipole elements (see section 6). In addition to the cell arrays, the code constructs face arrays to be used in the implementation of the finite-volume physics algorithms; for example, diffusion coefficients are evaluated on cell faces. Faces are called active if they connect two active cells. If a cell face lies on a processor boundary, the face is assigned only to one of the two processors, in order to avoid double counting in code that loops over faces. Most cell arrays for physical quantities persist between computational cycles, except for indexing changes due to mesh refinement and load balancing operations, and are written out to restart dump files; face arrays, in contrast, are rebuilt anew at the beginning of each cycle.
Load balancing is not done every cycle, but only when the distribution of cells among processors has gotten out of balance by more than some preset threshold. The load balancing operation occurs via a simple repartition of the global cell index set, followed by the necessary movement of cell data among processors.

Communication
Communication in RAGE involves movement of data between adjacent cells, between cells and cell-faces, and between cells at different levels (mothers to daughters and daughters to mothers). In a multiprocessor environment, this will involve movement of data between processors. RAGE provides a small set of routines that execute these communication operations in such a way that the details of the inter-processor communication are hidden, as much as possible, from the application code.
Much of the inter-processor communication at the application code level employs the RAGE "clone" construct. This is a kind of ghost-cell mechanism, where clones are local (on-processor) copies of cells that reside on other processors and share a common face with a cell on the local processor.
The following example shows how the communication process looks from the perspective of the application code. In this example, we compute the divergence of a heat conduction flux, F = −D∇T , integrated over each active cell, using Gauss's theorem to convert the volume integral to a sum over faces of the normal component of the flux at each face. In the following pseudocode fragment, N (Numcell) is the number of active cells that reside on the local processor, N + (Numcell clone) is equal to N plus the number of cells cloned from adjacent processors, T(1 : N + ) is the temperature array on the local processor, and divF(1 : N + ) is an array to hold the computed result, i.e., the integrated divergence of the flux on the local processor. The arrays T and divF are dimensioned to have N + elements: elements (1 : N ) for local cells (those that reside on the local processor), and elements (N + 1 : N + ) for the clones (corresponding to cells on adjacent processors). The example contains three typical steps, namely, a gather step (clone get), a compute step (the loop over cell-faces), and a scatter step (clone put): call clone_get(T) ! Get off-processor temperature data ! i.e. T(Numcell+1:Numcell_clone) divF(1:Numcell_clone) = 0.0 ! Initialize all elements of divF do (for all active on-processor faces) ! in all directions ! llo = local index of cell or clone on the lo-side of face ! lhi = local index of cell or clone on the hi-side of face A face enddo call clone_put("add",divF) ! Scatter data to adjacent processors ! cumulate into on-processor elements divF(1:Numcell) = divF(1:Numcell)/Volume(1:Numcell) Here llo and lhi refer, respectively, to the cell or clone on the low-side and to the cell or clone on the high-side of the current cell-face in the loop over faces. The map between the active local cell faces and the adjacent active cell or clone indices llo and lhi is built anew at the beginning of each computational cycle, together with the associated machinery used by the two communication routines clone get and clone put. D f ace is the diffusion coefficient evaluated at the current cell face via an appropriate interpolation using data from the adjacent cells (e.g., a harmonic or geometric average) and A f ace is the area of the face. x lo and x hi are the sizes of the cells referred to by llo and lhi. Once the off-processor contributions have been cumulated into the on-processor data via the clone put, they are divided by the volume of the cell in order to ensure that divF has the correct dimension to be the divergence of a flux (∇ · F). The calculation of the gradient, by using half the sum of the widths of the cells (∆x llo + ∆x lhi ) allows this expression to be used without modification at the edges of a periodic boundary, whereas using the difference of the positions (x lhi − x llo ) of the two zones would generate nonsense in that case; otherwise, the two expressions would be equivalent.
The clone data structure is built by the token library routines at the beginning of the computational cycle and includes a list of the off-processor cells adjacent to the boundary of the local processor, a list of the on-processor cells adjacent to the local processor boundary, and buffers to hold corresponding data to be sent to and received from other processors.

Parallel I/O
When it is necessary to write (or read) restart dump files, RAGE can use either the MPI-IO or specialized bulk I/O libraries at the lowest level. 3 At the code developer's level, a subroutine is invoked to write either a scaler or vector which may be the same or different on different processors. Multidimensional mesh-arrays, e.g. cell momentum(1:numcell, 1:numdim), are written one mesh-array at a time with an outer loop(s) over the non-mesh index(es): call pio_io('same', 'write', 0, 'do_hydro', do_hydro, pio_err) do n = 1,numdim call pio_io('diff', 'read', n, 'vel', vel(:,n),numcell,pio_err) enddo The PIO package has both serial and parallel I/O capabilities. Serial I/O is the default, although a flag allows one to switch to MPI I/O or another parallel I/O package called bulkio which is based on asynchronous MPI calls that move data to a set of I/O processors. Both packages can be tuned for optimal performance using input file parameters, e.g. striping factors, configuration lists, etc. The serial I/O package uses MPI sends and receives to move the data from each processor to the I/O processor, overlapping the writing of one processor's data with the receiving of another's.
When parallel I/O is activated, the lengths of the array on the individual processors are gathered and cumulated so that each processor knows at what location in the dump file 4 to begin processing its contribution to the total. All of this data management and the calls to specific I/O library routines occur at a low level that the (physics) code developer never sees. In the case that the data is the 'same' or when serial I/O is selected, only the I/O processor (usually pe=0) does the actual work. Regardless of the type of I/O, a return flag (e.g. pio err, above) indicates whether or not the data has been successfully read into the array or written onto the disk.
Restart dump files are formatted to be valid Scientific Data Sets in the HDF conventional terminology; everything but character strings are written as IEEE 64 bit Reals and characters are written as ASCII character strings. On input, RAGE swaps the bytes on the fly when it detects that the "endian-ness" of the machine is different from the big or little "endian-ness" of the file it is reading, as well as converting Cray 64 bit Reals to IEEE 64 bit Reals when such dumps are encountered.

Parallel Scaling Performance
In running simple timing problems on various architectures (e.g. hydro-only or conduction-only on a fixed level mesh 5 ), we have found relatively good scaling with processor count, in the sense of constant work (fixed zones) per processor, as shown in figure 4. The metric used for performance is the processing rate in terms of the number of cell-cycles that are processed per second on each processor, and is the inverse of grind-time that is often referred to. Almost a decade of different systems are shown in figure 4, the oldest being the ASCI Red system installed at Sandia National Laboratory in 1996 and the newest, the Blue Gene/L system, installed at Lawrence Livermore National Laboratory in 2005. The performance over this time has increased by approximately a factor of 20 resulting from both improvements in processor speeds, and network speeds.
Ideal scaling performance is a straight horizontal line at the level of the performance of a single processor when a constant amount of work is mapped to each processor; very good scaling performance of SAGE is observed on these large-scale systems. The performance on almost all of these systems degrades only by about a factor of 2 at 1,000s of processors. The degradation results from the need for communication between logically neighboring processors as well as the need for collective operations. Note that the performance was measured on the largest sized systems installed except for the Blue Gene/L which will be upgraded from the 2K processors indicated in figure 4 to 128K processors during 2005 Further details on the performance of SAGE on the ASCI White (IBM), ASCI Blue Mountain (SGI), and T3E (Cray) systems can be found in [8]. Details on ASCI Q (Compaq) performance as well as how SAGE was actually used in the optimization of the system performance is described in [9]. The performance of Lightning (Linux) has been described in [10], and early details on the performance of Blue Gene/L (IBM) is given in [11].

Adaptive Mesh Refinement Overview
The adaptive mesh refinement (AMR) method used in RAGE was initially developed to solve multidimensional, time dependent shock hydrodynamics problems. The particular problem of interest at the time required high resolution at a shock front that could be kilometers away from other resolved features in the mesh. In choosing a method, it was desirable to have rectangular 5 The scaling problems are 3-d, consisting of three nested blocks of γ = 7/5 ideal gases at different temperatures and pressures (ρ = 0.1, 10, 0.1 g/cm 3 ; T = 10 3 , 0.1, 0.1 eV; p ∼ 10 13 , 10 11 , 10 9 erg/cm 3 , side length Lx,y,z = 30, 15, 6 cm respectively), using the same number of fixed zones on each processor. The same logic is executed in all zones whether the advected quantities will be zero or non-zero. The problem can be run in hydro only mode (as shown in Figure 4, or in a pure conduction mode, not shown but giving similar scaling results -each face calculates a conductivity, and adds elements to a row and column of the matrix for subsequent inversion). Because of indirect addressing, the work per zone is the same whether zones are adapted or not, and it would be extremely difficult to force the same amount of adaption on every processor for a scaling study.  (square) cells where integration methods with well understood convergence properties could be used.
A primary consideration was to ensure global conservation (of mass, momentum and energy); because the arbitrarily-orientated patch-based methods such as those of Berger and Oliger [12] make conservation difficult, we chose to have the underlying coarse cell boundaries coincide exactly with the finer mesh boundaries [13], which greatly simplifies the process of maintaining global conservation.
The "patch-based" AMR of Berger et al. is not, however, used in RAGE. Our "cell-based" AMR method is most similar to that of Young et al. [14] -it is a block-structured mesh with cell refinement limited to bisection in each dimension where adjacent cells differ by at most one refinement level.
RAGE uses a Continuous Adaptive Mesh Refinement (CAMR) algorithm; it is continuous in space and time in the sense that mesh refinement occurs on a cell-by-cell and cycle-by-cycle basis. The mesh is evaluated for refinement on every cycle throughout the entire domain. The overhead associated with this CAMR technique has been empirically determined to be about 20% of the total runtime, which is small compared to the gains of using CAMR instead of a uniform mesh.
The CAMR algorithm has two phases. The first or setup phase performs refinement at the start of a problem and the second or dynamic phase performs refinement during each timestep. It is convenient to separate the algorithm in this way to allow for different refinement criteria in each phase.
In the setup phase the base grid is refined by checking whether a cell in the grid is composed of a pure material or contains an interface. It is frequently desirable to allow material interfaces to be refined to a different maximum level than in the bulk interior of the material (e.g. during the passage of a shock). In the case of RAGE, material interfaces are determined either by changes in density or changes in material volume fraction as determined by the user. Additionally, if a cell and its neighbor cell differ significantly in density, pressure, or velocity component then refinement is performed on that cell. The user can also specify levels of refinement at predefined points or regions.
Once the grid has been generated and initialized, the governing equations are integrated one timestep and the entire mesh is evaluated for refinement before any data is moved. This evaluation is a two-part process and occurs on every cycle for every cell, not just the top level cells. The reason for this will be discussed shortly. The first part of the process determines which cells are to be flagged for refinement, coarsening, or neither. The second part of the process handles the addition of new cells, the removal of others, and their associated linkage to neighboring cells. In RAGE, activity tests are performed to refine the mesh ahead of shock waves. These activity tests use small changes in velocity and pressure to trigger refinement. Cells are also refined based on first-order estimates of truncation error in a Taylor series expansion of cell properties. Because error estimates for several cell properties are considered, the maximum error is saved for a cell. This method favors higher refinement in the interest of better accuracy.
When all cells have been evaluated, the second, data-movement, phase takes place. Cells are added and removed according to how they are flagged in the first phase. However, there are two rules that can override that phase: (i) no cell is allowed to be more than one level different from its neighbor (including diagonal neighbors), and (ii) cells will only be removed if it is unlikely that in the following cycle that cell would be refined. The last condition is accomplished by basing cell coarsening on the error estimate in the mother cell as well as the daughters; thus the need to update mother cells as mentioned above. The mother cells are used because they would become the top level cells on the next cycle if the cell was coarsened. Algorithmically it is easier to update all cells rather than restrict the update operation to top-level cells and their mother cells. The number of cells unnecessarily updated ("grand-mothers") is a small fraction of the total.
When cells are added to the mesh, they are added in groups of known quantity, therefore only a reference to the first daughter cell needs to be stored, since the remaining cells follow in memory. Inter-processor communications are minimized by inserting the newly created daughter cells immediately after the mother cell in the global cell list. Load leveling is only done on those cycles where an individual processor's cell count exceeds the average cell count by some tolerance (approximately 10 percent).

Hydro
The hydrodynamic package in RAGE uses a two-shock or single-intermediate-state approximate Riemann solver to calculate the impulse and work on each face of a cell, using an alternatingdirection-explicit (ADX) [15] technique. The Riemann solution is used to determine the trajectory in the donor zone along which a Lagrangian mass particle would move through the unperturbed state (and possibly the "contact" state), in order to just reach the face at the end of the timestep. This trajectory distance is then used to integrate all quantities that will be advected from the donor zone, and allows us to finesse some details in the typical Eulerian Riemann-solver methods.
For each direction, we test on overall "volume" changes within zones, re-evaluating the sound speeds where appropriate and iterating the Riemann solver. This semi-implicitization allows us to handle collapsing bubbles or strongly shocked zones even though the solver itself uses only a "weak shock" approximate Riemann solver. This goes beyond strong-shock or general Riemann solvers to allow a coupling from logically separated zones through their common neighbor, preventing unphysical results from occurring due to the independent updates. We will demonstrate that this method is 2 nd -order accurate in space and time on smooth problems and 1 st -order accurate on discontinuous or shock problems. In general, our philosophy has been to develop a method that has a demonstrable convergence rate under mesh refinement, and to then work at reducing the error at a given level of refinement.

Conservation Equations
The hydrodynamical equations solved in the hydro package are: where RAGE stores extensive arrays of mass M = V ρdV , momentum M U = V ρudV , and total energy M E m = V ρE m dV ≡ V ρ(e + 1 2 u 2 )dV , with u the specific momentum (velocity), and E m the specific total energy. The treatment of gravity will be deferred to section 6; for the purposes of this section, we will assume g = 0.
RAGE has a strength-of-materials package, which includes a simple failure model and deviatoric strain-dependent stresses (e.g. Steinberg-Guinan), but for this paper, we will only consider isotropic stresses, P = pI. RAGE also has an operator-split heat conduction package that solves ∂(ρEm) ∂t + ∇ · (−κ∇T m ) = 0, but this is infrequently used compared to the radiation diffusion package discussed in section 7.
The conservation equations can be converted into a Lagrangian form which will be used to calculate U and P at the half-timestep, for input into the Harten-Lax-van Leer (HLL) [16] approximate Riemann solver: Using Eqs. (4) and (6) and the thermodynamic identity, c 2 s = dp dρ e + p ρ 2 dp de ρ , we can also write: Aside from E m , and unless otherwise noted, in this section we will generally use the convention of upper case letters for Lagrangian-frame variables, and lower case for Eulerian-frame variables.

Boundary Conditions
The hydrodynamic package in RAGE has only two real boundary conditions available to the user: reflecting and outflow; it has a "treatment" to allow inflow.
The reflecting boundary condition (BC) is self-explanatory, and was the original BC in the code. The justification was that with AMR technology, one could use very large zones to put the boundaries far away from the region of interest, and adapt that region appropriately. All (normal) derivatives are set to zero in zones adjacent to a vacuum face, and all (normal) components of velocity are set to zero on such faces.
Since then, we have implemented outflow boundary conditions, individually specifiable on the high and low faces of the mesh in each direction. In this case, zonal values of derivatives normal to the face, instead of being the limited average of the gradient on the far face and zero at the vacuum face (i.e. therefore zero), are just the value calculated on the face of the zone not adjacent to vacuum. The velocity of fluid on this face is calculated by the Riemann solver and is not forced to zero, as with reflecting BC's.
In order to meet the demand for some type of in-flow boundary condition, we have developed a crude "freeze boundary" condition. In this case, the user specifies a geometric extent of the problem space in which all properties are to be frozen. At the top of the timestep, the list of zones falling in this space are recorded, along with the values of all quantities subject to advection (in those zones); the hydro then operates as usual, and at the end of the timestep, those zones are re-set to their beginning of timestep values. This has been adequate to handle the class of problems that have a constant inflow, and can be extended to allow simple temporal ramps (albeit only first order accurate in time).

Limiters and Linear Reconstructions
Although users have the choice of a number of different limiter methods with which to calculate gradients, the two most commonly used are the minmod and the modified van Leer. Slopes, s ± , are calculated within a zone by first linearly interpolating (between neighboring zones) a "facevalue" of the quantity on the faces normal to a given direction on the mesh. These face-values are then broadcast into cell-centered "hi" and "lo" arrays, at which point one can calculate the slope of the quantity between between the high side and the centroid of mass (s + ) or between the centroid and the low low side (s − ).
The minmod limiter generates the most conservative slopes for reconstruction, allowing monotonic discontinuities at a face between zones without introducing any subzonal saw-tooth patterns in the reconstructed field. This ensures that an approximate Riemann solver will satisfy the correct entropy conditions. The van Leer slope and limiter [17] can be defined in terms of a generalized minmod: where the average, s++s− 2 , is the 2 nd -order estimate of the slope between the three zones and van Leer's definition of this limiter has F ≡ 2 (F = 1 reduces this to the ordinary minmod). Our modified version uses F = 3 2 , which van Leer [17] has noted gives about the limiting power of his harmonic form, s harm = s+s− s++s− [(sign(s − ) + sign(s + )]. We have found that the use of the standard F = 2 in this limiter introduces a small amount of ringing behind some shock fronts in problems we have run, while using F = 3 2 does not produce any such visible anomalies. We have also found that the problems associated with stronger-than-minmod limiters (i.e. the saw-teeth in the reconstruction) is accentuated in two-dimensional problems, even if it is imperceptible in one dimension, and these problems are alleviated by use of F = 3 2 . Our goal is to have a spatially and temporally second order accurate hydrodynamic method, physics permitting. For such temporal accuracy, we use a half-step, whole-step method to advance the states for the Riemann solver. Second order spatial accuracy can be regarded as that accuracy that will exactly preserve linear solutions across the mesh, which means that a minmod limiter will converge to the same order as a van Leer limiter, although in general the van Leer offers improvements in accuracy at a given level of refinement.
Linear profiles are constructed in each zone with monotonically limited slopes with the zonal average located at the center of mass of the zone. For example, Kinetic energy is calculated with a 3-point Gaussian quadrature, so that the worst case, 5 thorder polynomial (in 1d spherical geometry), r 2 ρ(r)u 2 (r)dr, is integrated exactly, allowing us to always calculate an "exact" internal energy, The velocity, u(r), is determined by a two-step method. Beginning with the specific momentum, u 0 = M U/M , we calculate its gradient, ∇u 0 , and integrate an estimate of the momentum, We then recalculate the velocity to minimize the error between this estimate and the known value of momentum, We have not found it necessary to converge the velocity any further, but we do calculate ∇u for use later.

Riemann Solver
A single-intermediate-state approximate Riemann solution [16] takes values that have been reconstructed on the Left and Right sides of a cell's face (U L , P L , U R , P R ), and combines them into the intermediate state, (U * , P * ), shown for example in Fig. 5 lying between the lines U L − S L and U R + S R . This state is a weak (i.e. integral) solution to the momentum equation on both sides of a face (we use the "mass" coordinate, dm = ρdx = ρc s dt in the following formulae):  On the assumption that a signal propagates from the face at a "signal speed" of ρc into each cell, these two equations in two unknowns can be solved for the intermediate state: where the impedances, ρc, are the fastest signal in each zone [16], found by iterating Eqs. (9) together with the "L" and "R" versions of Eq. (10) to pressure convergence. We have not, however, found it necessary to perform such iterations, given the "anti-cavitation" logic described below in section 4.3.3, and use just the "weak-shock" solution for the intermediate state, solving Eqs. (8) and (9) with the zone-centered, beginning-of-timestep values of c L and c R . We do use P * in a single pass through Eq. (10) before our advection logic, because this (ρc) * can be shown (via Rankine-Hugoniot relations) to be equal to ρ 0 U s , where U s is the speed of the shock/rarefaction in the unshocked fluid (of density ρ 0 ). The characteristic labeled "U L −S L " in Fig. 5 (and Fig. 6) moves in the lab frame at a speed of U L −U s(L) (similarly, the one labeled "U R + S R " moves at U R + U s(R) ). We want to ensure that advection is done as accurately as possible, and use certain Lagrangian concepts to effect that end. Note that a Lagrangian test particle -in the absence of any gradients -will start at the position "−D" in Fig. 5, and will drift at a constant velocity (U L ) until it hits (or is hit by) the characteristic wave that sweeps past it, jerking it to a new drift velocity, U * , with which it will drift to the face of the zone at the end of the timestep. If we take our linearly reconstructed density profiles and integrate them between −D and 0, we will be able to advect the "exact" amount of mass that needs to be moved from the donor cell to the donee cell.
The direct-Eulerian method in RAGE is a hybrid that uses a Lagrangian Riemann solution to infer the correct Eulerian fluxes, in particular, a "manifestly obvious" prescription for the multi-material continuity equation's mass fluxes. Using a Riemann Solver in the solution of any set of (conservative) hyperbolic equations involves four steps [18]: • The calculation of interpolated profiles for the dependent variables, • The construction of time-centered left and right states, • The solution of the Riemann problem with those states to generate the evolved state, • the conservative differencing of the fluxes in the original equations.
Rather than calculate the left and right states by averaging at t n over the domain that will be causally connected to the face [18,19] during the timestep, we calculate point values on either side of the face: p L + 1 2 ∆x L ∇p L and p R − 1 2 ∆x R ∇p R . As in van Leer [20], it is understood that 1 2 ∆x is actually the distance from the face-center to the centroid of mass. For T-cell faces, it really is ∆x · ∇p. These (p, u) ± point values are then advanced to the half time using the Lagrangian expressions (Eqs. (5), (7)), so that the point value Riemann solutions calculated from them will be second-order accurate in time.
RAGE goes beyond the standard Riemann solver logic to make a further attempt to include gradients of velocity and pressure in its reconstruction of the amount of mass to flux across a face. Clearly, if there were no gradients, the jump in velocity across the leftward traveling characteristic in Fig. 5 would always be U * − U L . However, if there is a gradient in velocity, then our Lagrangian test particle's initial velocity will nominally be U L − Du x , where u x is the gradient of velocity at the beginning of the timestep; it is reasonable to also assume that the final-state velocity will not be exactly U * either. Thus, we regard the Riemann solution as giving us "point values" at the half-time rather than time-averaged values, and make the further assumption that the velocity jump across the left-bound characteristic in Fig. 5 is a variable, v j (f ), where f is the fraction of the timestep that elapses before the Lagrangian particle-path hits the characteristic. We define v j as: where we denote the gradient in velocity at the beginning of the timestep with u x = U x , and the gradient at the end of the timestep (i.e. the gradient of the Riemann solutions) with U * x . The characteristic Lagrangian speed, S L , is calculated after solving our approximate Riemann problem, with (ρc) * as in Eq. (10). The sign of S R is positive, meaning that all relevant characteristics (those in a donor cell) move at "U + S". These expressions for U (f ) in v j (f ) are not immediately obvious: in principle, the shock or the tail of a rarefaction moves at a speed of U * + S L , while the head moves at U L + c L . In order to be consistent with the single-intermediate state approximation, we assume that the distance from the contact to the tail is S L f ∆t (measured from the contact moving at U * ), and we assume that the distance from the face to the head of the rarefaction is also S L f ∆t (measured from the face in the unshocked fluid frame moving at U L ). Thus U * (f ∆t) and U L (f ∆t) at Eq. (11) and the jump across them.
If there were no gradients in the problem, then the characteristic in the donor cell will be moving at speed U L + S L , and the Lagrangian point will have been drifting from the coordinate, −D, with a velocity of U L . Their intersection at a given time gives a constraint, −D + U L f ∆t = (U L + S L )f ∆t. It will be noted that the actual velocity of the fluid drops out of this constraint to give, −D = S L f ∆t. When gradients exist, the characteristic, U L + S L , ought to show some curvature, but we assume that it is locally straight at the point of intersection, traveling at some unspecifiedẊ + S L . We have a second constraint from the definition of the total distance that Lagrangian point moves to get to the face. In general we have which simplifies to Solving the second equation for D, and plugging it into the first, we have an equation for f alone: Once this equation is solved for f (e.g. by Newton-Raphson iteration), the second equation, (15), can be immediately solved: Clearly, if we divide D by ∆t, we will have an average velocity, u, that should compare to what would have resulted from a more standard Eulerian direct solution, u n+1 . Building on that insight, we construct the analogous value of pressure, where the pressure jump across the characteristic, is the analog to v j (f ), with P * x the derivative of the Lagrangian Riemann pressures across the cell (the "point" values on the two faces divided by the size of the cell).
Consider now the case when U * > 0 and (U L + S L ) > 0, so that no characteristic reaches into the donor zone (Fig. 6). In this (transonic/supersonic) case, "D" is determined trivially from the expression for the straight path in space-time, where u L,x = du dx n L , and U n+1/2 L is the time-advanced value of the velocity on the left-side of the face.
In order to calculate the value of p for the impulse in this supersonic case, recall that any Lagrangian derivative can be written as an Eulerian derivative: The Eulerian value at fixed position at the end of a timestep is given by: Given that p x = P x to the order of accuracy we need, and since p n ≡ P n at the beginning of the timestep, where the second line follows from the definition at equation (18). The average velocity, u, could have been obtained by similar reasoning, and agrees with that from equation (20).

Impulse and Work Terms
As we have seen, the volume of fluid that is moved from the donor to donee zone is based on the integral from −D to 0 (for rightward flowing fluid). In slab geometry, this is A f ace v∆t, but in cylindrical and spherical geometry, we can calculate the ± r δ ∓ (r − D) δ (δ = 1, 2) exactly.
Using the average values of p and u on the face, we calculate the impulse at the face f i (face normal in the i direction): R + I f , and similarly, the work done on the face is given by, R + W f , (the area, A fi , has dimensions of solid angle in spherical geometry, length in cylindrical geometry, and area in slab geometry). These changes to mass, momentum and energy are temporarily stored as increments in case the anti-cavitation logic requires us to repeat the Riemann step with improved sound speeds. The state variables are updated before the next directional sweep in the alternating direction explicit (ADX) logic.

Artificial Viscosity
It has been noted that Godunov solvers in multi-dimensions have to be augmented by some form of artificial viscosity [21] or by a more dissipative low-order solver [22]. RAGE too has an artificial viscosity, specifically in order to symmetrize numerical errors.
In particular, we want symmetric (comparable) errors on problems involving shear, independent of the orientation of the shear. Consider a 2-d box with (e.g., 4) alternating bands of fluid (aligned with mesh boundaries) moving to the left and right with periodic boundary conditions (e.g., Kelvin-Helmholtz without the perturbations). Because the velocity field contains discontinuities, Strang splitting will no longer be second-order accurate in general; but given our direction-split solver, such bands decouple completely, giving zero numerical error in this orientation. Now consider the flow field rotated by 45 • : on each face that straddles the discontinuity, the advection logic on each face will calculate an average velocity of zero, resulting in no change to those border zones; but the intermediate state Riemann logic will find a non-background pressure (the two velocities are in opposite directions and limiters turn on) on the face, and that will flux some negative momentum into a positive momentum zone (or vice versa), thereby smearing the discontinuity. This smearing is due to the 1 st -order (slope limited) Riemann solution being Strang-split in this orientation.
Our approach to such rotational asymmetry is not to attempt to remove the error along the diagonal, but rather to add some error along the axes (but in such a way that in smooth flows, it will have zero magnitude and thus contribute no error). To that end, RAGE has developed a tensor-like artificial viscosity of the form for each direction, d = f , parallel to a face f . Thus the y and z (d) components of velocityjumps contribute to those components of momentum fluxed across a face normal to the x (f ) direction. (Note that in smooth flows, ∆v ||d → 0 because the (unlimited) slopes will now construct continuous velocities on each side of the face, turning off the viscosity.) The energy fluxed across the face f by these components is calculated by dotting this component of the tensor into the velocity parallel to the face, Figure 7 shows an example of the effect of such a viscosity on the vortical flow induced by a jet along along the horizontal axis, in this case for the double Mach reflection test problem [21] (960 x 240 mesh, ∆x = ∆y = 1/240, t = 0.2). This same viscosity usually eliminates the "carbuncle" growth [22] seen at shock fronts along the axes of a multi-dimensional problem.

Anti-Cavitation Iteration Logic
Consider the case that a single, almost empty zone is sandwiched between two high-pressure zones -for example, a collapsing bubble. Iterating a standard Riemann solver will ensure that no more than half the zone will fill with material from any face, but on the next cycle -if the fluid was relatively incompressible, like water, the density in that former bubble might become 1.01, in which case the over-pressure would cause the now over-dense zone to explode on the next timestep, sending an unphysical shock throughout the fluid.
To prevent this, each face records the effective (1-d) change of volume, u ∆t A f ace , which is added and subtracted from the appropriate zone-centered volumes, and a new volume is calculated: V new = V 2 old /(V old + ∆V ). If the "zone" shrinks sufficiently, updated densities and specific energies are used to invoke the equation of state for a new sound speed in those zones (and those zones only), and the Riemann solver is re-invoked (for the entire mesh), until the new volumes' changes are within some small tolerance.
Each such iteration requires a new inter-processor communication to pass the sound speed information, and the Riemann solver has to re-communicate the intermediate state values on the faces back to the zones on either side of the faces (for the f -finder logic). Otherwise, this outer iteration loop costs little more than a typical Riemann iteration loop. However, by doing the iteration here, we can avoid the over-pressures that occur from causally disconnected faces updating their common zone with too much material.
van Leer [20] found that the weak solution is more than adequate for rarefactions, and Colella [23] found that only 1-2 iterations were necessary in a large class of problems he examined. In those few cases, we can also handle strong shocks with our beginning of timestep sound speeds in the weak-shock Riemann solver, since the few zones where a shock occurs will force one of these outer iterations with the shock-modified sound speeds in those zones. (The biggest problem we have found occurs when a very strong rarefaction would straddle the face if we allowed it. Then our advection/donor logic will make a large first order error, as seen in Fig. 13 on LeBlanc's shock tube problem.)

Advection
When our problems involve more than one material, it is not enough to just integrate the linearized fractional mass-densities in order to update the equations (e.g. ∆M j = D 0 ρ j r δ dr). If we wish to ensure that materials in equilibrium remain in equilibrium, we have to advect a certain volume of each material (based on the linearized gradient of volume fraction), multiplied by the partial mass density or partial internal energy density -for the continuity or energy equation, respectively. This method works because the partial energies as well as partial volumes of materials are returned by the Multi-Material EOS package (which finds the partial volume that each material has to occupy in order that the weighted partial densities sum to the total density, and the weighted fractional energies sum to the total energy of the zone; see section 5, below).
Because RAGE is a single-fluid code, all specific momenta are identical, and nothing fancy needs be done to update the momentum equation -it suffices to integrate the bulk momentum profile.
The original RAGE philosophy was to minimize the complexity of the various physics packages in the code. In the context of advection, this meant using the same slope-limiter logic to advect fractional masses that was used in the hydro to advect momentum and energy; if sharper "boundaries" were desired, one should increase the level of adaption, rather than impose a "subzonal" model of an interface that wasn't part of the equations being solved. In the context of an ASCI production code, not all problems run could afford to refine the mesh to the degree required -even ASCI computers have finite memories -and the Interface Preserver technique was developed to sharpen contact surfaces without further adaption. In a limited number of cases, it was found that even this is not sharp enough, and work is on-going at LANL to implement a Volume of Fluid method [24] that will confine mixed cells to a single interfacial zone.

The Interface Preserver (IP) Method
Having taken into consideration the existing data structures and approach to adaptive mesh refinement in RAGE, an efficient algorithm has been chosen to preserve the steep slope or profile of a given volume fraction. The approach makes use of the compact stencil currently used in the code. Compact means that only those mesh cells directly adjacent to the mesh cell at which a slope is being calculated are used in the reconstruction. Early renditions of this technique were developed in the late seventies and early eighties by Harten [25] and subsequently modified by Yang [26]. Although these methods are commonly termed slope steepeners, they were originally called artificial compression methods. The motivation was purely for contact discontinuities.
Here we have adopted the approach of Yang with modifications to capture material interfaces. For ease and clarity, consider volume fraction φ distributed on a one dimensional Cartesian mesh. If the volume fraction φ i is to be reconstructed using the IP method, an initial slope S i is always evaluated using minmod limiting (see section 4.2), a standard option currently available in RAGE. With the slope so determined, left and right face values can be calculated by where φ L i and φ R i are the left and right interpolated cell face values of φ i . To modify the initial minmod slope, we first determine the difference in the face values of φ on either side of the interfaces, and the change in slope is then evaluated as The final, but not necessarily monotone slope, S * , is then In Eq. (21), α i is a "discontinuity" or "material detector" that takes on values between 0 and 1 and is evaluated as This "detector" is effective in the following manner. At a material interface there will be a variation in volume fraction over a small number of cells giving a value of α i of order 1 while away from material interfaces, such as within a pure material or where material fractions vary over a larger number of mesh cells, its value will be near or equal to zero. The idea at this point is that this detector will not allow a compressing of materials in regions where material has been either intentionally mixed or the interface is no longer resolvable by the current mesh spacing. The factor of 33 in front of the discontinuity detector is a somewhat empirical constant controlling the onset of the preserver. Notice that if α i is 1.0, then ∆S is rather large (due to this factor of 33) and most certainly makes the final slope non-monotone. This aspect requires a final adjustment to S * to ensure that it produces no new extrema relative to the surrounding data. This operation is always performed at the end of the slope calculation. This final limiting gives our final slope S i from S * i . Although the presentation here is for the simplest of computational grids the actual implementation in RAGE does work with adaptive mesh refinement. Currently the IP method is used on volume fractions for all materials, all of the time, with α i being the sole "switch" for turning it on or off.

Implementation and Integration of IP into RAGE
The implementation into RAGE is more complicated than the above presentation due to the need to contend with variable cell sizes, the existence of T-cells involved with the AMR, and the use of one and two-dimensional meshes with non-Cartesian symmetry. However relying on the coding already in place has allowed for a relatively straight forward implementation of the method.
It turns out that for many applications outside of the area of ideal gases, the IP feeds back detrimentally on the hydrodynamics and yields unphysical results. We have implemented a technique where the IP is turned off when a shock is passing through a cell, it being preferable to drop the hydro variables to truly first order to alleviate some of the problems. Unfortunately, it has also been found that turning off the IP even for short amounts of time will allow enough material diffusion that the material detector in the IP will no longer see the material interface as sharp and never come back on even after the shock has passed. A first-order error can be seen to advect in from the pseudo-inflow boundary condition.
Future work will consider additional control switches. The IP is appropriate, in particular, for interfaces undergoing pure compression or expansion. However, interfaces where vorticity is large indicate an amount of physical mixing in which its use would be inappropriate. Thus an additional control switch may need to key off the level of vorticity within a given cell. Also, it has been suggested that as soon as material strength is turned off for a material that IP be turned off as well.

Accuracy and Limitations of the Hydro Algorithm
RAGE is second-order accurate in space and time when used to calculate "smooth" problems (and first-order accurate on discontinuous or shocked problems). We will demonstrate this by running problems on fixed uniform meshes with fixed timesteps, halving and quartering the zoning and timesteps in order to ensure known conditions (the adaption algorithm is not designed to be self-similar or scale invariant, so refined adapted meshes might not have exactly twice as many zones as coarse adapted meshes).
Godunov [27] has shown that conservation of entropy can be obtained by a linear combination of the mass, momentum and energy conservation equations, which means that since RAGE conserves mass, momentum and energy, it also conserves entropy. In order to demonstrate this, we refer to an adiabatic compression test problem (Fig. 8), where u(r) = −r, (ρ 0 = 1, T 0 ≈ 0), and show that the L 1 error of density hovers around machine accuracy independent of mesh resolution (until a first-order error propagates in from the outer boundary's pseudoinflow condition). Surprisingly, the internal energy (not shown, not surprisingly) is first-order accurate under mesh refinement. Figures 9 and 10 show RAGE results on another smooth problem, a "wave in a bucket", where a sinusoidal 10 −6 density perturbation is applied to a unit density fluid initially at rest. In Figure 9 the zone size and timestep were both halved for each curve; figure 10 shows why we usually halve zone size and timestep simultaneously -the coefficient of the ∆x 2 error term is much larger than that of the ∆t 2 or ∆x ∆t terms -it save us a few runs. Not shown are similar results obtained with 2-d waves in a square bucket, which are also second-order accurate in space and time -given our Strang splitting and the one dimensional results, this is to be expected. Figures 11 and 12 show results on a weak (Sod) shock tube problem, and Figs.13 and 14 show results on LeBlanc's strong shock tube problem, both of which are demonstrably first-order accurate as expected.
The order of accuracy is determined by reference to the L 1 error norms shown in Figs. 8-14. If we assume that the theoretical value is given by the computational value plus a truncation error, θ = C(∆x) + α∆x r + · · · , then the ratio of L 1 errors on two different meshes allows us to determine r: |θ − C(∆x)|/|θ − C(∆x/2)| = 2 r . The upper portions of Figs. 11, 12, and 14 show the result of calculating r, for adjacent sets of L 1 errors in the lower parts of the figures (the cycle-0 results have no error, and r = 0). In the problems that involve shocks, we have seen that this order of accuracy, or convergence ratio, is closer to three quarters than one, for reasons we will suggest presently. The same method applied to smooth shockless problems does generate convergence ratios very close to 2 (e.g. 1.95-1.98), when applied to the point-wise and mesh integrated results of Figs. 9 and 10.
As can be seen in Fig. 10, the L 1 error is sometimes unchanged when ∆t is refined at a given ∆x. For this reason, we quote convergence ratios based on runs where both ∆t and ∆x are simultaneously halved. Then the calculation of r is insensitive to whether the leading error term is α∆x r , β∆x r/2 ∆t r/2 or γ∆t r . 4.5.1 Why Shock tube convergence rates are less than Unity.
The monotonicity constraints of a finite-difference scheme will smear a sharp shock over 3 or 4 zones, no matter what the spatial zone size, so that those zones will always contribute about the same discrepancy, |X i − C i | 0, (X=exact, C=Code, i ∈ (3 − 4 zones) no matter what the zone size. If this dominates the error, then on mesh refinement, one expects that the new error will be half as much, resulting in a convergence ratio ≈ 1.
However, it does not seem to be commonly known that the typical order of accuracy in the treatment of a contact is O( n n+1 ) [28], where n is the order of the difference scheme. This means that a contact in a shockless (otherwise 2nd-order accurate) problem would be 2/3 order accurate, while a contact (in an otherwise 1st-order accurate) shock tube problem, would be 1/2 order accurate. While this doesn't dominate the overall error, it does drag our shock tube convergence rates down to approximately 3/4. Similar results are seen for MUSCL and WENO schemes involving shocks [29,  We have run a further problem, a 1-d Mach 10 shock driven into resting material (based on the 2-d double Mach reflection problem of Woodward and Colella, [21]) to show that a shock without any contacts really is first order accurate in space and time. The table inside Figure 15 shows that the L 1 errors in both density and velocity drop almost exactly by factors of 2 under mesh and timestep refinement; the figure itself shows the L 1 errors in density (multiplied by 10) compared to the highest resolution's density profile. The same error noted by Woodward & Colella, caused by specifying the time-dependent boundary conditions based on an infinitely  sharp shock instead of a numerically broadened one, show up here as two dips in the density behind the shock front. 6 We present in Figure 16 Figure 8, top 2 plots], although the last is a PPM (piecewise parabolic) method with 4 th order steepening, and our scheme is closer to PLM (piecewise linear method) with 2 nd order steepening). In this vu-graph norm, one can see the sharpening of contours by a factor of two as expected, although the magnitude of the error is such that, for example, Woodward and Colella's ∆x = 1/60 resolution calculation looks closer to our ∆x = 1/120 resolution.

Grid Imprinting and Grid Seeding
It may have been remarked (e.g., see Figure 1) that the square zones of RAGE do not provide a perfect representation of circular or spherical objects -the zoning results in a sawtooth representation of what should be a smooth curve (or inclined plane). There are a number of different cases in which this could present a problem, and we briefly discuss them. In the case of smooth  shockless flow, we have established that the Strang splitting results in a 2nd-order accurate hydro scheme; thus, while the grid may "seed" some fluctuations in volume-fractions around the circumference, the overall motion will be as correct as the hydrodynamics, and Rayleigh-Taylor perturbations in 2-d slabs, for example, do grow at the correct rates. In the case of shocked flow, there are two cases -a round shell may be converging (as in an ICF capsule implosion), in which case no problem arises with a shock going through such a shell and converging at the origin; however, in the case that a shock reflects off the origin and re-shocks the shell, the shell will exhibit significant instability growth due to the originally grid-seeded perturbations and the fact that gradients of pressure and density are oppositely directed (Rayleigh-Taylor unstable). This growth can be reduced by resolving the surface -at some cost -with finer zoning. If this re-shock occurs within a small number of zones of the origin (i.e., < 30 − 50 zones), then the shock itself will not have achieved self-similarity (roundness) after its bounce, and will hit the (constant radius) shell at different times at different angles. This can be ameliorated by resolving the section of the mesh at the origin as well. In the case that an outbound shock overtakes an outbound shell, there is some Richtmyer-Meshkov instability growth, but it is not as severe a problem since the gradients of pressure and density now tend to be aligned and this configuration is unlikely to be reshocked (or reshocked very strongly) by reflected shocks from larger radii.

Multi-Material Equations of State and Mixed Cells
In any Eulerian code, regardless of the degree of adaption, the majority of cells will eventually become mixed with more than one material. Given that observation, RAGE has developed a multi-material equation of state package (MMEOS) that assumes that all cells are mixed zones.
Given RAGE's time-explicit hydrodynamic algorithm, sound waves will cross the smallest zone in two to three timesteps, relaxing any subzonal pressure differentials in that time. While it    is not possible to hand-wave temperature equilibrium in as cavalier a fashion, RAGE nonetheless makes the assumption that all mixed cells are in pressure and (material) temperature equilibrium. This assumption has the advantage that there is a unique solution to the question of the mixedcell temperature and pressure, and has the further advantage of being consistent with the idea that if adaption (based on pressure gradients) has stopped in a region of zones, there is nothing further to resolve.

Inverted Pressure-Temperature Tables
In order to quickly invert mixed-cell equations of state, a RAGE preprocessor inverts the equations of state for every pure material before beginning a calculation. RAGE can do this preprocessing for (ρ, T )-based SESAME tables [32] , or a large variety of analytic (either (ρ, T )or (ρ, e)-based) equations of state (e.g. Mie -Grüneisen [33], JWL [34]), at which point only the inverted (P, T ) tables need be saved for use in multiple executions of RAGE problems.
The philosophy of the multi-material equation of state routine(s) is to take these inverted tables -v m (P, T ) and e m (P, T ), where v m ≡ ρ m −1 is the specific volume and e m the specific energy of a "chunk" of material m -and construct the total volume and energy by weighting the individual tables by means of the known fractional masses. In addition to the fractional masses, M m , and their sum, M , the volume of a zone, V , and the total energy, E, in a zone are also known quantities.
In outline, the user inputs to the preprocessor a range of temperature and pressure on which to build each material's inverted v m (P, T ) and e m (P, T ) tables. This is done by making Maxwell constructions in the original form of the tables, if they are not already there, and searching along isotherms (in (ρ, T )-based tables) to find the desired pressure. The values of 1/ρ and e are then saved. In order for tables to be invertible, it is necessary that various thermodynamic derivatives be positive: and these are enforced in the preprocessing step that creates the tables. When RAGE calls the equation of state to find a pressure, it knows the total volume V , mass M , internal energy E tot , and mass fractions, x m ≡ Mm M , within each cell. MMEOS begins an iteration process by searching for two adjacent isotherms that bound the specific total internal energy of the cell corresponding to the lowest pressure in the table and calculating the fractional distance between them that would give the correct total energy, such that x m e m (T hi , P lo ) , Using the weighted values of various thermodynamic derivatives, MMEOS also calculates so that this can be used to estimate the next pressure in the table at which to repeat the process. In general, one interpolates along an isotherm at the desired pressure to calculate the internal energy, until one finds the two bounding isotherms. Although there are various limits designed to ensure that one does not shoot too far through the table too fast, this in broad outline describes the method of finding the point at which a mixed cell exists in pressure and temperature equilibrium.

Speed and Accuracy
The procedure outlined above is quite fast, since the expensive table inversion process occurred once during a pre-processing step. Given the inverted tables, the time to look up a mixed equation of state compared to a pure equation of state is a linear function of the number of materials, where various table properties have to be mass-weighted (such loops would be short-vector or superscalar loops).
Since some other methods (e.g., pairwise balancing of pressures) scale quadratically with materials, they would be significantly inferior applied to some of the ASCI-scale problems (e.g., 50 million zones) that have upwards of 50 materials in the problem. Even so, this linear algorithm, on small problems with few materials ( 50k zones, 5 materials) can consume half the execution time for a pure hydro problem.
The accuracy of this process does depend on the models used in the original tables; if some SESAME tables have limited domains, then the extrapolations off those tables to construct the global inverted table will be problematic: densities on the lowest isobar, for example, are overridden so that they will extrapolate to zero at zero pressure for any temperature, while energies are reset to the next higher isobar's value so that de dp | T = 0. Similarly, de dT | p > 0 is ensured by requiring that e (T n−1 , P ) = min(0.999 · e(T n , P ), e(T n−1 , P )), starting at the penultimate temperature and working downward along each isobar.

MMEOS Limitations
The MMEOS algorithm was designed in the context of solid mechanics problems, and works equally well in the higher-temperature (LTE) regime at Los Alamos. 7 This does mean, however, that the algorithm is not easily adapted to the "2-T" regime of separate ion and electron temperatures -one would have to construct 3-dimensional arrays, e.g. v m (T i , T e , P tot ), assuming that it is the total pressure of each chunk of material that equilibrates. Moreover, one would need 4 such arrays: v m , e ion m , e elec m and g elec m ≡ p elec m (v m , T e )/P tot , in order to satisfy all the energy and volume constraints, as well as to partition pdV work among the ions and electrons. This is a subject of research at present.

Gravity
Because gravity is a long range force, the large spatial scales characteristic of astrophysical systems make necessary its inclusion in many problems. Three different types of gravity have been incorporated into the RAGE code: constant acceleration, analytic, and self-gravity. They may be used in any combination. The simplest is a constant acceleration, which may be applied in any coordinate direction. The magnitude and direction of the acceleration are specified by the user. The analytic gravity routines accessible through user inputs are spatially varying but constant in time. They include single and binary star potentials. Also included is a power-law mass distribution to model the potential due to a star cluster. The gravity routines are written to allow users to add other analytic potentials with little effort.

Self-Gravity
For problems involving self-gravity, the potential and acceleration are determined through a global consideration of masses. For a matter distribution characterized by a spatial function of density ρ(x, y, z), the specific gravitational potential (φ) is determined everywhere in space by solving Poisson's equation, where G is the gravitational constant, G = 6.67259 × 10 −8 cm 3 g −1 s −2 . The acceleration is in turn given by the negative gradient of the potential, g = −∇φ. The equations of hydrodynamics are hyperbolic equations and can be solved by a local (explicit) method, but Poisson's equation is an elliptic equation which must be solved by a global (implicit) method. For this reason, the calculation of gravitational energies and accelerations typically requires as much computing time as all other hydrodynamic processes combined. In order to model self-gravitating systems with sufficient spatial resolution within reasonable run times, a fast gravity solver is desired.
While many methods for solving Poisson's equation exist, the adaptive grid algorithm used by RAGE restricts the possibilities. The chosen method is the multipole method of Salmon, Warren, & Winckelmans [35]. The SWW method was developed for N-body particle codes, but the hierarchical tree structure employed is also compatible with the unstructured mesh of RAGE. At any moment in time, the discrete integral representation of the Poisson equation is written as This is a direct sum over q cells, each of mass m q , in order to determine the potential at a chosen location x. To solve equation (22)  The key to the multipole method lies in combining distant cells into groups and then solving for the potential due to the group rather than each individual cell. As SWW point out, this is analogous to treating the Earth as one extended object rather than a collection of individual atoms when calculating the gravitational force beyond its surface. The right hand side of equation (22) is approximated by a multipole expansion for a group of q cells within the region S: In this expression, M S is the total mass contained in S, x cm the center of mass of S, and Q ij the quadrupole moment of S about its center of mass. The dipole term of the expansion vanishes when taken about the center of mass. The quadrupole moment Q ij , is determined by where δ ij is the Kronecker delta function, and repeated indices indicate a summation. A multipole acceptability criterion (MAC) is required in order to quantify the error introduced through the use of equation (23). The multipole expansion requires the use of data on all cells, but RAGE typically computes data on leaf cells only. (A leaf cell is one that is not sub-divided into daughter cells.) So before the gravity solver can begin, the mass, center of mass, and quadrupole moment for all mother cells must be determined through the aggregation of data from the daughter cells.
Once the needed data are in place, cell-centered values of φ and −∇φ are determined for each leaf cell with a treewalk through the grid hierarchy. The treewalk starts with a loop over the cells on the coarsest grid level. If the cell in question is a leaf cell, its contribution to the potential is computed from equation (22). If the cell is a mother cell, its MAC is determined. Should the MAC prove to be satisfied, equation (23) is used to calculate the contribution to the potential. If the approximation error exceeds the allowable value, the source cell is divided into its daughters and the test is repeated. This continues until either the MAC is satisfied or only leaf cells remain. Note that in the case that AMR is not used, this procedure becomes the direct method embodied by equation (22). In the manner of SWW, we consider the MAC to be satisfied should the MAC error be lower than a user specified limit e: In this expression, p represents the highest order retained in the multipole expansion, d = x − x cm , and b = max q∈S x q − x cm . The moments B n are defined as Embodied in the MAC is the notion that the multipole approximation is more accurate if the distance to the measurement point is large (large d), the sources are scattered over a small region (small b), higher order approximations are used (large p), and the multipoles which have been neglected are small (small B p+1 ). For the RAGE implementation, the quadrupole approximation is used, so p = 2. As expressed by equation (25), the MAC is an absolute criterion. The error tolerance must be set according to the smallest contribution of ∇φ. However, applying the most stringent error tolerance everywhere on the computational grid will result in a higher accuracy than needed and a significant slow down of the code. The solution to this dilemma is to turn the MAC into a relative error criterion. This is accomplished by noting that the expression in equation (25) has the dimensions of a mass divided by distance squared. Therefore, a dimensionless relative MAC is defined by dividing the left hand side of equation (25) by M s /d 2 to get (for p = 2), This new relative MAC (named mac tol in RAGE) represents the relative size of the error in using the quadratic approximation instead of the direct sum. Typically a value of mac tol=.01 proves sufficient. It is useful to recast the relative MAC into a form in which it becomes a property of the source cell [36]. The MAC can then be calculated for all mother cells prior to the treewalk. A critical radius r c is defined by solving equation (27) for d. This problem is simplified by noting that by the definition of the moments in equation (26), B 4 is always greater than zero. Therefore, that term can be dropped from equation (27) without violating the inequality. The equation can then be re-written to produce a cubic expression for d: When treated as an equation, this expression can be solved analytically. The result is the critical radius r c , defined as the closest distance to the source for which the quadrupole approximation meets the criterion set by mac tol. Note that the computational simplicity gained by dropping the B 4 term comes at the price of a critical radius which is larger than is strictly necessary.

Integration of Gravity into RAGE
Once the gravitational potential and acceleration have been determined, the results must be incorporated into RAGE. The conservation equations for mass, momentum, and energy take the modified form: where φ is the gravitational potential, and ∇φ ≡ g the gravitational acceleration. Adding gravity into a Godunov code is a multi-step process. During each direction-sweep through a hydro cycle, gravity affects the dynamics both before and after the Riemann solver. Before being passed to the Riemann solver, the gravitational acceleration is combined with the acceleration due to the pressure gradient to form half-time values of u L and u R . After the Riemann solution is formed and the state variables are advected, the energy and momentum are then modified by the changes in gravitational potential and acceleration such that for a cell of mass m with a change in mass due to advection of ∆m.

Limitations
A significant decrease in the run time of the self-gravity solver can be achieved by taking advantage of the manner in which RAGE indexes cells among different grid levels. Setting up the best combination of grid size, AMR levels, and gravity input parameters requires some vigilance on the part of the user. As of this writing, the RAGE self-gravity solver operates only for three-dimensional problems with non-periodic boundary conditions. Two-dimensional solutions and periodic boundary conditions are under development.

Radiation-Matter Energy Exchange in RAGE
Many authors [37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60] have addressed the topic of handling the coupling of energy and radiation at high temperatures. Our approach has three interesting features: 1) It produces an exponential relaxation of the difference between these two, which is more accurate in general. We explain this in detail and compare it to alternative treatments. 2) It can handle arbitrary variation of C V and "Planck" opacity within a timestep, and 3) it uses a novel technique to solve our resultant set of nonlinear equations.

Nomenclature
The primary dependent variables in RAGE are E, ρ, E m , and u, where E represents the radiation energy density, ρ the mass density, ρE m the material total energy density, and u the material velocity. The material specific internal energy is defined by The temperature θ, pressure p, and specific heat C V are given in terms of e and ρ by an (inverted) equation of state, which may be tabular or analytic. The radiation coupling µ = ρκ abs and mean free path λ = 1/ρκ tot are also given by opacity tables or analytic functions. The quantities Φ and T R will appear; they are given by where a is the radiation constant. T R is called the radiation temperature, but we do not mean to imply by this that the radiation has any particular spectrum; T R is just a measure of radiation energy density.

Energy Equations: General
Although RAGE has options that enable more ambitious physics treatments, here we assume local thermodynamic equilibrium, with electrons and ions each described locally by the same temperature θ. RAGE is a multifrequency code, but this article will treat radiation in the gray diffusion approximation, for which only its energy density E is of concern. Its spectrum is of interest only insomuch as it affects the calculation of µ and λ; it is not calculated by the equations described herein. The spectrum may be the result of a simultaneous multifrequency-gray calculation (as described by Winslow and others [53,61,62,63]) of which these equations form the gray pass. Or it may be assumed to be that of a Planckian at temperature θ, or even at temperature T R . The treatment described herein applies to all of these cases, with small suitable modifications.
RAGE is not a relativistic code, but it does consider terms in u/c to first order [53,54]. All the quantities involved with radiation are taken to be evaluated in the rest frame of whatever fluid element they reside in at the moment, the co-moving frame. To O(u/c), there is a simple transformation between the radiant fields in a comoving frame and those of an inertial frame [53]. Finally, by making the P 1 approximation (I = 1 4π (cE + 3Ω · F) or equivalently, p r = E/3), the hierarchy of moments of the transport equation close at the momentum equation, so that the relativistically correct equations (to O(u/c)) become

Gray Radiation Hydrodynamics
We now make the diffusion approximation, ∂F ∂t = 0, relieving us of the burden of solving the radiation momentum equation, (36), and dropping the last term from the radiation energy equation, (35), above. We further assume that velocity-dependent terms do not contribute to the steady-state flux so that we can maintain F ∝ ∇p r . The resulting radiation diffusion equations for the energy E and flux F are where the diffusion coefficient D is given by and L denotes a flux limiter, a common device [55,56,57,58,64] invoked to limit the diffusion flux so that it does not exceed the streaming limit, cE. In RAGE this limiter, if active, depends on the dimensionless quantity λ|∇(E)|/E, among other things. RAGE uses the Levermore-Pomraning flux limiter [64]; others have been tried at various times (by us and others [57]) and not found to make much difference. The material momentum equation, to O(u/c), is and the material energy equation, to the same order, is with the material heat flux F m , given by The material conductivity, χ, can be read from SESAME tables or calculated from simple analytic expressions.

Split Equations
At each timestep RAGE splits the time integration into three separately performed phases: hydro, material heat conduction, and radiation. The various relativistic terms are treated differently in RAGE. The advection term, line 2 of Eq. (37), is calculated explicitly in the hydro phase and the updated result becomes the initial value E − , for the radiation phase. The energy and momentum changes wrought by the radiation pressure and the flux, line 3, are explicitly treated at the beginning of the radiation phase. Line 1 is of course solved implicitly as the final step of the radiation phase. Since u and ρ are regarded as time-independent functions of position during the radiation phase, the matter equation during this part of the operator split simplifies, leading to the following coupled equations, where material properties (D, µ, C V , etc.) depend on t because of their dependence on θ, which itself depends on the time dependent internal energy e. Of course these properties depend on ρ too; we shall not bother to indicate this. Adding Eqs. (41) and (42) gives the basic diffusion equation to be solved in this phase:

Special Case: Equilibrium Diffusion
If the coupling coefficient µ is very large, in a sense which we will make precise below, Eq. (42) becomes stiff and implies From the definition of Φ, Eq. (34), we have and the equation of state, E, provides the material specific energy e, and specific heat C V : We see that Eq. (43) requires ∂(ρe)/∂t, which we can evaluate as where the last line follows from the result of the large µ, Eq. (44) (remember that ρ is fixed in this phase). Inserting this into Eq. (43) produces the equation of equilibrium diffusion (This is often expressed as an equation for θ: (ρC V + 4aθ 3 ) ∂θ ∂t = ∇4aLDθ 3 ∇θ, and may be referred to as 1-T radiation conduction.) In general, (47) is an implicit and nonlinear equation. (The exception is the unphysical case in which D is constant and C V is proportional to θ 3 , introduced by Pomraning [65] specifically for the purpose of generating a simple test problem and so used by Su and Olson [59,66], and others [58,60,67].) It is implicit because the material properties C V and D are related to θ through the equation of state and opacity function, and the flux limiter L (if active) is explicitly nonlinear. It is nonlinear in addition because θ 3 appears, and θ is a function of E, not just a parameter: θ = 4 E a .

Special Case: Linear Diffusion
If µ is very small, the matter and radiation decouple, Eq. (41) becomes and the material temperature never changes. Unless the flux limiter or other phenomenon (e.g., two-temperature opacities, see appendix B) causes LD to vary with E, this is a linear equation and can be solved by standard methods.

Time Differencing
We now return to the case of non-equilibrium diffusion. One of the unique features of RAGE's treatment is the asymmetric temporal treatment of the radiation energy equation and the matter energy equation. As we will see, the radiation energy (diffusion) equation is handled with a typical backward difference approach, but the matter equation is rearranged, treated as an ordinary differential equation (ODE), and integrated exactly, leading to what we term exponential differencing.

Solving for E Using Backward Differencing
At every timestep RAGE integrates the diffusion equation from a starting time t − to a finishing time t + , an interval ∆t. It assumes that E(t) throughout the timestep can be approximated by its final value, E + . This is the backward difference approximation, only first-order accurate in time but stable and very robust. Thus integration of Eq. (43) over a timestep results in where LD is and the evaluation of LD is left ambiguous at this point. This is not a closed system, because we do not yet have a value for e + , but with the E variation now fixed the matter energy equation, (42), becomes or equivalently, Integration of Eq. (50) will enable us to calculate e + needed in (49), effectively allowing us to solve those two equations simultaneously.

Solving for e + Using Backward Differencing
A very popular approach to finding e + is to use a backward difference scheme for the matter equation also. This approach is used by Freeman [68], Knoll et al. [69], Shestakov et al. [58], and many others. It leads to a replacement of (50) by Eqs. (49) and (51) form a closed system once LD and µ are specified. It is implicit and non-linear for all the same reasons as the equilibrium diffusion system, Eq. (47), and in addition involves the implicit (via θ and the opacity function) and (likely) nonlinear variation of µ.

Solving for e + Using Exponential Differencing
We choose not to use the backward difference scheme on the material energy equation, (50), because it can result in predicting noticeable temperature differences between matter and radiation in cases where the coupling is large enough that these should be in, or close to, equilibrium, as we shall see below. We are able to avoid this problem by treating Eq. (50) as an ODE -possible because the coordinate x is there only as a parameter; no spatial derivatives appear. This will let us find the actual time dependence of Φ (and θ and the variables which depend on it) during the timestep interval ∆t, allowing a more accurate evaluation of the energy transfered during the timestep. If we define a variable with the dimension of time, then Eq. (50) implies From this we see that τ is the scale time for radiation-matter coupling (strictly speaking, τ is the relaxation time for Φ and E coupling in a gradient-free medium) and it is the proper measure by which we can quantify our earlier statements about µ being large, leading to Eq. (44) or µ being small, leading to Eq. (48) (i.e., τ small compared to the time of interest or vice versa).
Viewing t as the dependent variable and E + as a simple fixed parameter we can rearrange Eq. (53): This implies that as the interval time increases monotonically from zero toward infinity, the difference between E + and Φ relaxes from its initial value, (E + − Φ − ), towards 0. τ , which is positive and finite on physical grounds, affects only the rate at which the equilibrium is approached. Equation (54) has the formal solution, and inversion of this equation gives Φ + as a function of E + . This in turn gives the temperature θ + at the end of the timestep and, via the equation of state, e + . Thus we have in principle:

Inverting ∆t(Φ + )
The method RAGE uses to invert Eq. (55) is the subject of the next few sections. We begin with an illuminating special case.

An Ideal Case
It is very common that µ (for ionized materials) varies approximately as the inverse cube of the temperature θ while C V varies slowly by comparison. Figure 17 shows, at left, a log-log plot of the Planck absorption opacity µ of aluminum at one tenth normal density, ρ = 0.27 g/cm 3 , as calculated by the TOPS opacity code at LANL [70]; at right is a plot of the specific heat C V of the same material, calculated using a Saha-based equation of state [71]. These plots suggest, at least in this case, a) that a 1/θ 3 approximation is roughly correct and b) C V is constant by comparison.  Let us idealize this behavior, for the sake of an example, to the case where µ varies as θ −3 and C V is constant. Then τ is also constant, and the integration and inversion of Eq. (54) becomes trivial: explicitly displaying the above-mentioned relaxation toward E + . (The same favorable circumstance of constant τ occurs for gray Pomeranium, a (mythical) element with a material specific heat C V ∝ θ 3 , but with a constant opacity [65].) The fully backward differenced matter equation, (51), with the same assumptions, would replace Eq. (53) with giving the solution, If we recall that the equilibrium diffusion model, Eq. (47), gives we see that all three schemes' solution can be written as with α = 0, α = 1 1+∆t/τ , or α = e −∆t/τ , for equilibrium diffusion, fully backward differencing, or exponential differencing, respectively.
To repeat, the decay toward equilibrium in our scheme proceeds as exp (−∆t/τ ) whereas the backward scheme goes like 1/(1 + ∆t/τ ). This means that, although both agree to first order in ∆t/τ , the exponential scheme always transfers more energy between matter and radiation. Although it does not insist on equilibrium, it gets very close as ∆t/τ gets large. We expect our method to do better for intermediate timesteps where ∆t > ∼ τ . Figure 18 shows these differing solutions over a range of a few time constants; Appendix A discusses the application of this technique to the case when the opacity has a power-law behavior other than µ ∼ θ −3 , and we discuss arbitrary opacity behavior in the next section.

Φ in the General Case
As can be seen from our example opacity plot, Fig. 17 (left), the inverse cube behavior of µ is common but not a law, nor, as seen on the right, is C V always a constant. For aluminum in the temperature range 10-30 eV, both actually increase with θ, and τ turns out to vary approximately as θ −5/2 . RAGE's strategy for handling such cases consists of a more careful treatment of Eq. (55). It is invoked if we find that the constant-τ approximation fails, as given by this criterion: evaluate τ at Φ − , use it with E − in Eq. (57) to find Φ + and hence θ + (and thence C V and µ) and reevaluate τ . If the two values of τ are substantially different, then the timestep ∆t is broken into M subdivisions ∆t (m) small enough that τ is reasonably constant over each one. Equation (55) is evaluated piece by piece with τ time lagged in the integral on each subinterval, Eq. 30, equilibrium Figure 18: Three models for α so that we can do the integration and inversion for that subinterval: Expanding Φ (m−1) in terms of Φ (m−2) , etc., then gives which lets us find the new time constant τ (m) : and we repeat until ∆t is reached. We emphasize that ∆t (m) is chosen first, then Φ (m) is found from it. Defining an unscripted α, and using Eq. (65) we have, as before, giving the desired Φ + in terms of known quantities, even though α is no longer a simple exponential. This subcycling need be done only for those cells that are undergoing large material property changes in a given timestep, for instance when the opacity does not approximate the common inverse cube law or the specific heat increases suddenly as the result of beginning to burn out an ionization stage. In practice this is a rare event. The overall cost in time and load balance is small anyway compared to the cost of a diffusion solve, at least in 2-D or 3-D, since inversion of the diffusion operator is not involved in this procedure. We also reap the benefit of treating non-ideal cells accurately without slaving the timestep of the entire mesh to the small values required by a few stiff zones. For example, a non-subcycled backward difference code might require a timestep based on a relative temperature change of 5-10%, while we can allow an order of magnitude change without adverse effect (if this is the only stiff effect).
Notice that we have separated the problem of handling material properties that have any, possibly tricky, variation with temperature from the problem of dealing with the so-called ratio of specific heats. This latter problem involves only the treatment of the behavior of θ 3 , simpler in principle than dealing with the vagaries of temperature variation of opacity or specific heat, and is with us even in the case of equilibrium diffusion, Eq. (47). We treat it in section 7.7.1.
Note that α is always between 0 and 1. For the simple case of constant τ this is obvious because α is an exponential of a negative real (α = e −∆t/τ ). For the more complicated case of general τ , α is a product of M factors, each of which is bounded by 0 and 1. We can test the accuracy of the approximation that τ is a constant (over the subdivision) by comparing successive values of τ (m) , and adapting the subcycle timesteps appropriately (e.g., given the logarithms of the ratios of the two τ 's and the two θ's, we can parameterize τ as a power law in temperature (τ ∼ θ p ). Requiring δτ Using Newton-Raphson, one can solve for ∆t (m) . )

A Practical Illustration of the Iterated Φ
A valuable side effect of this approach is that it has the capability to handle cases in which the opacity has a local maximum, and this approach can actually provide the correct solution even if some zones start the timestep below the opacity peak peak and end above the peak: one merely marches through the difficult region using Eqs. (61)- (66). An example is given by the gray opacity of oxygen at ρ = 1.3 · 10 −4 g/cm 3 , as it responds to an 80 eV radiation front. Figure 19 shows that the opacity increases by almost 10 3 between 1 and 10 eV before it drops back over the next decade of temperature (SESAME's EOS 5010 has C V constant over this interval). Without subcycling, one has to hope that τ is constant for Φ to be correct, and when µ is not proportional to θ −3 , this means limiting relative changes in τ to a small value, say 10%. However, as Fig. 19 shows, for temperatures between 1 and 10 eV, µ ∼ θ +3 , so τ ∼ θ −6 , and a 10% change in τ requires no more than a 1.66% relative change in temperature. Figure 20 shows a non-equilibrium Marshak wave in relatively low-density oxygen (idealized from the gas-filled hohlraum that raised this issue). The oxygen sits in a (relatively) high temperature radiation bath until enough ionization occurs to allow it to rapidly heat up. The solid black (θ(x)) and solid red (T R (x)) curves shows a converged calculation, run with 2000 timesteps to a final time of 20 picoseconds. The various dashed curves show profiles of T R (x) and θ(x) at t = 20 ps, when RAGE was altered to not do any opacity subcycling but to run on a timestep control that kept the relative temperature change below 20, 10, 5 and 2% (requiring 122, 254, 554 and 1493 timesteps, respectively). Only the last calculation has converged the matter temperature.
When RAGE has adapted an interface in order to get the correct energy deposition into a wall, say between air and gold in a hohlraum, it has to adapt in the other direction(s) as well.     Fig. 20, RAGE has very small zones in the direction that the Marshak wave is propagating, even though it does not need them in that direction, and would prefer to spend as little time as possible calculating the flow in that direction. Contextually, "little time" means flowing radiation on a radiation Courant timestep, ∆t ≈ ∆x/c, taking 1 or 2 timesteps to cross each zone. Since the level 1 mesh size is 0.01 cm in this example, this means ∆t ≈ 3 · 10 −13 s, requiring 67 timesteps to get to the displayed time of 20 ps. If one used traditional timestep controls, allowing 2% relative temperature changes from 1 eV (when opacity begins to change) to 40 eV, one would require 186 timesteps to cross each zone. That becomes prohibitive given the size of 2D and 3D AMR matrices and the time required to invert them. Figure 21 compares the results with and without the opacity subcycling when the timestep is forced to be the radiative Courant timestep, and both are compared to the 2000 timestep converged result. The width of the leading edge of the radiation temperature reflects the O(∆t) accuracy of the (flux-limited) diffusion solution; it can only be improved with smaller timesteps, and even c∆t/∆x = 1/2 will be much better. If RAGE had not been able to subcycle and thereby dispense with traditional timestep controls, it would not have been able to run this type of problem.

dΦ + /dE +
Recalling that the equation we need to solve involves ∂(ρe)/∂t, and using the chain rule, we can now write As we will see presently, our solution scheme will require the time-advanced dΦ + /dE + . For 1/θ 3 opacities and constant specific heat, τ is constant, α = e −∆t/τ , and dΦ + /dE + results trivially from Eq. (60): In the general case we still have Eq. (68), but since α depends upon E + , we must compute the derivative in a more roundabout way. By taking the derivative with respect to E + of Eq. (55) and doing a little manipulation, we have Again we break the integral into smaller pieces. The added cost is small because we can recycle the partial results from Eqs. (61)-(66), building and calculate this at the same time we calculate α. By working out the telescoping sums of the weights, we see that m W (m) = 1 − α, and rescaling to create normalized weights w m that sum to unity, we can define an average τ , enabling us to rewrite Eq. (71) as with α given by Eq. (67). This means that in the general case we have the result for a constant τ , Eq. (69), multiplied by the ratio of the average of τ over the timestep to its value at the end of the timestep. In all but the most perverse cases this ratio will be of order unity (RAGE does not depend on this supposition, using Eq. (71) directly).

The
Equations for E + and e + .
Equation (68) tells us that Φ + can always be written as αΦ − + (1 − α)E + , regardless of the behavior of µ and C V . What about e + , which is required to solve for the radiation energy at the new time, as seen in Eq. (49)?
We use exactly the same logic as in the case of equilibrium diffusion, described above in section (7.2.3). The equation of state E(θ) gives where Φ + has been used to find θ + , as per Eq. (34). Hence we have e + as a function of Φ + .

The Equations to be Solved.
Finally, we specify LD as evaluated at the forward time. Our equations then become the nonlinear coupled set, for which we invert a matrix equation for E + , and then solve for the remainder, where Eq. (73) represents the diffusion equation (49) and Eqs. (74) . . . (76) result from the exponential treatment of the radiation-matter heat exchange, Eq. (50). This set is to be solved at every timestep for every space point, thus advancing E and e from time t − to time t + . Their subordinate variables can then also be advanced.

Spatial Differencing
All cells in RAGE have rectilinear cross sections with unit aspect ratio in every direction, and all faces are orthogonal to the appropriate coordinate axes. Therefore only one component of a gradient will contribute to a divergence on any face in any RAGE geometry, and spatial differencing becomes quite simple. For example, if an interface normaln points in thex direction, then for any scalar s: − → ∇s ·n → ∂s ∂x , and this O(∆x) expression for the face gradient suffices to build O(∆x 2 ) zone-centered divergences.

Touching Cells
There are two possible configurations of the geometry of touching cells: either both cells are the same size or else one cell is half the size of the other (a situation we refer to as T-cells, given the form of the interface ([ ], as in Fig. 22). In 1D, adaption does not affect the size of the interface, but in 2D or 3D the larger cell communicates through two or four interfaces with the two or four smaller cells.

Differencing for Same Size Interface Cells
Same size interface cells in RAGE have the same size in the dimension(s) perpendicular to the line joining cell centers. In this sense all 1D cells qualify, as do all other geometries except across interfaces where cell adaption has occurred. RAGE's difference stencil is simple. For every interface not subject to a boundary condition, we require that the flux calculated between the cell centroid and the interface on one side be equal to the flux between the interface and cell centroid on the other side, and that they be equal to a net flux between the two centroids: ∆s Zf and ∆s f Z are the distances between the appropriate cell centroid (Z = L or R) and the interface f (half the zone size along Cartesian coordinates), and the τ are optical depths (D = λc/3 ⇒ ∆τ = ∆s/λ). Solving these equations gives the face quantities: implying that optical depths add: ∆τ LR = ∆τ L + ∆τ R .

Temperatures on a Face
Experience has shown that if the two diffusion coefficients (D L , D R ) are evaluated at the zonecenter temperatures (θ L , θ R ), it is next to impossible to use Eq. (77) to calculate properly the self-similar solution of a Marshak [72] wave. The cold material has a very small D Z (very large ∆τ Zf ), hence ∆τ RL is also large and thus very little radiation flows -according to the difference equations. The solution to this problem has been to find a common face temperature, use it to define D L (θ f ace ) and D R (θ f ace ), and then compute the ∆τ L , ∆τ R and ∆τ LR from those D's. Even though the temperature(s) are now the same, these D's will differ if the densities or constituents of the two cells differ. Flux equality requires that E vary linearly in optical depth, Eq. (78), and we note that if radiation and matter were in thermal equilibrium (E = Φ) then θ 4 would vary similarly. We adopt this ansatz in all cases, using a two step approximation: First we calculate D's using the current cell-center temperatures. This lets us define optical depths and so find θ 4 f ace by interpolation: We then use this value of θ f ace for the "real" lookup of opacities, D's, and ∆τ 's. Some alternate ways of estimating θ f ace are discussed in Appendix B; we have not found them to make enough of a difference to offer users a choice in the matter.

Differencing for Unequal Size Interface Cells
In the case where adaption has occurred, as depicted in Fig. 22, we "splat" (extend central values across the bigger cells) and continue to use Eq. (77), for F 03 and F 01 . This leads to a formal loss of spatial accuracy, but is mitigated by the fact that in our AMR implementation the gradients are relatively small at any place in the mesh where such scale size changes occur. Splatting has the advantage that the resulting matrix couples only cells which share an interface. Edwards [73] introduced a fix for this loss of accuracy, in the context of petroleum extraction, that he called the "Flux Continuous Approximation" (FCA). Lipkinov et al. and Olson et al. [74,75], re-derived it in the context of support-operators (such that div and grad are consistent) for diffusion. The argument is summarized briefly in Appendix C.
FCA was coded for RAGE in 1996 for 2D geometries. Unfortunately in 3D FCA requires a coupling between cells that touch only at an edge. This was not allowed by RAGE's data structure at that time. In order that the 3D code reproduce 2D results where expected, RAGE chose not to use this technique in 2D either.
This is not such a blemish on the code as might be expected from the resulting formal loss of accuracy. For one thing, RAGE's AMR logic places such scale changes only in regions with small gradients (gradients are the trigger for adaption, after all). In the case of isotropic coefficients and a 2:1 split in cell sizes -both always the case in RAGE -simple splatting, observed Edwards [73], "can produce extremely good results".
In any event, recent changes to the data structure should now allow such couplings, and we expect the FCA treatment to appear in the next major revision of RAGE.

Boundary Conditions
RAGE users have the option of four different types of boundary conditions for radiation: Neuman, Dirichlet, and two forms of Robin.
The default condition is to assume a Neuman (mirror) condition, ∇E = 0, at all edges of the mesh. (This is reasonable for heat transfer mediated by ions and electrons that one does not expect to leave the mesh, but less so for radiation, since one generally does expect photons to leave the problem if they get to the boundary.) All other boundary conditions are specified by nonzero values of temperatures on the appropriate face. The user can choose to impose either a Dirichlet or a Robin boundary condition on all boundaries with non-zero temperatures (e.g., in 3D, one could have all six boundaries at six different temperatures, as long as they all use the same condition). The Dirichlet condition holds the temperature to a fixed value, T f at face f ; of the two types of Robin condition, the Marshak/Milne holds it at that temperature but 2/3 of an optical depth beyond the face, and the Spillman variant [76] of Marshak/Milne holds it to that temperature at a variable optical depth beyond the face. The Marshak/Milne boundary condition is designed to reproduce the effect for a single energy equation that the transport equation would accomplish with two equations; it requires optically thin boundary zones to succeed. Spillman's variant is conceptually similar to, but considerably simpler to implement than, the Levermore-Pomraning boundary condition [64] that enabled their flux-limited diffusion to recover transport solutions in thin and thick limits. We compare and contrast the Robin conditions in Appendix D.

Solution Techniques for Non-Linear Systems
As mentioned above, Eqs. (73). . . (76) are very similar to the result of a backward difference scheme, so it is of interest to examine techniques for solving that scheme. Shestakov et al. [58] discuss the use of pseudo-transient continuation to stabilize a Picard-Newton series of linearizations and iterations, and Knoll et al. [69] discuss wrapping a Newton-Krylov iteration around an inner iteration. We might well profit from the use of such methods, but at present RAGE uses a different (nameless) method, also involving linearization and iteration.

RAGE's Technique
At each timestep we execute an iteration process, with iterates {E We identify this as a separate step in the iteration process because more elaborate methods for this step are currently under study, having the flavor of the prediction step in a predictorcorrector technique.
The second step consists of solving a linearization of the radiation diffusion equation, in which the coefficients are evaluated in terms of the estimated values E where the coefficient (LD) The RHS of Eq. (82)   Early versions of RAGE used a conjugate gradient linear solver with a point Jacobi (diagonal scaling) preconditioner. We have now converted to an algebraic multigrid preconditioner (LAMG) developed at LANL by Joubert [77]; this preconditioner has typically given, for moderate to large 2D and 3D meshes, a factor of five or better in speed-up of the linear solver, relative to the point Jacobi preconditioner.
Owing to the linearization in the second step, the values {E  }. In order to decide when to terminate the iterative process, we (currently) look at the total energy in the cell, exiting when changes in that quantity becomes sufficiently small. We are also investigating the more usual method of looking at the nonlinear residual. Normally, one or two iterations suffice.
Our experience has been that, for small to moderate timesteps, the dominant nonlinearities come from the radiative heating and cooling associated with the material energy equation. For large timesteps, the nonlinearity associated with the diffusion coefficient and flux limiter can become more important. In fact, for large enough timesteps ∆t, we may encounter diffusion fronts that move faster than ∆x/∆t, where ∆x is the size of a computational cell; in this case it is necessary to execute several iterations in order for the front to traverse the required number of cells.
The great bulk of the cpu time (at least for 2D and 3D meshes) is normally spent in executing the linear solve in the second step of the iteration (and possibly in the first step in future versions of the code). No other computations require inter-cell communication; they execute relatively quickly, particularly in a parallel environment.
The current version of the code makes the simplifying approximation that the specific heat C V is held fixed during the radiation diffusion computations. This simplification eliminates a number of equation of state calls, replacing Eq. (85) with the simple linear relation e + = e − + C V (θ + − θ − ). Future versions of the code will include as an option the elimination of this simplification. Figure 23 shows a comparison of the results of a test problem taken from Shestakov et al. [58].The RAGE results for Φ (aθ 4 ) at 20 ns are shown as histograms and the original calculation using pseudo-transient continuation (Ψtc) as a curve through cell centers. As one can see, the temperatures (or Φ's) calculated by the two codes (with the same zoning and timesteps) are in very close agreement. For this problem RAGE was modified to use the same hydrogenic Saha equation of state used by Shestakov et al. rather than SESAME tables, in part because the SESAME EOS includes the effect of dissociation of molecular H 2 , not just ionization of atomic hydrogen (as specified by the test problem). RAGE was also modified to use the n = 1 Larsen flux limiter [57] (RAGE's Levermore limiter agreed with the n = 2 Larsen limiter, and both disagreed with the n = 1 limiter by a little more than the thickness of the lines in Fig. 23). Figure 24 shows the temporal behavior of a related problem. The configuration is changed to an infinite medium with fixed C V but the same µ a ∝ T −7/2 opacity, so we can focus on the coupling of matter and radiation only (i.e. without diffusion complications). The plot shows the electron temperature as a function of time, up to 20 ns. The problem has been run three times, varying the integration timestep ∆t. The smoothest curve was generated by using 200 timesteps, each of size 10 −10 s. Another curve used only 20 timesteps of 10 −9 s, and finally the solid curve used only two(!) timesteps of 10 −8 s. Clearly one loses information as to what happens in the intervals between such monster timesteps, but we call attention to the fact that RAGE's exponential difference scheme gets the right answer at a given time regardless of the timestep -as long as opacity subcycling is allowed. Figure 25 shows temperature profiles for a self-similar problem due to Reinicke and Meyerter-Vehn [78] and popularized by Shestakov [79, his figure 4, "E0=235"]. A thermal wave and a hydrodynamic shock are created whose radii maintain a fixed ratio, generated by a point explosion in cold material with ρ ∼ r −19/9 . The solid black curve of the figure shows the resultant  Figure 23: Comparison of aT 4 spatial profiles at 20 ns from RAGE (thin black histogram) and Ψtc [58] (smooth blue curve).

Comparisons to Other Codes and Standard Problems
temperature profile using the heat conduction package in RAGE (with the appropriate analytic conductivity, κ = 31.62 · ρ −2 T 6.5 ); the red dot-dashed curve shows the radiation temperature T rad and the blue dashed curve the material temperature T mat calculated by RAGE's radiation diffusion package, all at a time of 5.15 · 10 −10 s. T mat and T rad do not overlay the conduction solver's temperature profile because the self-similar problem was defined with a constant C V [78]; unlike the conduction package, the radiation diffusion package cannot help including an extra 4aθ 3 in the total specific heat, a significant component at these temperatures that destroys the self-similarity of the solution. The more interesting point is that RAGE's exponential differencing has kept the two temperatures in equilibrium (at an early time of 1 · 10 −10 s, there is a slight disequilibrium near the origin); we expect standard backward differencing methods would have trouble keeping the temperatures clamped together (recall Fig. 18). We did, however, assume the freedom to define a Planck-weighted material absorption opacity κ abs = 5.433 · 10 12 ρT −3.5 , 31.32 times larger than the presumed Rosseland total opacity κ tot = 1.735 · 10 11 ρT −3.5 , making relaxation times τ ∼ 10 −12 −10 −11 s. Without this, the two temperatures were out of equilibrium everywhere, and the thermal front was another half millimeter retarded (τ ∼ 10 −11 − 10 −10 s). Finally, we present an example of the ability of RAGE to deliver reasonable results even under unreasonable conditions. In particular, we display solutions of a Marshak-like problem, with a left boundary held at 1000 eV. Cell 1 is initialized to that same temperature, but the rest of the (1D) grid is started at room temperature. We use a dummy material with a constant specific heat, C V = 10 12 erg/g/eV, constant density, ρ = 1 g/cm 3 , and variable opacity, Each cell is 0.02 cm wide. We run to 10 −11 s, at which time a simple constant-flux-approximation argument suggests that a wave would have burned through about 6 1 2 cells. RAGE normally  chooses its timesteps fairly carefully, beginning with a very small one and inching up as seems appropriate, but for this series we forced it to use a constant timestep. We ran the problem three times, taking one, two, and fifty timesteps ( c∆t ∆x ≈ 15, 7, 1/3) to get to the final time; the results are shown in Fig. 26 (radiation temperature vs. space) and Fig. 27 (electron temperature vs. space).
We do not claim that the results obtained with the ridiculously large timestep are accurate (after all, the backward difference scheme used for the diffusion equation is only first-order accurate), rather that they are not unreasonable; in particular, since we have no theorem that proves our iteration method will converge, these results are encouraging in that they do not diverge. Moreover, they show that when multiple (expensive) physics packages are turned on (e.g. the gravitational Poisson solver), we can iterate the matrix-solves in the radiation package without requiring the rest of the code to run on the smaller effective timestep.

Radiation Summary
We have described just a subset of RAGE's radiation capability: gray diffusion. RAGE also contains a multifrequency diffusion package, in which the spectrum of the radiation energy field is calculated by the multifrequency-gray method [61], and of which this algorithm is a subset, and a gray P 1 package which contains a number of novel features. The gray diffusion algorithm described here is of note primarily because of the exponential differencing of the material energy equation and of the novelties in the iteration scheme. Exponential differencing is particularly useful in a class of problems in which radiation floods a region of space and serves to heat a contained body. This situation occurs frequently in inertial confinement systems, where the radiation source is produced by e-beam or laser deposition in a hohlraum wall and the goal is to heat and compress a capsule containing the nuclear fuel.

Verification and Validation for SAGE and RAGE
During the course of the ASCI (now ASC) Project, substantial time has been devoted at LANL, either directly or indirectly, to the verification and validation (V&V) of the Project codes. Direct V&V results from a Project team member running specific simulations to quantitatively compare to analytic results (verification) or to results from experimental data (validation). Indirect validation results from the use of the Project codes by end-users for their own purposes, either in designing a new AGEX (Above Ground EXperiment) experiment, in the analysis of a previous experiment, or merely as a tool for analysis of a physical problem.

Verification
There exists a large number of verification test problems in our community, including the Sedov exploding blast wave [80, pp. 393-396], the Guderley converging shock, [81], [82, p. 793], anthologies of hydrodynamic [83] and radiation-hydrodynamic problems [84], both LTE and non-LTE [72,85,66], Marshak radiation wave problems, and many more (non-LTE in this context means that T r = T e ). A detailed report on the suite of verification problems is being prepared (see e.g. [86,87]), as well as work on automating such ongoing testing [88].
One example of such a verification problem is the Sedov problem, with unit density and zero initial velocity but 4.936 · 10 15 erg in a delta function at the origin (i.e. in one zone at the origin, so that r shock = 1 cm at t=10 ns). The results from the (2d) cylindrical Sedov problem are shown in Figure 28, comparing the SAGE adaptive mesh hydrodynamics to the analytic result,   and some error norms are provided in table 1, which incidentally shows that the same accuracy can be achieved in less time and fewer zones using adaption. The ratio of L 1 (ρ, u, p) errors for the fixed meshes is consistent with those errors being proportional to ∆x 0.6−0.75 , typical for most shock problems we have looked at (shockless problems without discontinuities do behave like ∆x 2 ).  Table 1: Comparisons of L 1 error norms for density, pressure and velocity on 4 fixed size meshes as well as a mesh that adapts ("A1-4") to the finest level. Also shown are measures of the time required to calculate the answer (on a Mac G5).
In general we are very pleased with the success of the SAGE and RAGE algorithms in comparing to analytic and semi-analytic results (see e.g., [89] for an example of comparisons to semi-analytic problems and Section 4.5 for four comparisons to analytic solutions). Many of these verification problems are included in a regression test suite (run on the LANL machines) as part of the software quality assurance process.

Validation
Our first validation effort, [90], "The Simulation of Shock Generated Instabilities" compared detailed code simulations of shock tube experiments to the experimental results in a quantitative manner; the results from this ("gas curtain") work are summarized (in the "vu-graph" norm) Figure 29: Comparison of experimental and calculated density profiles for a gas curtain shock-tube experiment. t = 0 profiles provided initial conditions. in figure 29. In order to agree with the data at t ∼ 450 ms, the code was initialized with the density profile measured at t = 0 ms (after background subtraction and some smoothing). Further Richtmyer-Meshkov validation was done using shocks generated by lasers [91] and pulsedpower [92]. RAGE has also been used to model more general laser-induced shock effects [93,94,95] and more recently, 3-D oceanic disturbances generated by asteroid impacts and tsunamis [96,97,98,99].
Code validation also formed the basis for a Ph.D. thesis for a SUNY Stony Brook student [100,101]. In this work, code simulations were done for a similar shock tube environment but in this case the Mach 1.2 shock interacted with a dense cylinder of gas, instead of a gas curtain. Both the code and the experimental techniques had changed by the time of this comparison (the code was rewritten from Fortran77 with a Lagrange plus Remap hydro to Fortran90 with the Direct Eulerian hydro described in this paper; the original experiments [90] used Rayleigh scattering as a diagnostic while the new ones used "fog" tracer particles, which were thought to be equivalent to Rayleigh images and were essential in order to perform particle velocity measurements) and the agreement between the two was not as good as it was in the earlier work. Having determined that changing the hydro algorithm didn't change the discrepancy with experiment, changes to the initial conditions were tried -it was hypothesized that the SF 6 diffused faster than the fog particles before the first image was taken, leading to incorrect initial density profiles -and much better agreement was obtained. The top part of figure 31 shows how much the initial RAGE density profile had to be diffused in order to match the (late-time) data (shown in figure 32), compared to the original profile based on the "fog" tracer particles (left half, figure 30). This prompted the experimenters to re-image their initial profiles using a Rayleigh scattering technique (right half, figure 30). The lower part of figure 31 shows how RAGE's diffused density profile (which gave good late-time agreement) compared to a (later) Planar Laser Rayleigh Scattering (PLRS) measurement. 8 In this case, one could say that the code validated the experimental technique. On the other hand, the advantage of the fog tracer particles is that they enable the experimenters to measure the velocity field in the gas cylinder (from two closely spaced images). When the calculated and simulated images of the density agreed (with the diffused density profiles), the velocity data were also found to agree, to with 10-15 percent, and one can say that the experiment validated the code in this case.
Another example of validation is the work being done as a part of a multi-lab, multinational, collaboration on the simulation of supersonic jets and shock interactions [102], where detailed quantitative inter-comparisons have been done among several computational tools (including RAGE) and the actual experimental data obtained from a series of AGEX ("Above Ground EXperiments") experiments. All the code techniques, including CAMR RAGE and ALE codes from both Livermore and AWE (the UK Atomic Weapons Establishment), show good agreement with the data. In this and other cases [103], results of code simulations are used to design as well as predict future experiments.
These examples touch on the breadth of validation efforts completed and underway in the ongoing effort to understand the range of capability for the CAMR Eulerian-based hydrodynamics and radiation energy flow as coded in the massively-parallel, modular version of the RAGE code. More detailed documentation efforts are underway. Two of the areas of active research in the Crestone Project currently are the treatment of interfaces between materials and the improvement of the material strength package.

Summary
We have presented a description of the SAIC/LANL Radiation Adaptive Grid Eulerian code, RAGE, describing some of its data structures, its parallelization strategy and performance, its hydrodynamic algorithm(s), and its (gray) radiation diffusion algorithm. We have also shown that all packages and the code as a whole are subject to a considerable amount of verification and validation efforts.
The Godunov-based hydrodynamics uses a second-order alternating-direction explicit (ADX) method, and this algorithm determines the maximum timestep the code can take. The iterative nature of the anti-cavitation logic handles stiffnesses in the equation of state in an implicit manner, meaning that we do not have to control the timestep based on, for example, relative changes in density or mass fraction.
The gray radiation package uses a temporally first-order and spatially second-order backward Euler method to solve the diffusion equation. The material energy equation is integrated exactly for opacities proportional to T −3 , and to any specified accuracy for all other behaviors by subcycling the numerical integration. This solution is then linearized to build the coupled total energy diffusion equation, and a subsequent Newton-Raphson step rebalances the material and radiation energy to correct for any linearization errors. Finally, an outer iteration is performed until the total radiation plus material energy in each zone is relatively unchanged between two iterations.
Because the radiation package has taken great pains to treat (1) stiff material coupling terms (µ abs or C V ) by subcycling, and (2) linearization errors by Newton-Raphson fixups, it can be considered to be maximally implicit. This means that its accuracy is completely determined by the relative changes in radiation energy, independent of the changes in material temperature.
RAGE by default attempts to take relatively large timesteps, based on controlling just the relative change in total energy over the entire timestep, δ(ρE M + E)/(ρE M + E), to some tolerance, e.g., sie pct = 0.2. It is also possible to control the material temperature change on shows how reproducible the experimental conditions were. Figure 30: Experimental images of a gas cylinder, using "fog" tracer particles on the left, and Rayleigh-scattering on the right. Figure 31: Top: original (measured) and "fitted" (diffused) initial density profiles for a gas cylinder experiment; the original was an (angularly) smoothed measured fog distribution. Bottom: comparison of that "fitted" profile to the measured Rayleigh-scattering profile (n.b., scales have changed). Figure 32: Experimental images of a gas cylinder, using "fog" tracer particles on the left; calculated images using sharp (middle) and diffuse (right) profiles from Fig. 31 top.
each step of the operator splitting to a tolerance (tev pct), but this is a legacy of Lagrangian technology. We intend to update the splitting logic to instead control relative changes of material and radiation energy on each step of the operator split, and may couple this with an overall control on the non-linear residual.
The implementation of these implicit treatments allow our AMR code to run on timesteps that are considerably larger than would be allowed by the traditional small-percentage-changes in material and radiation temperature on each step of an operator split, and this means that our AMR code can run more highly-resolved problems than might otherwise have been possible, within time and resource constraints.

A A Less Ideal Opacity
The case of constant τ is but one example of a whole class of situations in which exponential differencing (55) can be treated analytically: that for which τ has a power-law dependence on θ, and thus Φ. Given that Eq. (55) becomes where any φ represents the corresponding Φ scaled to E + : 2 F 1 (λ + 1, 1; λ + 2; φ) is of course the confluent hypergeometric function of the second kind [104]. This is a very nice form: since Φ − , λ, Φ * , and τ * are all known, the left side of Eq. (87) is just a constant multiple of ∆t (also known). The right side, on the other hand, is a function of both E and Φ + . This then provides all we need to find Φ + as a function of E; if need be we can resort to a numerical root-finding routine to extract the information.
This chore can be rather intimidating in general, but is much easier when λ is a multiple of 1/4, because then the hypergeometric function can be expressed in terms of standard simple functions.
For example, consider the case λ = − 1 2 , as occurs with a constant specific heat and 1/θ opacity. Then since 2 F 1 (l + 1, 1; l + 2; φ)| l=−1/2 = Eq. (87) becomes: Using this we can explicitly display φ + in terms of a simple parameter, X, where X is Note that we used Eq. (86), the scaling law for τ , to evaluate τ E . Earlier we discussed three different simple models for the parameter α, which governs the energy transferred between matter and radiation during a timestep, Eq. (60). It is of interest to examine α for this new case: Figure 33 shows the results for three different values of φ − , a parameter heretofore irrelevant. The bottom curve, for φ − = 0, is the same as the curve labeled "exponential" in Fig. 18. Again this more sophisticated treatment yields seemingly more reasonable results. Instead of assuming that θ 4 varies linearly with optical depth across an interface, another choice would be to make the minimalist assumption that the zone-centered temperatures vary linearly between zones and interpolate a material temperature on the face, unless the problem is in radiative equilibrium, in which case we might expect E and θ 4 to vary linearly in space, leading to If one knows that if the problem is in a completely ionized regime, so that the opacity is dominated by free-free (Bremsstrahlung) physics, the mono-frequency opacity goes like T −1/2 e /ν 3 . When this opacity is averaged over frequency using as a weight the Rosseland at the radiation temperature (since we are talking about gray diffusion, the only questions are whether the flux is a Planckian or Rosseland, and whether it is at T e or T r ; we make the latter choices), we end up with a "two-temperature" opacity with mean free path λ ∼ T 1/2 e T 3 r (a two-temperature opacity is the result of using a non-equilibrium Planckian or Rosseland radiation field to calculate a gray average: e.g., κ(θ e , T r ) ∼ ∞ 0 κ ν (θ e )B ν (T r )dν). Ignoring the weak dependence on material (electron) temperature, we can write this flux as F ∼ T 3 r ∇E ∼ E 3/4 ∇E. If we ignore quarter powers of things, we can say that the flux behaves approximately as F ∝ ∇E 2 , whose finite difference approximation would be (For electron conduction, the Spitzer conductivity varies as θ 5/2 e , leading to F ∼ θ 5/2 e ∇θ e ≈ θ 2 e ∇θ e . Since ∇θ 3 e /3 = 1 3 (θ 2 L + θ L θ R + θ 2 R )∇θ e , we have found that the intermediate temperature, θ f ace = (θ 2 L + θ L θ R + θ 2 R )/3, in Spitzer's formula gives good results on this variant of a Marshak wave.) From this point of view, one would argue that the face temperature ought to be related to the average of the two energy densities without any linear interpolation (in space or optical depth): It has been our experience that whether one uses the assumption of linearity in optical depth space for θ 4 (Eq. 88), or either of the linear in space methods (Eqs. 89,90), all Marshak-wave types of problems can be well modeled. Although we have not looked at the last averaging method, Eq. (91), in any detail, it ought to give results similar to the linear-in-space average of θ 4 in most cases, and might be the best approach to use with two-temperature opacities in general.

C Improved Difference Scheme for T-Cells
When adjacent zones are equal-sized, as for the case of vertical flow between zones "1" and "3" in Fig. 34, the assumption of piecewise continuity of E along with the requirement of continuity of flux from zone 1 to the (13) face and from that face to zone 3 gave the two equations that defined E f ace as the optical-depth weighted average of the zonal E's and τ f ace as the sum of the zonal optical depths. (In principle, E above 13 = E below 13 , but piecewise continuity means that dV ∇E → 0 around a tiny volume straddling the (13) interface as dV → 0, and that forces E above 13 = E below 13 = E f ace .) For T-cells (reverting to the horizontal direction of flow in Fig. 34), piecewise continuity of E means that as dV ∇E = dA f E f → 0, the three non-coincident values of E at the respective faces of the zones are related by: A 0 · E A f ace = A 3 · E C f ace + A 1 · E B f ace , where A 0 is the area of the large zone's face at the T-cell interface, and A 1 and A 3 are the areas of the interface as seen by zones 1 and 3 respectively.
In 2-D Cartesian geometry, the A i = A 0 /2; in 3-D, they are a A 0 /4. In 2-D cylindrical geometry, for face-normals in ther direction, the A i are also A 0 /2, but for normals in theẑ direction, they vary because dA = rdr, and r varies (the smaller radius face varies from 1 4 A 0 at the origin to 1 2 A 0 approaching infinity). Let us define a i ≡ A i /A 0 . By requiring continuity of the fluxes at T-cells and using the relation E A = a 1 E B + a 3 E C , we have the following two constraints: and the mean flux is similarly given by The Spillman variant [76] is an attempt to cover more bases. Assume that the user specifies a (dimensional) vacuum boundary length, λ V (default = 1 cm), and recall that λ Z was used to calculate the optical depth: ∆τ Z = ∆s Zf /λ Z . Defining the ratio, Spillman sets the contribution to so that when λ Z is large compared to both λ V and ∆x, ξ → 3 2 and the denominator goes to 1 3 + ∆τ Z . When λ Z is small compared to λ V but still large compared to ∆x, ξ → 2 and the denominator goes to 2 3 + ∆τ Z , the standard Milne result. Even if the user makes a poor choice of λ V , in optically thin zones, Spillman will always give Milne-ish behavior.
However, if the user intended to have inflow via a Milne condition, but made the zone abutting the boundary so large that its optical depth was enormous, Milne's 2/3 + ∆τ is indistinguishable from Dirichlet's τ , and the flux → 0. In this case the numerics allows essentially no flow across the boundary, contrary to what would have been allowed if the zone had been thinner, and what physics presumably always requires.
In such a case of large optical depth (∆s Zf λ Z ), Spillman's exponential term vanishes and the denominator varies between unity and 4 3 (depending on the exact value of λ V relative to λ Z ). The numerical fore-factor (∼ 3 4 · · · 1) will now be within a factor of two of the optically thin Milne result ( 3 2 ), instead of vanishing. These contributions from boundary faces do not contribute to any off-diagonal elements of the matrix to be solved, but the factor multiplying (E − aT 4 ) does add (actually, subtract) from the diagonal element of the matrix for the zone adjacent to the face. This coefficient also contributes to calculating the source term of the matrix equation (the "vacuum" "zone" shows up as an explicit source for the adjacent real zones ). In the case of thermal conduction, one calculates the divergence twice: once to form the source for the matrix solver, and a second time afterwards to calculate a conservative energy transfer between all faces based on κ∇θ fluxes; for radiation, E is a conserved quantity already, so only the first, matrix-source divergence is necessary.
Thus we can see that the optically thick limit of Spillman's variant allows the vacuum to radiate like a black body (more correctly, it radiates like the difference of two black bodies). F · dS ≈ ± 3 4 A f ace c(E Z −aT 4 b ), and compares well to the optically thin Milne result, ± 3 2 A f ace c(E Z − aT 4 b ), and disagrees with the optically thick Marshak/Milne (or Dirichlet) result, ±A f ace c(E Z − aT 4 b )/τ Z → 0.