Finite difference numerical schemes provide solutions of the
discretized version of the original system of partial
differential equations. Therefore, convergence properties under
grid refinement must be enforced on such schemes to ensure that
the numerical results are correct (i.e. the global error of the
numerical solution has to tend to zero as the cell width tends to
zero). For hyperbolic systems of conservation laws, schemes
written in
*conservation form*
are preferred as they guarantee that the convergence, if it
exists, is to one of the
*weak solutions*
of the original system of equations (*Lax-Wendroff theorem*
[116]). Such weak solutions are generalized solutions that satisfy
the integral form of the conservation system. They are classical
solutions (continuous and differentiable) in regions where they
are continuous and have a finite number of discontinuities.

Let us consider an initial value problem for a one-dimensional scalar hyperbolic conservation law,

Introducing a discrete numerical grid of space-time points , an algorithm written in conservation form reads:

where
and
are the time-step and cell width respectively,
is a consistent numerical flux function (i.e.,
) and
is an approximation to the average of
*u*
(*x*,
*t*) within the numerical cell
:

The class of all weak solutions is too wide in the sense that
there is no uniqueness for the initial value problem. The
numerical method should guarantee convergence to the
*physically admissible solution*
. This is the vanishing-viscosity limit solution, i.e., the
solution when
, of the ``viscous version'' of Eq. (48):

Mathematically, this solution is characterized by the
so-called
*entropy condition*
(in the language of fluids, the condition that the entropy of
any fluid element should increase when running into a
discontinuity). The characterization of the
*entropy-satisfying solutions*
for scalar equations was given by Oleinik [162]. For hyperbolic systems of conservation laws it was developed
by Lax [115].

The Lax-Wendroff theorem [116] cited above does not establish whether the method converges. To
guarantee convergence, some form of stability is required, as Lax
proposed for linear problems (*Lax equivalence theorem*
; see, e.g., [182]). Along this direction, the notion of total-variation stability
has proven very successful although powerful results have only
been obtained for scalar conservation laws. The total variation
of a solution at
, TV
, is defined as

A numerical scheme is said to be TV-stable if TV
is bounded for all
at any time for each initial data. In the case of non-linear
scalar conservation laws it can be proved that for numerical
schemes in conservation form with consistent numerical flux
functions, TV-stability is a sufficient condition for
convergence [117]. Current research has focused on the development of
high-resolution numerical schemes in conservation form satisfying
the condition of TV-stability, e.g., the
*total variation diminishing*
(TVD) schemes. This issue is further discussed in Section
3.1.2
. Additionally, an important property that a numerical method
must satisfy is to be monotone. A scheme of the form
, with
and
two non-negative integers, is said to be monotone if all
coefficients
are positive or zero. It has been shown that, for scalar
conservation laws, monotone methods are TVD and satisfy a
discrete entropy condition. Therefore, they converge in a
non-oscillatory manner to the unique entropy solution. However,
monotone methods are at most first order accurate [117].

The hydrodynamic equations constitute a non-linear hyperbolic system and, hence, smooth initial data can give rise to discontinuous data (crossing of characteristics in the case of shocks) in a finite time during the evolution. As a consequence classical finite difference schemes present important deficiencies when dealing with such systems. Typically, first order accurate schemes are too dissipative across discontinuities (excessive smearing), and second order (or higher) schemes produce spurious oscillations near discontinuities which do not disappear as the grid is refined. Standard finite difference schemes have been conveniently modified in two ways to obtain high-order, oscillation-free accurate representations of discontinuous solutions as we discuss next.

Von Neumann and Richtmyer derived the following expression for the viscosity term:

with
,
*v*
being the fluid velocity,
the density,
the spatial interval, and
*k*
a constant parameter whose value is conveniently adjusted in
every numerical experiment. This parameter controls the number of
zones in which shock waves are spread.

This type of recipe, with minor modifications, has been used in all numerical simulations employing May and White's formulation, mostly in the context of gravitational collapse, as well as Wilson's formulation. So, for example, in May and White's code [132] the artificial viscosity term, obtained in analogy with the one originally proposed by von Neumann and Richtmyer [223], is introduced in the equations accompanying the pressure, in the form

Further examples of equivalent expressions for the artificial viscosity terms, in the context of Wilson formulation, can e.g. be found in [226, 95].

The main advantages of the artificial viscosity approach are:
(i) It is straightforward to implement (compared to the HRSC
schemes which need the characteristic fields of the equations),
and (ii) it is computationally very efficient. Experience has
shown, however, that this procedure is (i) problem dependent and
(ii) inaccurate for ultrarelativistic flows [159]. Additionally, the artificial viscosity approach has the
implicit difficulty of finding the appropriate form for
*Q*
that introduces the necessary amount of dissipation to reduce
the spurious oscillations and, at the same time, avoids
introducing excessive smearing in the discontinuities. In many
instances both properties are difficult to achieve
simultaneously. A comprehensive numerical study of artificial
viscosity induced errors in strong shock calculations in
Newtonian hydrodynamics (including also proposed improvements)
was presented by Noh [158].

To focus the discussion let us consider the system as formulated in Eq. (37). Let us consider a single computational cell of our discrete spacetime. Let be a region (simply connected) of the four-dimensional manifold , bounded by a closed three-dimensional surface . We take the 3-surface as the standard-oriented hyper-parallelepiped made up of two spacelike surfaces plus timelike surfaces that join the two temporal slices together. By integrating system (37) over a domain of a given spacetime, the variation in time of the state vector within is given - keeping apart the source terms - by the fluxes through the boundary . The integral form of system (37) is

which can be written in the following conservation form, well-adapted to numerical applications:

where

An important property of writing a numerical scheme in conservation form is that, in the absence of sources, the (physically) conserved quantities, according to the partial differential equations, are numerically conserved by the finite difference equations.

The ``hat'' symbol appearing on the fluxes of Eq. (54) indicates the numerical fluxes. These are recognized as approximations to the time-averaged fluxes across the cell-interfaces, which depend on the solution at those interfaces during a time step. At the cell-interfaces the flow conditions can be discontinuous and, following the seminal idea of Godunov [84], the numerical fluxes can be obtained by solving a collection of local Riemann problems. This is depicted in Fig. 2 . The continuous solution is locally averaged on the numerical grid, a process which leads to the appearance of discontinuities at the cell-interfaces. Physically, every discontinuity decays into three elementary waves: a shock wave, a rarefaction wave and a contact discontinuity. The complete structure of the Riemann problem can be solved analytically (see [84] for the solution in Newtonian hydrodynamics and [127] in special relativistic hydrodynamics) and, accordingly, used to update the solution forward in time.

For reasons of efficiency and, particularly, in multidimensions, the exact solution of the Riemann problem is frequently avoided and linearized (approximate) Riemann solvers are preferred. These solvers are based on the exact solution of Riemann problems corresponding to a linearized version of the original system of equations and, after extensive experimentation, they are found to produce results comparable to those obtained with the exact solver (see [126] for a summary of such approximate solvers in special relativistic hydrodynamics).

In the frame of the local characteristic approach the numerical fluxes are computed according to some generic flux-formula which makes use of the characteristic information of the system. For example, in Roe's approximate Riemann solver it adopts the following functional form:

where and represent the values of the primitive variables at the left and right sides, respectively, of the corresponding interface. They are obtained from the cell-centered quantities after a suitable monotone reconstruction procedure. The way these variables are computed determines the spatial order of the numerical algorithm and controls the local jumps at every interface. If these jumps are monotonically reduced, the scheme provides more accurate initial guesses for the solution of the local Riemann problems. A wide variety of cell reconstruction procedures is available in the literature. Among the most popular slope limiter procedures for TVD schemes [88] are the second order piecewise linear reconstruction, introduced by van Leer [219] in the design of the MUSCL scheme (Monotonic Upstream Scheme for Conservation Laws), and the third order piecewise parabolic reconstruction developed by Colella and Woodward [52] in their Piecewise Parabolic Method (PPM). High order piecewise polynomial functions are also available for Essentially Non-Oscillatory (ENO) schemes [89].

The last term in the flux-formula represents the numerical viscosity of the scheme, and it makes explicit use of the characteristic information of the Jacobian matrices of the system. This information is used to provide the appropriate amount of numerical dissipation to obtain accurate representations of discontinuous solutions without excessive smearing, avoiding, at the same time, the growth of spurious numerical oscillations associated with the Gibbs phenomenon. Hence, are, respectively, the eigenvalues and right-eigenvectors of the Jacobian matrix, and quantities are the jumps of the characteristic variables across each characteristic field. They are obtained by projecting the jumps of the state-vector variables with the left-eigenvectors matrix:

The ``tilde'' indicates that the corresponding fields are averaged at the cell interfaces from the left and right (reconstructed) values.

During the last few years most of the
*classical*
Riemann solvers developed in fluid dynamics have been extended
to relativistic hydrodynamics: Eulderink [64], as discussed in Section
2.2.1, has explicitly derived a relativistic Roe Riemann solver [183], Schneider et al. [190] carried out the extension of Einfeldt's HLLE method [90,
61], Martí and Müller [128] extended the PPM method of Woodward and Colella [233], Wen et al. [224] extended Glimm's method, Dolezal and Wong [55] put into practice Shu-Osher ENO techniques, Balsara [16] extended Colella's two-shock approximation and, Donat et
al. [56] extended Marquina's method [57]. The interested reader is referred to [126] for a comprehensive description of such solvers in special
relativistic hydrodynamics.

Numerical Hydrodynamics in General Relativity
José A. Font
http://www.livingreviews.org/lrr-2000-2
© Max-Planck-Gesellschaft. ISSN 1433-8351 Problems/Comments to livrev@aei-potsdam.mpg.de |