Basics of HRSC methods

9.5 Basics of HRSC methods

In this section we introduce the basic notation of finite differencing, and summarize the foundations of HRSC methods for hyperbolic systems of conservation laws. The content of this section is not specific to SRHD, but applies to hydrodynamics in general.

In order to simplify the notation and taking into account that most powerful results have been derived for scalar conservation laws in one spatial dimension, we will restrict ourselves to the initial value problem given by the equation

with the initial condition $u (x,t = 0) = u0(x)$ .

In hydrodynamic codes based on finite difference or finite volume techniques, Equation (92 ) is solved on a discrete numerical grid $(xj,tn)$ with

and

where $Δt$ and $Δx$ are the time step and the zone size, respectively. A difference scheme is a time-marching procedure allowing one to obtain approximations to the solution at the new time, $un+1 j$ , from the approximations in previous time steps. The quantity $n u j$ is an approximation to $n u (xj,t )$ , but, in the case of a conservation law, it is often preferable to view it as an approximation to the average of $u (x, t)$ within a zone $[xj−1∕2,xj+1∕2]$ (i.e., as a zone average), where $xj±1∕2 = (xj + xj±1)∕2$ . Hence

which is consistent with the integral form of the conservation law.

Convergence under grid refinement implies that the global error $||E Δx||$ , defined as

tends to zero as $Δx → 0$ . For hyperbolic systems of conservation laws, methods in conservation form are preferred as they guarantee that if the numerical solution converges, it converges to a weak solution of the original system of equations (Lax–Wendroff theorem [156 ]). Conservation form means that the algorithm can be written as

where $q$ and $r$ are positive integers, and $fˆ$ is a consistent (i.e., $ˆf(u,u, ...,u) = f(u)$ ) numerical flux function.

The Lax–Wendroff theorem cited above does not establish whether the method converges. To guarantee convergence, some form of stability is required, as for linear problems (Lax equivalence theorem [244 ]). In this context the notion of total-variation stability has proven to be very successful, although powerful results have only been obtained for scalar conservation laws. The total variation of a solution at $n t = t$ , TV( $n u$ ), is defined as

A numerical scheme is said to be TV-stable, if TV( $un$ ) is bounded for all $Δt$ at any time for each initial data. One can then prove the following convergence theorem for nonlinear, scalar conservation laws [157

]: For numerical schemes in conservation form with consistent numerical flux functions, TV-stability is a sufficient condition for convergence.

Modern research has focussed on the development of high-order, accurate methods in conservation form, which satisfy the condition of TV-stability. The conservation form is ensured by starting with the integral version of the partial differential equations in conservation form (finite volume methods). Integrating the PDE over a finite spacetime domain $[xj−1∕2,xj+1∕2] × [tn,tn+1]$ and comparing with Equation (97 ), one recognizes that the numerical flux function $ˆ fj+1∕2$ is an approximation to the time-averaged flux across the interface, i.e.,

Note that the flux integral depends on the solution at the zone interface, $u(xj+1∕2,t)$ , during the time step. Hence, a possible procedure is to calculate $u(xj+1∕2,t)$ by solving Riemann problems at every zone interface to obtain

This is the approach followed by an important subset of shock-capturing methods, called Godunov-type methods [122 , 81 ] after the seminal work of Godunov [105 ], who first used an exact Riemann solver in a numerical code. These methods are written in conservation form and use different procedures (Riemann solvers) to compute approximations to $n n u(0;uj,u j+1)$ . The book of Toro [279 ] gives a comprehensive overview of numerical methods based on Riemann solvers. The numerical dissipation required to stabilize an algorithm across discontinuities can also be provided by adding local conservative dissipation terms to standard finite difference methods. This is the approach followed in the symmetric schemes developed in [68 , 249 , 306 , 164 ].

High order of accuracy is usually achieved by using conservative monotonic polynomial functions to interpolate the approximate solution within zones. The idea is to produce more accurate left and right states for the Riemann problem by substituting the mean values $un j$ (that give only first-order accuracy) by better representations of the true flow near the interfaces, let’s say $L uj+1∕2$ , $R u j+1∕2$ . The FCT algorithm [33 ] constitutes an alternative procedure where higher accuracy is obtained by adding an anti-diffusive flux term to the first-order numerical flux. The interpolation algorithms have to preserve the TV-stability of the scheme. This is usually achieved by using monotonic functions which lead to the decrease of the total variation (total-variation-diminishing schemes, TVD [120 ]). High-order TVD schemes were first constructed by van Leer [282 ], who obtained second-order accuracy by using monotonic piecewise linear slopes for cell reconstruction. The piecewise parabolic method (PPM) [60 ] provides even higher accuracy. The TVD property implies TV-stability, but can be too restrictive. In fact, TVD methods degenerate to first-order accuracy at extreme points [216 ]. Hence, other reconstruction alternatives have been developed where some growth of the total variation is allowed. This is the case for the total-variation-bounded (TVB) schemes [259 ], the essentially non-oscillatory (ENO) schemes [121 ] and the piecewise-hyperbolic method (PHM) [175 ].

There are several strategies to extend HRSC methods to more than one spatial dimension. A brief summary of these strategies can be found in LeVeque’s book [157 ] (see also [160 ]). The simplest strategy is dimensional splitting, where the differential operators along the spatial directions are applied in successive steps (fractional step methods). Second order in time is achieved when one permutes cyclically the order in which the directional (i.e., 1D) operators are applied (Strang splitting [270 ]). In semi-discrete methods (method of lines), the process of discretization proceeds in two stages. First only operators involving spatial derivatives are discretized, leaving the problem continuous in time. This gives rise to a system of ordinary differential equations (in time) which can be integrated by any ODE solver. In the method of lines approach, the numerical fluxes across cell interfaces are computed in all two or three spatial directions, before they are simultaneously applied to advance the equations. Particularly of interest are TVD Runge–Kutta time discretization algorithms [260 , 261 ], which preserve the TVD properties of the algorithm at every substep. A third approach relies on unsplit methods, where the different spatial directions are also advanced simultaneously as in the semi-discrete methods. However, the extension of unsplit methods to second-order accuracy requires incorporating not only slopes in the normal direction (as in one-dimensional or split algorithms), but also cross-derivatives arising from the multi-dimensional Taylor series expansion. Good examples of genuinely multi-dimensional upwind methods for hyperbolic conservation laws (using slightly different strategies) are those described in [58 , 158 ]. In [58 ] the algorithm proceeds in two steps. First, interface values are interpolated, using information from all orthogonal directions. Secondly, the Riemann problems defined by these interface values are solved. The algorithm proposed in [158 ] first solves the Riemann problem, and then distributes the information to the appropriate directions.