#equation33#
the conservative form of the relativistic Euler equation is
#equation36#
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
#equation45#
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_inline177# and #tex2html_wrap_inline179# 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_inline185# 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:
#equation65#
#equation67#
#equation71#
#equation75#
#equation79#
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:
#eqnarray88#
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_inline215#. 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.