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 the state vector tex2html_wrap_inline4902 on an initial hypersurface, together with a suitable choice of EOS, followed by a recovery of the primitive variables, leads to the computation of the fluxes and source terms. Through this procedure 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.

The hydrodynamic equations (either in Newtonian physics or in general relativity) constitute a nonlinear hyperbolic system and, hence, smooth initial data can transform into discontinuous data (the crossing of characteristics in the case of shocks) in a finite time during the evolution. As a consequence, classical finite difference schemes (see, e.g., [152Jump To The Next Citation Point In The Article, 287Jump To The Next Citation Point In The Article]) present important deficiencies when dealing with such systems. Typically, first-order accurate schemes are much 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. To avoid these effects, standard finite difference schemes have been conveniently modified in various ways to ensure high-order, oscillation-free accurate representations of discontinuous solutions, as we discuss next.

3.1.1 Artificial viscosity approach 

The idea of modifying the hydrodynamic equations by introducing artificial viscosity terms to damp the amplitude of spurious oscillations near discontinuities was originally proposed by von Neumann and Richtmyer [295Jump To The Next Citation Point In The Article] in the context of 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 becomes smooth, extending over a small number of intervals tex2html_wrap_inline4924 of the space variable. In their original work, von Neumann and Richtmyer proposed the following expression for the viscosity term:


with tex2html_wrap_inline4926, v being the fluid velocity, tex2html_wrap_inline4640 the density, tex2html_wrap_inline4924 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 generic recipe, with minor modifications, has been used in conjuction with standard finite difference schemes 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 original code [172Jump To The Next Citation Point In The Article] the artificial viscosity term, obtained in analogy with the one proposed by von Neumann and Richtmyer [295], is introduced in the equations accompanying the pressure, in the form


Further examples of similar expressions for the artificial viscosity terms, in the context of Wilson formulation, can be found in, e.g., [300Jump To The Next Citation Point In The Article, 123Jump To The Next Citation Point In The Article]. A state-of-the-art formulation of the artificial viscosity approach is reported in [13Jump To The Next Citation Point In The Article].

The main advantage of the artificial viscosity approach is its simplicity, which results in high computational efficiency. Experience has shown, however, that this procedure is both problem dependent and inaccurate for ultrarelativistic flows [208Jump To The Next Citation Point In The Article, 13Jump To The Next Citation Point In The Article]. Furthermore, the artificial viscosity approach has the inherent ambiguity 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 [207].

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

In finite difference schemes, convergence properties under grid refinement must be enforced to ensure that the numerical results are correct (i.e., if a scheme with an order of accuracy tex2html_wrap_inline4738 is used, the global error of the numerical solution has to tend to zero as tex2html_wrap_inline4940 as the cell width tex2html_wrap_inline4924 tends to zero). For hyperbolic systems of conservation laws, schemes written in conservation form are preferred since, according to the Lax-Wendroff theorem [150Jump To The Next Citation Point In The Article], they guarantee that the convergence, if it exists, is to one of the so-called weak solutions of the original system of equations. Such weak solutions are generalized solutions that satisfy the integral form of the conservation system. They are tex2html_wrap_inline4944 classical solutions (continuous and differentiable) in regions where they are continuous and have a finite number of discontinuities.

For the sake of simplicity let us consider in the following an initial value problem for a one-dimensional scalar hyperbolic conservation law,


and let us introduce a discrete numerical grid of space-time points tex2html_wrap_inline4946 . An explicit algorithm written in conservation form updates the solution from time tex2html_wrap_inline4628 to the next time level tex2html_wrap_inline4630 as


where tex2html_wrap_inline4952 is a consistent numerical flux function (i.e., tex2html_wrap_inline4954) of p + q +1 arguments and tex2html_wrap_inline4634 and tex2html_wrap_inline4924 are the time step and cell width respectively. Furthermore, tex2html_wrap_inline4962 is an approximation of the average of u (x, t) within the numerical cell tex2html_wrap_inline4966 tex2html_wrap_inline4968 :


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, hence, guarantee convergence to the physically admissible solution . This is the vanishing-viscosity limit solution, that is, the solution when tex2html_wrap_inline4970, of the ``viscous version'' of the scalar conservation law, Equation (39Popup Equation):


Mathematically, the solution one needs to search for 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 these entropy-satisfying solutions for scalar equations was given by Oleinik [212]. For hyperbolic systems of conservation laws it was developed by Lax [149].

The Lax-Wendroff theorem [150] cited above does not establish whether the method converges. To guarantee convergence, some form of stability is required, as Lax first proposed for linear problems (Lax equivalence theorem ; see, e.g., [241]). 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 time tex2html_wrap_inline4972, TV tex2html_wrap_inline4974, is defined as


A numerical scheme is said to be TV-stable if TV tex2html_wrap_inline4974 is bounded for all tex2html_wrap_inline4634 at any time for each initial data. In the case of nonlinear, scalar conservation laws it can be proved that TV-stability is a sufficient condition for convergence [152], as long as the numerical schemes are written in conservation form and have consistent numerical flux functions. Current research has focused on the development of high-resolution numerical schemes in conservation form satisfying the condition of TV-stability, such as the so-called total variation diminishing (TVD) schemes [115Jump To The Next Citation Point In The Article] (see below).

Let us now consider the specific system of hydrodynamic equations as formulated in Equation (30Popup Equation), and let us consider a single computational cell of our discrete spacetime. Let tex2html_wrap_inline4980 be a region (simply connected) of a given four-dimensional manifold tex2html_wrap_inline4982, bounded by a closed three-dimensional surface tex2html_wrap_inline4984 . We further take the 3-surface tex2html_wrap_inline4984 as the standard-oriented hyper-parallelepiped made up of two spacelike surfaces tex2html_wrap_inline4988 plus timelike surfaces tex2html_wrap_inline4990 that join the two temporal slices together. By integrating system (30Popup Equation) over a domain tex2html_wrap_inline4980 of a given spacetime, the variation in time of the state vector tex2html_wrap_inline4902 within tex2html_wrap_inline4980 is given - keeping apart the source terms - by the fluxes tex2html_wrap_inline4888 through the boundary tex2html_wrap_inline4984 . The integral form of system (30Popup Equation) is


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





A numerical scheme written in conservation form ensures 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 computation of the time integrals of the interface fluxes appearing in Equation (45Popup Equation) is one of the distinctive features of upwind HRSC schemes. One needs first to define the so-called numerical fluxes, which are recognized as approximations to the time-averaged fluxes across the cell interfaces, which depend on the solution at those interfaces, tex2html_wrap_inline5020, during a time step,


At the cell interfaces the flow can be discontinuous and, following the seminal idea of Godunov [108Jump To The Next Citation Point In The Article], the numerical fluxes can be obtained by solving a collection of local Riemann problems, as is depicted in Figure  2 . This is the approach followed by the so-called Godunov-type methods [117Jump To The Next Citation Point In The Article, 75Jump To The Next Citation Point In The Article]. Figure  2 shows how the continuous solution is locally averaged on the numerical grid, a process that 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 [108] for the solution in Newtonian hydrodynamics and [165Jump To The Next Citation Point In The Article, 231Jump To The Next Citation Point In The Article] 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,, 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 tex2html_wrap_inline4628 these discontinuities decay into three elementary waves, which propagate the solution forward to the next time level tex2html_wrap_inline4630 (top panel). The time step of the numerical scheme must satisfy the Courant-Friedrichs-Lewy condition, being small enough to prevent the waves from advancing more than tex2html_wrap_inline4632 in tex2html_wrap_inline4634 .

For reasons of numerical efficiency and, particularly in multi-dimensions, 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. After extensive experimentation, the results achieved with approximate Riemann solvers are comparable to those obtained with the exact solver (see [287Jump To The Next Citation Point In The Article] for a comprehensive overview of Godunov-type methods, and [164Jump To The Next Citation Point In The Article] for an excellent summary of approximate Riemann solvers in special relativistic hydrodynamics).

In the frame of the local characteristic approach, the numerical fluxes appearing in Equation (45Popup Equation) are computed according to some generic flux-formula that makes use of the characteristic information of the system. For example, in Roe's approximate Riemann solver [242Jump To The Next Citation Point In The Article] it adopts the following functional form:


where tex2html_wrap_inline5036 and tex2html_wrap_inline5038 are the values of the primitive variables at the left and right sides, respectively, of a given cell 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 accuracy of the numerical algorithm and controls the amplitude of the local jumps at every cell interface. If these jumps are monotonically reduced, the scheme provides more accurate initial guesses for the solution of the local Riemann problems (the average values give only first-order accuracy). A wide variety of cell reconstruction procedures is available in the literature. Among the slope limiter procedures (see, e.g., [287Jump To The Next Citation Point In The Article, 153Jump To The Next Citation Point In The Article]) most commonly used for TVD schemes [115] are the second order, piecewise-linear reconstruction, introduced by van Leer [290Jump 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 [58Jump To The Next Citation Point In The Article] in their Piecewise Parabolic Method (PPM). Since TVD schemes are only first-order accurate at local extrema, alternative reconstruction procedures for which some growth of the total variation is allowed have also been developed. Among those, we mention the total variation bounded (TVB) schemes [268] and the essentially non-oscillatory (ENO) schemes [116].

Alternatively, high-order methods for nonlinear hyperbolic systems have also been designed using flux limiters rather than slope limiters, as in the FCT scheme of Boris and Book [46Jump To The Next Citation Point In The Article]. In this approach, the numerical flux consists of two pieces, a high-order flux (e.g., the Lax-Wendroff flux) for smooth regions of the flow, and a low-order flux (e.g., the flux from some monotone method) near discontinuities, tex2html_wrap_inline5040 with the limiter tex2html_wrap_inline5042 (see [287Jump To The Next Citation Point In The Article, 153Jump To The Next Citation Point In The Article] for further details).

The last term in the flux-formula, Equation (49Popup Equation), 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. In Equation (49Popup Equation), tex2html_wrap_inline5044 are the eigenvalues and right-eigenvectors of the Jacobian matrix of the flux vector, respectively. Correspondingly, the quantities tex2html_wrap_inline5046 are the jumps of the so-called 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'' in Equations (49Popup Equation) and (50Popup Equation) 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 standard Riemann solvers developed in Newtonian fluid dynamics have been extended to relativistic hydrodynamics: Eulderink [78Jump To The Next Citation Point In The Article], as discussed in Section  2.2.1, explicitly derived a relativistic Roe Riemann solver [242Jump To The Next Citation Point In The Article]; Schneider et al. [250] carried out the extension of Harten, Lax, van Leer, and Einfeldt's (HLLE) method [117, 75]; Martí and Müller [166Jump To The Next Citation Point In The Article] extended the PPM method of Woodward and Colella [306]; Wen et al. [297] extended Glimm's exact Riemann solver; Dolezal and Wong [68] put into practice Shu-Osher ENO techniques; Balsara [19] extended Colella's two-shock approximation, and Donat et al. [69Jump To The Next Citation Point In The Article] extended Marquina's method [70Jump To The Next Citation Point In The Article]. Recently, much effort has been spent concerning the exact special relativistic Riemann solver and its extension to multi-dimensions [165, 231, 238, 239Jump To The Next Citation Point In The Article]. The interested reader is addressed to [164] and references therein for a comprehensive description of such solvers in special relativistic hydrodynamics.

3.1.3 High-order central schemes 

The use of high-order non-oscillatory symmetric (central) TVD schemes for solving hyperbolic systems of conservation laws emerged at the mid 1980s [61Jump To The Next Citation Point In The Article, 243Jump To The Next Citation Point In The Article, 311Jump To The Next Citation Point In The Article, 205Jump To The Next Citation Point In The Article] (see also [312] and [287Jump To The Next Citation Point In The Article] and references therein) as an alternative approach to upwind HRSC schemes. Symmetric schemes are based on either high-order schemes (e.g., Lax-Wendroff) with additional dissipative terms [61Jump To The Next Citation Point In The Article, 243, 311], or on non-oscillatory high-order extensions of first-order central schemes (e.g., Lax-Friedrichs) [205]. One of the nicest properties of central schemes is that they exploit the conservation form of the Lax-Wendroff or Lax-Friedrichs schemes. Therefore, they yield the correct propagation speeds of all nonlinear waves appearing in the solution. Furthermore, central schemes sidestep the use of Riemann solvers, which results in enhanced computational efficiency, especially in multi-dimensional problems. Its use is, thus, not limited to those systems of equations where the characteristic information is explicitly known or to systems where the solution of the Riemann problem is prohibitively expensive to compute. This approach has gradually developed during the last decade to reach a mature status where a number of straightforward central schemes of high order can be applied to any nonlinear hyperbolic system of conservation laws. The typical results obtained for the Euler equations show a quality comparable to that of upwind HRSC schemes, at the expense of a small loss of sharpness of the solution at discontinuities [287Jump To The Next Citation Point In The Article]. An up-to-date summary of the status and applications of this approach is discussed in [287Jump To The Next Citation Point In The Article, 145, 281].

In the context of special and general relativistic MHD, Koide et al. [141Jump To The Next Citation Point In The Article, 142Jump To The Next Citation Point In The Article] applied a second-order central scheme with nonlinear dissipation developed by Davis [61Jump To The Next Citation Point In The Article] to the study of black hole accretion and formation of relativistic jets. One-dimensional test simulations in special relativistic hydrodynamics performed by the author and coworkers [92Jump To The Next Citation Point In The Article] using the SLIC (slope limiter centred) scheme (see [287Jump To The Next Citation Point In The Article] for details) showed its capabilities to yield satisfactory results, comparable to those obtained by HRSC schemes based on Riemann solvers, even well inside the ultrarelativistic regime. The slopes of the original central scheme were limited using high-order reconstruction procedures such as PPM [58Jump To The Next Citation Point In The Article], which was essential to keep the inherent diffusion of central schemes at discontinuities at reasonable levels. Very recently, Del Zanna and Bucciantini [63Jump To The Next Citation Point In The Article] assessed a third-order convex essentially non-oscillatory central scheme in multi-dimensional special relativistic hydrodynamics. Again, these authors obtained results as accurate as those of upwind HRSC schemes in standard tests (shock tubes, shock reflection test). Yet another central scheme has been assessed by [13Jump To The Next Citation Point In The Article] in one-dimensional special and general relativistic hydrodynamics, where similar results to those of [63Jump To The Next Citation Point In The Article] are reported. These authors also validate their central scheme in simulations of spherical accretion onto a Schwarzschild black hole, and further provide a comparison with an artificial viscosity based scheme.

It is worth emphasizing that early pioneer approaches in the field of relativistic hydrodynamics [208, 54] used standard finite difference schemes with artificial viscosity terms to stabilize the solution at discontinuities. However, as discussed in Section  2.1.2, those approaches only succeeded in obtaining accurate results for moderate values of the Lorentz factor (tex2html_wrap_inline5048). A key feature lacking in those investigations was the use of a conservative approach for both the system of equations and the numerical schemes. A posteriori, and in the light of the results reported by [63, 13, 92], it appears that this was the ultimate reason preventing the extension of the early computations to the ultrarelativistic regime.

The alternative of using high-order central schemes instead of upwind HRSC schemes becomes apparent when the spectral decomposition of the hyperbolic system under study is not known. The straightforwardness of a central scheme makes its use very appealing, especially in multi-dimensions where computational efficiency is an issue. Perhaps the most important example in relativistic astrophysics is the system of (general) relativistic MHD equations. Despite some progress in recent years (see, e.g., [20, 143]), much more work is needed concerning their solution with HRSC schemes based upon Riemann solvers. Meanwhile, an obvious choice is the use of central schemes [141Jump To The Next Citation Point In The Article, 142Jump To The Next Citation Point In The Article].

3.1.4 Source terms 

Most ``conservation laws'' involve source terms on the right hand side of the equations. In hydrodynamics, for instance, those terms arise when considering external forces such as gravity, which make the right hand side of the momentum and energy equations no longer zero (see Section  2). Other effects leading to the appearance of source terms are radiative heat transfer (accounted for in the energy equation) or ionization (resulting in a collection of non-homogeneous continuity equations for the mass of each species, which is not conserved separately). The incorporation of the source terms in the solution procedure is a common issue to all numerical schemes considered so far. Since a detailed discussion on the numerical treatment of source terms is beyond the scope of this article, we simply provide some basic information in this section, addressing the interested reader to [287, 153Jump To The Next Citation Point In The Article] and references therein for further details.

There are, essentially, two basic ways of handling source terms. The first approach is based on unsplit methods by which a single finite difference formula advances the entire equation over one time step (as in Equation (45Popup Equation)):


The temporal order of this basic algorithm can be improved by introducing successive sub-steps to perform the time update (e.g., predictor-corrector, Shu and Osher's conservative high order Runge-Kutta schemes, etc.)

Correspondingly, the second approach is based on fractional step or splitting methods . The basic idea is to split the equation into different pieces (transport + sources) and apply appropriate methods for each piece independently. In the first-order Godunov splitting, tex2html_wrap_inline5050, the operator tex2html_wrap_inline5052 solves the homogeneous PDE in the first step to yield the intermediate value tex2html_wrap_inline5054 . Then, in the second step, the operator tex2html_wrap_inline5056 solves the ordinary differential equation tex2html_wrap_inline5058 to yield the final value tex2html_wrap_inline5060 . In order to achieve second-order accuracy (assuming each independent method is second order), a common fractional step method is the Strang splitting, where tex2html_wrap_inline5062 . Therefore, this method advances by a half time step the solution for the ODE containing the source terms, and by a full time step the conservation law (the PDE) in between each source step.

We note that in some cases the source terms may become stiff, as in phenomena that either occur on much faster timescales than the hydrodynamic time step tex2html_wrap_inline4634, or act over much smaller spatial scales than the grid resolution tex2html_wrap_inline4924 . Stiff source terms may easily lead to numerical difficulties. The interested reader is directed to [153] (and references therein) for further information on various approaches to overcome the problems of stiff source terms.

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