3.2 Other techniques3 Numerical Schemes3 Numerical Schemes

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 tex2html_wrap_inline3757 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  [116Jump To The Next Citation Point In The Article]). 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 tex2html_wrap_inline3777, an algorithm written in conservation form reads:


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


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 tex2html_wrap_inline3795, of the ``viscous version'' of Eq. (48Popup Equation):


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 tex2html_wrap_inline3797, TV tex2html_wrap_inline3799, is defined as


A numerical scheme is said to be TV-stable if TV tex2html_wrap_inline3799 is bounded for all tex2html_wrap_inline3507 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 [117Jump To The Next Citation Point In The Article]. 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 tex2html_wrap_inline3805, with tex2html_wrap_inline3807 and tex2html_wrap_inline3809 two non-negative integers, is said to be monotone if all coefficients tex2html_wrap_inline3811 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 [223Jump To The Next Citation Point In The Article] 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 tex2html_wrap_inline3781 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. tex2html_wrap_inline3819 .

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


with tex2html_wrap_inline3821, v being the fluid velocity, tex2html_wrap_inline3575 the density, tex2html_wrap_inline3781 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 [132Jump To The Next Citation Point In The Article] 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 [226Jump To The Next Citation Point In The Article, 95Jump To The Next Citation Point In The Article].

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 [125Jump To The Next Citation Point In The Article] it has been gradually demonstrated [72, 64Jump To The Next Citation Point In The Article, 184Jump To The Next Citation Point In The Article, 68, 17Jump To The Next Citation Point In The Article, 224Jump To The Next Citation Point In The Article, 178Jump To The Next Citation Point In The Article] 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. (37Popup Equation). Let us consider a single computational cell of our discrete spacetime. Let tex2html_wrap_inline3833 be a region (simply connected) of the four-dimensional manifold tex2html_wrap_inline3835, bounded by a closed three-dimensional surface tex2html_wrap_inline3837 . We take the 3-surface tex2html_wrap_inline3837 as the standard-oriented hyper-parallelepiped made up of two spacelike surfaces tex2html_wrap_inline3841 plus timelike surfaces tex2html_wrap_inline3843 that join the two temporal slices together. By integrating system (37Popup Equation) over a domain tex2html_wrap_inline3833 of a given spacetime, the variation in time of the state vector tex2html_wrap_inline3757 within tex2html_wrap_inline3833 is given - keeping apart the source terms - by the fluxes tex2html_wrap_inline3743 through the boundary tex2html_wrap_inline3837 . The integral form of system (37Popup Equation) is


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





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. (54Popup Equation) 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 tex2html_wrap_inline3855 during a time step. At the cell-interfaces the flow conditions can be discontinuous and, following the seminal idea of Godunov [84Jump To The Next Citation Point In The Article], 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.


Click on thumbnail to view image

Figure 2: Godunov's scheme: local solutions of Riemann problems. At every interface, tex2html_wrap_inline3495, tex2html_wrap_inline3497, and tex2html_wrap_inline3499, 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 tex2html_wrap_inline3501 these discontinuities decay into three elementary waves which propagate the solution forward to the next time level tex2html_wrap_inline3503 (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 tex2html_wrap_inline3505 in tex2html_wrap_inline3507 .

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 [126Jump To The Next Citation Point In The Article] 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 tex2html_wrap_inline3871 and tex2html_wrap_inline3873 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 [219Jump To The Next Citation Point In The Article] 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 [52Jump To The Next Citation Point In The Article] 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, tex2html_wrap_inline3875 are, respectively, the eigenvalues and right-eigenvectors of the Jacobian matrix, and quantities tex2html_wrap_inline3877 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 [64Jump To The Next Citation Point In The Article], as discussed in Section  2.2.1, has explicitly derived a relativistic Roe Riemann solver [183Jump To The Next Citation Point In The Article], Schneider et al. [190] carried out the extension of Einfeldt's HLLE method [90, 61], Martí and Müller [128Jump To The Next Citation Point In The Article] 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. [56Jump To The Next Citation Point In The Article] extended Marquina's method [57Jump To The Next Citation Point In The Article]. The interested reader is referred to [126] for a comprehensive description of such solvers in special relativistic hydrodynamics.

3.2 Other techniques3 Numerical Schemes3 Numerical Schemes

image Numerical Hydrodynamics in General Relativity
José A. Font
© Max-Planck-Gesellschaft. ISSN 1433-8351
Problems/Comments to livrev@aei-potsdam.mpg.de