next up previous
Next: Vectorization or Parallelization Implementation Up: Performance Report Previous: Performance Report

Algorithmic Description

Our current calculations assume an inviscid and compressible gas, and an ideal equation of state with constant adiabatic index. We use a Godunov-type solver which is a relativistic generalization of the method due to Harten, Lax, & Van Leer (1983), and Einfeldt (1988), in which the full solution to the Riemann problem is approximated by two waves separated by a piecewise constant state. We evolve mass density R, the three components of the momentum density tex2html_wrap_inline151 , tex2html_wrap_inline153 and tex2html_wrap_inline155 , and the total energy density E relative to the laboratory frame. Earlier, 2D axisymmetric calculations computed the axial and radial components of momentum density, tex2html_wrap_inline159 and tex2html_wrap_inline155 , in a cylindrical coordinate system.

Defining the vector


and the three flux vectors




the conservative form of the relativistic Euler equation is


The pressure is given by the ideal gas equation of state tex2html_wrap_inline163 The Godunov-type solvers are well known for their capability as robust, conservative flow solvers with excellent shock capturing features. In this family of solvers one reduces the problem of updating the components of the vector U, averaged over a cell, to the computation of fluxes at the cell interfaces. In one spatial dimension the part of the update due to advection of the vector U may be written as


In the scheme originally devised by Godunov (1959), a fundamental emphasis is placed on the strategy of decomposing the problem into many local Riemann problems, one for each pair of values of tex2html_wrap_inline169 and tex2html_wrap_inline171 to yield values which allow the computation of the local interface fluxes tex2html_wrap_inline173 . In general, an initial discontinuity at tex2html_wrap_inline175 due to tex2html_wrap_inline169 and tex2html_wrap_inline171 will evolve into four piecewise constant states separated by three waves. The left-most and right-most waves may be either shocks or rarefaction waves, while the middle wave is always a contact discontinuity. The determination of these four piecewise constant states can, in general, be achieved only by iteratively solving nonlinear equations. Thus the computation of the fluxes necessitates a step which can be computationally expensive. For this reason much attention has been given to approximate, but sufficiently accurate, techniques. One notable method is that due to Harten, Lax, & Van Leer (1983; HLL), in which the middle wave, and the two constant states that it separates, are replaced by a single piecewise constant state. One benefit of this approximation, which smears the contact discontinuity somewhat, is to eliminate the iterative step, thus significantly improving efficiency. However, the HLL method requires accurate estimates of the wave speeds for the left- and right-moving waves. Einfeldt (1988) analyzed the HLL method and found good estimates for the wave speeds. The resulting method combining the original HLL method with Einfeldt's improvements (the HLLE method), has been taken as a starting point for our simulations. In our implementation we use wave speed estimates based on a simple application of the relativistic addition of velocities formula for the individual components of the velocities, and the relativistic sound speed tex2html_wrap_inline181 , assuming that the waves can be decomposed into components moving perpendicular to the two coordinate directions.

In order to compute the pressure p and sound speed tex2html_wrap_inline181 we need the rest frame mass density n and energy density e. However, these quantities are nonlinearly coupled to the components of the velocity as well as to the laboratory frame variables via the Lorentz transformation:






where tex2html_wrap_inline191 is the Lorentz factor and tex2html_wrap_inline193 . When the adiabatic index is constant it is possible to reduce the computation of n, e, tex2html_wrap_inline199 , tex2html_wrap_inline201 and tex2html_wrap_inline203 to the solution of the following quartic equation:


where tex2html_wrap_inline205 . This quartic is solved at each cell several times during the update of a given mesh using Newton-Raphson iteration.

Our scheme is generally of second order accuracy, which is achieved by taking the state variables as piecewise linear in each cell, and computing fluxes at the half-time step. However, in estimating the laboratory frame values on each cell boundary, it is possible that through discretization, the lab frame quantities are unphysical - they correspond to rest frame values v;SPMgt;1 or p;SPMlt;0. At each point where a transformation is needed, we check that certain conditions on M/E and R/E are satisfied, and if not, recompute cell interface values in the piecewise constant approximation. We find that such `fall back to first order' rarely occurs.

The relativistic HLLE (RHLLE) method discussed above constitutes the basic flow integration scheme on a single mesh. Our past and planned work also utilizes adaptive mesh refinement (AMR) in order to gain spatial and temporal resolution.

The AMR algorithm used is a general purpose mesh refinement scheme which is an outgrowth of original work by Berger (1982) and Berger and Colella (1989). The AMR method uses a hierarchical collection of grids consisting of embedded meshes to discretize the flow domain. We have used a scheme which subdivides the domain into logically rectangular meshes with uniform spacing in the three coordinate directions, and a fixed refinement ratio of tex2html_wrap_inline147 . The AMR algorithm orchestrates i) the flagging of cells which need further refinement, assembling collections of such cells into meshes; ii) the construction of boundary zones so that a given mesh is a self-contained entity consisting of the interior cells and the needed boundary information; iii) mechanisms for sweeping over all levels of refinement and over each mesh in a given level to update the physical variables on each such mesh; and iv) the transfer of data between various meshes in the hierarchy, with the eventual completed update of all variables on all meshes to the same final time level. The adaption process is dynamic so that the AMR algorithm places further resolution where and when it is needed, as well as removing resolution when it is no longer required. Adaption occurs in time, as well as in space: the time step on a refined grid is less than that on the coarser grid, by the refinement factor for the spatial dimension. More time steps are taken on finer grids, and the advance of the flow solution is synchronized by interleaving the integrations at different levels. This helps prevent any interlevel mismatches that could adversely affect the accuracy of the simulation.

In order for the AMR method to sense where further refinement is needed, some monitoring function is required. We have used a combination of the gradient of the laboratory frame mass density, a test that recognizes the presence of contact surfaces, and a measure of the cell-to-cell shear in the flow, the choice of which functions to use being determined by which part of the flow is of most significance in a given study. Since the tracking of shock waves is of paramount importance, a buffer ensures the flagging of extra cells at the edge of meshes, ensuring that important flow structures do not `leak' out of meshes during the update of the hydrodynamic solution. The combined effect of using the RHLLE single mesh solver and the AMR algorithm results in a very efficient scheme. Where the RHLLE method is unable to give adequate resolution on a single coarse mesh the AMR algorithm places more cells, resulting in an excellent overall coverage of the computational domain.

next up previous
Next: Vectorization or Parallelization Implementation Up: Performance Report Previous: Performance Report

G. Comer Duncan
Sat. April 3, 1999