Go to previous page Go up Go to next page

4.1 Finite difference schemes

Any system of equations presented in the previous sections 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 U 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 timestep constrained by the Courant–Friedrichs–Lewy (CFL) condition.

The hydrodynamic and MHD equations (either in Newtonian physics or in general relativity) constitute nonlinear hyperbolic systems and, hence, smooth initial data can transform into discontinuous data (cross of characteristics in the case of shocks) in a finite time during the evolution. As a consequence, conventional finite-difference schemes (see, e.g., [218Jump To The Next Citation Point219Jump To The Next Citation Point398Jump To The Next Citation Point]) 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 under grid refinement. 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.

4.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 [407Jump To The Next Citation Point] 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 Δx of the space variable. In their original work von Neumann and Richtmyer proposed the following expression for the viscosity term:

{ − α ∂vif ∂v< 0 or ∂ρ > 0, Q = ∂x ∂x ∂t 0 otherwise,

with α = ρ (kΔx )2∂∂vx, v being the fluid velocity, ρ the density, Δx the spatial interval, and k a constant parameter whose value needs to be 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 conjunction 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 [247Jump To The Next Citation Point] the artificial viscosity term, obtained in analogy with the one proposed by von Neumann and Richtmyer [407], is introduced in the equations accompanying the pressure, in the form:

{ 2 ρ(aΔRm2-)2∂R∂muβˆ•Γ if ∂∂ρt > 0, Q = 0 otherwise.

Further examples of similar expressions for the artificial viscosity terms, in the context of Wilson formulation, can be found in, e.g., [411Jump To The Next Citation Point171Jump To The Next Citation Point417]. In particular, a state-of-the-art formulation of the artificial viscosity approach is reported in [19Jump To The Next Citation Point]. Correspondingly, the interested reader is also directed to the recent work by [86Jump To The Next Citation Point20Jump To The Next Citation Point108Jump To The Next Citation Point] for details on diverse modern implementations of artificial viscosity terms in nonconservative formulations of the GRMHD equations.

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 [291Jump To The Next Citation Point19Jump To The Next Citation Point]. 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 at 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 [290].

4.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 α is used, the global error of the numerical solution has to tend to zero as π’ͺ (Δx )α as the cell width Δx tends to zero). For hyperbolic systems of conservation laws, schemes written in conservation form are preferred since, according to the Lax–Wendroff theorem [216Jump To The Next Citation Point], they guarantee that the convergence, if it exists, is to one of the 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 1 π’ž classical solutions (continuous and differentiable) away from discontinuities 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,

∂u- + ∂f-(u) = 0, u(x,t = 0) = u (x), (70 ) ∂t ∂x 0
and let us introduce a discrete numerical grid of spacetime points (xj,tn). An explicit algorithm written in conservation form updates the solution from time n t to the next time level n+1 t as:
n+1 n Δt n n n n n n u j = uj − Δx--(fˆ(uj−p,uj−p+1,⋅⋅⋅,uj+q) − ˆf(uj−p− 1,u j− p,⋅⋅⋅,uj+q−1)), (71 )
where ˆ f is a consistent numerical flux function (i.e., ˆ f(u,u, ⋅⋅⋅,u) = f(u)) of p + q + 1 arguments and Δt and Δx are the timestep and cell width respectively. Furthermore, unj is an approximation of the average of u(x,t) within the numerical cell [xj− 1βˆ•2,xj+1βˆ•2] (xj±1βˆ•2 = (xj + xj±1)βˆ•2):
n -1--∫ xj+1βˆ•2 n uj ≈ Δx x u(x, t )dx. (72 ) j−1βˆ•2

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 η → 0, of the “viscous version” of the scalar conservation law, Equation (70View Equation):

∂u- ∂f(u-) ∂2u- ∂t + ∂x = η∂x2 . (73 )

Mathematically, the solution one needs to search for is characterized by the entropy condition (in the language of fluids, the condition that the entropy of any fluid element should increase when crossing a discontinuity). The characterization of these entropy-satisfying solutions for scalar equations was given by Oleinik [301]. For hyperbolic systems of conservation laws it was developed by Lax [215].

The Lax–Wendroff theorem [216] 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., [336]). 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 n t = t, TVn (u ), is defined as

+∞ TV (un) = ∑ |un − un|. (74 ) j=0 j+1 j

A numerical scheme is said to be TV-stable if TV(un) is bounded for all Δt at any time for each initial data. In the case of nonlinear, scalar conservation laws it can be proven that TV-stability is a sufficient condition for convergence [218Jump To The Next Citation Point], 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 total variation diminishing (TVD) schemes [162Jump To The Next Citation Point] (see below).

Let us now consider the specific system of hydrodynamic equations as formulated in Equation (32View Equation) and let us consider a single computational cell of our discrete spacetime. Let Ω be a region (simply connected) of a given four-dimensional manifold β„³, bounded by a closed three-dimensional surface ∂ Ω. We further take the three-surface ∂Ω as the standard-oriented hyper-parallelepiped made up of two spacelike surfaces {Σ 0,Σ 0 0} x x +Δx plus timelike surfaces {Σ i,Σ i i} x x +Δx that join the two temporal slices together. By integrating system (32View Equation) over a domain Ω of a given spacetime, the variation in time of the state vector U within Ω is given – keeping apart the source terms – by the fluxes Fi through the boundary ∂ Ω. The integral form of system (32View Equation) is

∫ 1 ∂√ γU ∫ 1 ∂ √−-gFi ∫ √---- ----0-d Ω + √---------i---dΩ = Sd Ω, (75 ) Ω − g ∂x Ω − g ∂x Ω
which can be written in the following conservation form, well adapted to numerical applications:

-d-- -- dx0 (U ΔV ) = ( ∫ √ --- ∫ √ --- ) − − gF1dx2dx3 − − gF1dx2dx3 ( Σx1+ Δx1 Σx1 ) ∫ √ --- ∫ √ --- − − gF2dx1dx3 − − gF2dx1dx3 ( Σx2+ Δx2 Σx2 ) ∫ √ --- 3 1 2 ∫ √ --- 3 1 2 − Σ − gF dx dx − Σ − gF dx dx ∫ x3+ Δx3 x3 + Sd Ω, (76 ) Ω
∫ -- --1- √ -- 1 2 3 U = ΔV ΔV γUdx dx dx , (77 )
∫ x1+ Δx1∫ x2+Δx2 ∫ x3+ Δx3√ -- ΔV = γdx1dx2dx3. (78 ) x1 x2 x3

Besides its convergence properties, 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 (76View Equation) is one of the distinctive features of upwind HRSC schemes. One needs first to define the numerical fluxes, which are recognized as approximations to the time-averaged fluxes across the cell interfaces, which depend on the solution at those interfaces, j j 0 U (x + Δx βˆ•2,x ), during a timestep,

∫ n+1 ˆ -1- t j+1βˆ•2 0 Fj+ 12 ≈ Δt tn F (U (x ,x )). (79 )

At the cell interfaces the flow can be discontinuous and, following the seminal idea of Godunov [155Jump To The Next Citation Point], the numerical fluxes can be obtained by solving a collection of local Riemann problems, as is depicted in Figure 2View Image. This is the approach followed by the Godunov-type methods [164Jump To The Next Citation Point112Jump To The Next Citation Point]. Figure 2View Image 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 of any of the following type: shock waves, rarefaction waves, and contact discontinuities (see, e.g., [398Jump To The Next Citation Point]). The complete structure of the Riemann problem can be solved analytically (see [155] for the solution in Newtonian hydrodynamics, [238Jump To The Next Citation Point322Jump To The Next Citation Point] in special-relativistic hydrodynamics, and [340Jump To The Next Citation Point150Jump To The Next Citation Point] in special-relativistic MHD) and, accordingly, used to update the solution forward in time.

View Image

Figure 2: Godunov’s scheme: local solutions of Riemann problems. At every interface, xj− 1 2, xj+ 1 2 and xj+3 2, 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 n t these discontinuities decay into three elementary waves, which propagate the solution forward to the next time level tn+1 (top panel). The timestep of the numerical scheme must satisfy the Courant–Friedrichs–Lewy condition, being small enough to prevent the waves from advancing more than Δx βˆ•2 in Δt.

For reasons of numerical efficiency and, particularly in multiple 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. The spectral decomposition of the flux-vector Jacobian matrices is on the basis of all solvers (extending ideas used for linear hyperbolic systems). After extensive experimentation, the results achieved with approximate Riemann solvers are comparable to those obtained with the exact solver (see [398Jump To The Next Citation Point] for a comprehensive overview of Godunov-type methods, and [240Jump To The Next Citation Point] 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 (76View 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 [337Jump To The Next Citation Point] it adopts the following functional form:

( ∑5 ) ˆFj+ 1= 1- F (wR ) + F(wL ) − |^λn|Δ ^ωn^rn , (80 ) 2 2 n=1
where wL and wR 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.

For a purely linear system, Equation (80View Equation) provides the exact expression for the flux in terms of the conserved variables. Therefore, the above expression for the Roe flux of conserved variables is the natural extension (after linearization) of the upwind flux for characteristic variables (see, e.g., [218]), once the quasilinear system is written in diagonal form, namely

∂ ω ∂ω --- + Λ ---= 0, (81 ) ∂t ∂x
with −1 − 1 ω = R U, Λ = R AR and Λ = diag (λ1, λ2,⋅⋅⋅), where A is the Jacobian matrix of the quasilinear system, λi are its eigenvalues, and R is the right-eigenvectors matrix. The upwind flux of the characteristic variables is given by 1(λiωi L + λiωi R − |λi|Δ ωi) 2, which yields the Roe flux given by Equation (80View Equation) when written in terms of the conserved variables.

The last term in the flux formula, Equation (80View 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 (80View Equation), ^ {λn, ^rn}n=1...5 are, respectively, the eigenvalues and right eigenvectors of the Jacobian matrix of the flux vector. Correspondingly, the quantities {Δ ^ωn}n=1...5 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-eigenvector matrix:

5 U (wR ) − U (wL ) = ∑ Δ ^ωn^rn. (82 ) n=1
The “tilde” in Equations (80View Equation) and (82View Equation) indicates that the corresponding fields are averaged (local linearization) at the cell interfaces from the left and right (reconstructed) values.

The way the cell-reconstructed 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., [398Jump To The Next Citation Point219Jump To The Next Citation Point]) most commonly used for TVD schemes [162] are the second-order, piecewise-linear reconstruction, introduced by van Leer [403Jump To The Next Citation Point] 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 [80Jump To The Next Citation Point] 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 [377] and the essentially nonoscillatory (ENO) schemes [163] and extensions thereof.

Alternatively, high-order methods for nonlinear hyperbolic systems have also been designed using flux limiters rather than slope limiters, as in the Flux-Corrected Transport (FCT) scheme of Boris and Book [60Jump To The Next Citation Point]. 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, ˆ ˆ ˆ ˆ F = Fh − (1 − Φ )(Fh − Fl) with the limiter Φ ∈ [0,1] (see [398Jump To The Next Citation Point219Jump To The Next Citation Point] for further details).

During the last few years most of the standard Riemann solvers developed in Newtonian fluid dynamics have been extended to relativistic hydrodynamics: Eulderink [117Jump To The Next Citation Point], as discussed in Section 2.2.1, explicitly derived a relativistic Roe’s Riemann solver [337]; Schneider et al. [345] carried out the extension of Harten, Lax, van Leer, and Einfeldt’s (HLLE) method [164Jump To The Next Citation Point112]; Martí and Müller [239Jump To The Next Citation Point] extended the PPM method of Woodward and Colella [420]; Wen et al. [410] extended Glimm’s exact Riemann solver; Dolezal and Wong [100] put into practice Shu–Osher ENO techniques; Balsara [33] extended Colella’s two-shock approximation; Donat et al. [101Jump To The Next Citation Point] extended Marquina’s method [102Jump To The Next Citation Point]. Recently, much effort has been spent concerning the exact special relativistic Riemann solver and its extension to multiple dimensions [238322333334Jump To The Next Citation Point]. The interested reader is directed to [240Jump To The Next Citation Point] and references therein for a comprehensive description of such solvers in special-relativistic hydrodynamics.

Upwind HRSC schemes are general enough to be applicable to any hyperbolic system as long as the wave structure of the equations is known. As discussed before, the relativistic (and classical) MHD system of equations show degeneracies in the eigenvectors of the flux-vector Jacobian matrices. This fact makes it hazardous to solve them with linearized Riemann solvers based on the full spectral decomposition of the flux-vector Jacobians. For the case of special-relativistic MHD, approaches based on such full-wave decomposition have been put forward by [199Jump To The Next Citation Point34Jump To The Next Citation Point198]. Advances in this front have been significant, as even the exact solution of the Riemann problem in special-relativistic MHD has been obtained recently [340150]. In addition, there has been a successful attempt by [24Jump To The Next Citation Point] to develop a GRMHD code in which a full-wave decomposition (Roe-type) Riemann solver based on a single, renormalized set of right and left eigenvectors has been implemented. This set of eigenvectors is regular for any physical state, including degeneracies (see [23Jump To The Next Citation Point] for details). Such a Riemann solver is invoked in the code of [24Jump To The Next Citation Point] after a (local) linear coordinate transformation based on the procedure developed by [320Jump To The Next Citation Point] that allows one to use special-relativistic Riemann solvers in general relativity, and which has been properly extended to include magnetic fields (see [24Jump To The Next Citation Point] for details on such extension and Section 6.1 for details of the procedure in the purely hydrodynamic case). A similar approach is followed in the GRMHD code of Komissarov [200Jump To The Next Citation Point].

Due to the numerical degeneracies of the GRMHD system most approaches to solve those equations use, however, incomplete Riemann solvers, i.e., simpler alternative approaches to compute the numerical fluxes such as the Harten–Lax–van Leer (HLL) single-state Riemann solver of [164Jump To The Next Citation Point] (see also [174] for an improved HLL scheme at contact discontinuities (HLLC) for ideal special-relativistic MHD), or high-order central (symmetric) schemes, which entirely sidestep the use of the Jacobian matrix eigenvector structure. Such symmetric schemes, which have proven recently to yield results with an accuracy comparable to those provided by full-wave decomposition Riemann solvers in simulations involving both hydrodynamic and magnetohydrodynamic flows, are briefly discussed next.

4.1.3 High-order central schemes

The use of high-order nonoscillatory symmetric (central) TVD schemes for solving hyperbolic systems of conservation laws emerged in the mid 1980s [83Jump To The Next Citation Point338Jump To The Next Citation Point427Jump To The Next Citation Point287Jump To The Next Citation Point] (see also [428] and [398Jump To The Next Citation Point] 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 [83Jump To The Next Citation Point338427], or on nonoscillatory high-order extensions of first-order central schemes (e.g., Lax–Friedrichs) [287226Jump To The Next Citation Point]. 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 multidimensional problems. Its use is, thus, not limited to those systems of equations in which the characteristic information is explicitly known or in which the Riemann problem is prohibitively expensive to compute.

For illustrative purposes let us write the numerical flux function resulting in the fully-discrete central scheme of Kurganov and Tadmor [211Jump To The Next Citation Point], to use in the conservation form algorithm given by Equation (76View Equation):

ˆ 1-[ L R ] ai+-12[ L R] Fi+1βˆ•2 = 2 F (Ui+1) + F(U i ) − 2 U i+1 − U i . (83 )
This numerical flux depends only on the local propagation speeds ai+1βˆ•2, given by the maximum of the eigenvalues of the Jacobian matrix at the cell interface i + 1βˆ•2. Due to its simple form, it can be implemented and extended to any spatial order straightforwardly through the appropriate choice of monotone cell-reconstruction algorithms to compute the left and right states at cell interfaces.

Conservative central schemes have been gradually developed during the last few years to reach a mature status where a number of characteristic-information-free 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 [398Jump To The Next Citation Point]. An up-to-date summary of the status and applications of this approach is discussed in [398Jump To The Next Citation Point211Jump To The Next Citation Point392].

In recent years there have been various successful attempts to apply high-order central schemes to solve the relativistic hydrodynamics equations as well, including those by [196Jump To The Next Citation Point] (Lax–Wendroff scheme with conservative TVD dissipation terms), [89Jump To The Next Citation Point] (Lax–Friedrichs or HLL schemes with third-order ENO reconstruction algorithms), [19Jump To The Next Citation Point] (nonoscillatory central differencing), and [230Jump To The Next Citation Point361Jump To The Next Citation Point] (semidiscrete central scheme of Kurganov–Tadmor [211Jump To The Next Citation Point]). On the other hand, symmetric schemes (such as HLL, Kurganov–Tadmor, or Liu–Osher schemes) are being currently employed by a growing number of groups in GRMHD [196Jump To The Next Citation Point197Jump To The Next Citation Point149Jump To The Next Citation Point108Jump To The Next Citation Point366Jump To The Next Citation Point20Jump To The Next Citation Point24Jump To The Next Citation Point286Jump To The Next Citation Point].

In the context of special and general-relativistic MHD, Koide et al. [196Jump To The Next Citation Point197Jump To The Next Citation Point] applied a second-order central scheme with nonlinear dissipation developed by [83Jump To The Next Citation Point] to the study of black-hole accretion and formation of relativistic jets, investigating issues such as the magnetic extraction of rotational energy of the black hole [194Jump To The Next Citation Point]. More recently [89Jump To The Next Citation Point] assessed a state-of-the-art third-order, convex, essentially nonoscillatory, central scheme [226] in multidimensional special-relativistic hydrodynamics, later extended to relativistic MHD by [90Jump To The Next Citation Point]. 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 considered by [19Jump To The Next Citation Point20Jump To The Next Citation Point] in one-dimensional special and general-relativistic hydrodynamics and MHD, where results similar to those reported by [89Jump To The Next Citation Point90] are discussed. The scheme of Kurganov–Tadmor (see Equation (83View Equation)) was assessed by [230Jump To The Next Citation Point] in the context of special-relativistic hydrodynamics, using standard numerical experiments such as shock tubes, the shock reflection test, and the relativistic version of the flat-faced step test. As for the other central schemes analyzed by [89Jump To The Next Citation Point19Jump To The Next Citation Point] the results were comparable to those obtained by HRSC schemes based on Riemann solvers, even well inside the ultrarelativistic regime. Lucas–Serrano et al. [230Jump To The Next Citation Point] used high-order reconstruction procedures such as those provided by the PPM scheme [80Jump To The Next Citation Point] and the piecewise hyperbolic method (PHM) [102Jump To The Next Citation Point], which proved essential for keeping the inherent diffusion of central schemes at discontinuities at reasonably low levels. The scheme also produced accurate results in the case of full general-relativistic hydrodynamic simulations involving dynamic spacetimes, as shown by [361Jump To The Next Citation Point] for oscillations of rapidly rotating neutron stars and the merger of neutron-star binaries. The scheme of Kurganov–Tadmor is also the adopted choice in the GRMHD simulations in dynamic spacetimes of [108Jump To The Next Citation Point366Jump To The Next Citation Point]. Finally, a MUSCL-type scheme with HLL numerical fluxes is used in the code of [149Jump To The Next Citation Point], which allows a systematic investigation of GRMHD processes in accretion tori around black holes. Such HLL fluxes are also the choice used in the recent approaches of [396Jump To The Next Citation Point91Jump To The Next Citation Point], which also implement a weighted, essentially nonoscillatory (WENO) method to build up fifth-order convergent numerical codes. Such high-order schemes may be suitable to solve the total (kinetic, thermal, and magnetic) energy equation in GRMHD codes, which deal with the challenging regime posed by high-Mach number flows.

It is worth emphasizing that early pioneer approaches in the field of relativistic hydrodynamics [29173Jump To The Next Citation Point] 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 (W ∼ 2). A key feature lacking in those investigations was the use of a conservative approach for both the system of equations and, accordingly, for the numerical schemes. In light of the results reported in recent investigations employing central schemes [8919Jump To The Next Citation Point230], it appears that this was the ultimate reason preventing the extension of early computations to the ultrarelativistic regime.

The alternative of using high-order component-wise central schemes instead of upwind HRSC schemes becomes apparent when the spectral decomposition of the hyperbolic system under study is unknown. The straightforwardness of a central scheme makes its use very appealing, especially in multiple dimensions where computational efficiency is an issue. Perhaps the most important example in relativistic astrophysics is the system of general-relativistic MHD equations for which, despite some progress, has been achieved in recent years (see, e.g., [34199Jump To The Next Citation Point24Jump To The Next Citation Point23]), much more work is needed concerning their solution with full-wave-decomposition HRSC schemes based upon Riemann solvers. Meanwhile, an obvious choice is the use of central schemes.

4.1.4 Source terms

Most “conservation laws” involve source terms on the right-hand side (RHS) of the equations. In hydrodynamics, for instance, those terms arise when considering external forces such as gravity, which make the RHS 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) and ionization (resulting in a collection of nonhomogeneous 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, directing the interested reader to [398219Jump To The Next Citation Point] and references therein for further details.

There are, essentially, two 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 (76View Equation)):

( ) dUj- = -1-- ˆFj− 1 − ˆFj+ 1 + Snj. (84 ) dt Δx 2 2
The temporal order of this basic algorithm can be improved by introducing successive substeps 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, Un+1 = β„’ Δstβ„’ΔftUn, the operator β„’ Δt f solves the homogeneous PDE in the first step to yield the intermediate value ∗ U. Then, in the second step, the operator Δt β„’s solves the ordinary differential equation ∗ ∗ dU βˆ•dt = S (U ) to yield the final value n+1 U. In order to achieve second-order accuracy (assuming each independent method is second order) a common fractional-step method is the Strang splitting, where Un+1 = β„’ Δstβˆ•2β„’Δftβ„’ Δstβˆ•2Un. This method advances by a half timestep the solution for the ODE containing the source terms, and by a full timestep 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 time scales than the hydrodynamic timestep Δt, or act over much smaller spatial scales than the grid resolution Δx. Stiff source terms may easily lead to numerical difficulties. The interested reader is directed to [219] (and references therein) for further information on various approaches to overcome the problems of stiff source terms.

  Go to previous page Go up Go to next page