## 3.1 Finite difference schemes

Any system of equations presented in the previous section can be solved numerically by replacing the partial derivatives by finite differences on a discrete numerical grid, and then advancing the solution in time via some time-marching algorithm. Hence, specification of on an initial hypersurface, together with a suitable EOS, followed by a recovery of the primitive variables, leads to the computation of the fluxes and source terms. In doing so, the first time derivative of the data is obtained, which then leads to the formal propagation of the solution forward in time, with a time-step constrained by the Courant-Friedrichs-Lewy (CFL) condition.

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.

### 3.1.1 Artificial viscosity schemes

High-resolution schemes based on the idea of modifying the equations by introducing some terms providing artificial viscosity to damp the amplitude of spurious oscillations near discontinuities were originally proposed by von Neumann and Richtmyer [223] to solve the (classical) Euler equations. The basic idea is to introduce a purely artificial dissipative mechanism whose form and strength are such that the shock transition is smooth, extending over a small number of intervals of the space variable. In this approach the equations are discretised using a high order finite difference method (e.g., a second order Lax-Wendroff scheme), which gives good results in smooth parts of the flow solution, adding an artificial viscosity term Q to the hyperbolic system which only acts where discontinuities arise. The term Q must vanish as the grid is refined and the time step is reduced, i.e. .

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].

### 3.1.2 High-resolution shock-capturing (HRSC) schemes

Since [125] it has been gradually demonstrated [72, 64, 184, 68, 17, 224, 178] that conservative methods exploiting the hyperbolic character of the hydrodynamic equations are optimally suited for accurate integrations, even well inside the ultrarelativistic regime. The explicit knowledge of the characteristic fields (eigenvalues) of the equations, together with the corresponding eigenvectors, provides the mathematical (and physical) framework for such integrations, by means of either exact or approximate Riemann solvers. As further discussed below these solvers compute, at every interface of the numerical grid, the solution of local Riemann problems (i.e., the simplest initial value problem with discontinuous initial data). Hence, a HRSC scheme automatically guarantees that the physical discontinuities appearing in the solution, e.g. shock waves, will be treated consistently (the shock-capturing property). HRSC schemes also produce stable and sharp discrete shock profiles, while providing a high order of accuracy, typically second order or more, in smooth parts of the solution.

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.

Figure 2: Godunov's scheme: local solutions of Riemann problems. At every interface, , , and , a local Riemann problem is set up as a result of the discretization process (bottom panel), when approximating the numerical solution by piecewise constant data. At time these discontinuities decay into three elementary waves which propagate the solution forward to the next time level (top panel). The time step of the numerical scheme must satisfy the Courant-Friedrichs-Lewy condition, i.e. be small enough to prevent the waves from advancing more than in .

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