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
In hydrodynamic codes based on finite difference or finite volume techniques, Equation (92) is solved on a discrete numerical grid with
Convergence under grid refinement implies that the global error , defined as). Conservation form means that the algorithm can be written as
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 ). 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 , TV(), is defined as: 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 and comparing with Equation (97), one recognizes that the numerical flux function is an approximation to the time-averaged flux across the interface, i.e.,[122, 81] after the seminal work of Godunov , 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 . The book of Toro  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 (that give only first-order accuracy) by better representations of the true flow near the interfaces, let’s say , . The FCT algorithm  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 ). High-order TVD schemes were first constructed by van Leer , who obtained second-order accuracy by using monotonic piecewise linear slopes for cell reconstruction. The piecewise parabolic method (PPM)  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 . 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 , the essentially non-oscillatory (ENO) schemes  and the piecewise-hyperbolic method (PHM) .
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  (see also ). 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 ). 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  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  first solves the Riemann problem, and then distributes the information to the appropriate directions.
© Max Planck Society and the author(s)