A residual-based artiﬁcial viscosity ﬁnite difference method for scalar conservation laws

In this paper, we present an accurate, stable and robust shock-capturing ﬁnite difference method for solving scalar non-linear conservation laws. The spatial discretization uses high-order accurate upwind summation-by-parts ﬁnite difference operators combined with weakly imposed boundary conditions via simultaneous-approximation-terms. The method is an extension of the residual-based artiﬁcial viscosity methods developed in the ﬁnite- and spectral element communities to the ﬁnite difference setting. The three main ingredients of the proposed method are: ( i ) shock detection provided by a residual-based error estimator; ( ii ) ﬁrst-order viscosity applied in regions with strong discontinuities; ( iii ) additional dampening of spurious oscillations provided by high-order dissipation from the upwind ﬁnite difference operators. The method is shown to be stable for skew-symmetric discretizations of the advective ﬂux. Accuracy and robustness are shown by solving several benchmark problems in 2D for convex and non-convex ﬂuxes.


Introduction
Conservation laws are used to describe flow processes in many different areas of physics, e.g., in fluid dynamics, electromagnetics and geophysics, and are therefore also of relevance in many engineering applications. Deriving numerical methods for the treatment of non-linear conservation laws that are accurate and robust is a vibrant field of research that has received considerable attention in the last decades. Already in the 1950s, Lax [21] showed that first-order convergence is readily achieved by using the upwind flux. Later, high-resolution methods based on modifications of central upwind schemes such as MUSCL, ENO and WENO were developed and successfully used in solving conservation laws, see e.g., [46,16,23,47,3].
Despite their success in solving many complex problems, high-resolution methods based on central upwind schemes have been shown to fail when solving conservation laws with non-convex fluxes, see Kurganov et al. [19]. Therefore, the development of stable high-order schemes that are convergent to the physically realizable entropy solution has been, and still is, an active research topic. One technique used to ensure that a method is entropy stable is to construct the scheme such that it satisfies a discrete entropy inequality, mimicking the underlying continuous entropy inequality, see e.g., Fisher and Carpenter [7] for an entropy stable high-order WENO finite difference method (FDM).
Viscous regularization is an alternative numerical technique used for solving conservation laws that was initially introduced by Richtmayer and Von Neumann in the 1950s. The idea consists of adding a mesh-dependent parabolic operator, the so-called artificial viscosity, to the equation. In the original artificial viscosity method of Richtmayer and Von Neumann, the viscosity was constructed to be proportional to the gradient of the solution, which is inconsistent with the differential equation. Moreover, the viscosity loses regularity close to shocks and sharp discontinuities. To overcome this issue Szepessy [43] proposed constructing the artificial viscosity term proportional to the residual of the differential equation. The resulting term is then consistent with the partial differential equation (PDE), since the residual is zero for an exact solution. Moreover, Szepessy proved that by augmenting the residual-based artificial viscosity term to the well-known streamline diffusion method one obtains a scheme for which the solution is convergent to the entropy solution, see e.g., [43,17].
The direct extension of Szepessy's results to the finite difference counterpart was not possible due to the least-squares construction of the streamline diffusion term. Recently, Nazarov [32] proved that by disregarding least-squares entirely and using only the residual-based artificial viscosity term in the standard Galerkin finite element approximation one can prove that the approximation is uniformly bounded in the L ∞ -norm, weakly consistent with all entropy inequalities and strongly consistent with the initial data, and therefore converges to the unique entropy solution. From here on we refer to this method as the RV method, i.e., the method is stabilized using only a residual-based artificial viscosity term. Extensions of the RV method to systems of conservation laws, in the framework of continuous and discontinuous, high-order Lagrange, and spectral element methods, are presented in [33,26,24].
The RV method is closely related to the entropy viscosity method (see e.g., [10,12,11,34,45]), where the artificial viscosity coefficient is proportional to the entropy inequality residual. There are many entropies for scalar conservation laws and in particular, the PDE solution can be one of the entropies. For scalar conservation laws RV is therefore a specific case of entropy viscosity. However, this is not the case for systems of conservation laws. It should also be noted that although entropy viscosity has shown to be a highly accurate and robust method, a theoretical convergence proof, equivalent to that presented in [32] for RV, is not yet known.
The focus of the present work is on presenting a new shock-capturing method in the framework of summation-by-parts (SBP) finite differences, paired with the simultaneous-approximation-term (SAT) technique for the imposition of boundary conditions. For the derivation and detailed presentations of the SBP-SAT methodology, we refer the reader to [18,40,29,30,4,6,42]. The main contribution of the present work is extending the RV method to the finite difference SBP-SAT framework, focusing on scalar conservation laws. Since the usual projection techniques used to compute the residual in finite element or spectral collocation methods do not carry over to FDMs, this extension is not straightforward. We present novel techniques for computing the residual in the finite difference framework. As pointed out in Guermond et al. [12,Sec. 2.4] and Upperman and Yamaleev [45,Sec. 10.2], the residual has a highly oscillatory behavior, which can result in small scale oscillations in smooth regions of the numerical solution. To overcome this issue, we combine the upwind SBP finite difference operators presented in Mattsson [28] with the proposed RV method. The upwind SBP operators naturally provide high-order accurate dissipation, suppressing high-frequency oscillations when paired with standard flux-splitting methods, see e.g., [28,39,25,22]. Our numerical tests show that the proposed method captures shocks accurately, is high-order accurate, does not invoke any small scale oscillations in smooth regions, and can be applied to both convex and non-convex fluxes.
The structure of the paper is as follows: In Section 2 the SBP FDM is introduced, together with relevant definitions. We then continue to present the RV method in the SBP-SAT finite difference setting in Section 3. In Section 4 the method is analyzed in 1D, where energy and stability estimates are developed, and conservation properties of the method are analyzed. The method is then applied to Burgers' equation in 1D demonstrating how boundary conditions are imposed weakly via SAT, in a stable and consistent fashion. The section is ended by showing that a bound on the CFL condition is obtained under certain assumptions on the discrete advective flux. The theoretical results are corroborated by computations for a nonsmooth problem presented in Section 5. Next, the method is extended to 2D in Section 6, and it is demonstrated how to weakly impose boundary conditions leading to linear stability estimates, for problems where no skew-symmetric split form of the advective flux is available. In Section 7 the method is tested against a set of benchmark problems, which include both smooth and non-smooth cases with convex and non-convex fluxes, showing that the method is shock-capturing while still preserving high-order accuracy outside of shock locations. Finally, the paper is concluded and future developments are discussed in Section 8.

The finite difference method
Traditionally, SBP finite difference operators are central difference stencils constructed with careful boundary closures, such that the operators yield a SBP property with respect to a discrete norm, mimicking the underlying integration-by-parts (IBP) property of the L 2 -norm. For the derivation of traditional SBP operators see e.g., [18,40,29,27]. The FDM in the present work is based on diagonal-norm SBP operators. Since only diagonal-norm operators are utilized the term 'diagonal-norm' will be omitted from here on, definitions excluded.

Notation
Before introducing the formal definitions of SBP operators, we first establish some notation. To start with, the discretization of the computational domain and related concepts are defined. The definitions are given in 1D. For the extension to 2D see Section 6.1.1. Consider the domain = [x l , x r ]. The domain is discretized using m equidistant grid points as For a scalar function v(x, t) ∈ L 2 ( ) × [0, T ], let the grid function v be its restriction onto the grid, i.e., v = . Let ·, · denote the L 2 inner product on and · the corresponding norm. Similarly, in the discrete setting, let H be a symmetric positive definite matrix, defining a quadrature rule. Then the discrete inner product ·, · H is defined as v, w H = v T Hw for grid functions v, w. A discrete norm is then defined as w H = w, w H . The grid function u will solely be used as the approximate solution to an initial-boundary value problem (IBVP) with unknown u.
For a scalar conservation law in 1D the continuous flux function is denoted f (u), while the discrete counterpart is denoted f (u), i.e. f (u) is the grid function restriction of f (u). The notation f (u) and f (u) will be used to denote the continuous and discrete flux Jacobian. When the dimension of the problem is greater than 1, the continuous flux function is vector valued which is indicated using boldface notation f(u), with the discrete counterpart being denoted f(u). Similarly f (u) and f (u) are used to denote the continuous and discrete flux Jacobian in higher dimensions. The discrete advective flux, approximating ∇ · f(u) is denoted D f (u). Note that D f (u) could also be based on some other form of the flux (e.g., quasi-linear, or split forms), as discussed in Section 4, and Section 6.
Finally, the following boundary value operators will be used in the following definitions and analysis where e l and e r are m element vectors. For a grid function v, the left and right boundary elements are then given by v l = e T l v, and v r = e T r v respectively.

Central difference SBP operators
First-derivative SBP finite difference operators are introduced by the following definitions given in Mattsson [28]: In particular, for v = w = u (2) reduces to u, Second-derivative SBP finite difference operators are introduced by the following definition, also presented in Mattsson [28]:  [18,40,27]. In the present work we will primarily make use of wide-stencil second-derivative operators. Let A be the diagonal matrix obtained by projecting a smooth function a ≥ 0 onto the diagonal. Wide-stencil second-derivative SBP operators can then be constructed by applying a first-derivative SBP operator twice, e.g., where √ A is the diagonal matrix with elements √ a ii , i = 1, . . . , m. In particular, for v = w = u (4) reduces to

Upwind SBP operators
As previously mentioned, the present work utilizes the upwind SBP operators derived in Mattsson [28]. The upwind SBP operators are SBP finite difference operators with a non-central interior stencil, which naturally introduce dissipation for high-frequency modes in such a way that linear stability and accuracy of the approximation is preserved without introducing unnecessary stiffness. Recent examples of realistic settings where the upwind SBP operators have been successfully used include [39,25,22]. Narrow-stencil upwind SBP first-derivative operators are introduced through the following definition: accurate narrow-stencils, are said to be pth-order narrow-diagonal upwind SBP operators if the diagonal matrix H defines a discrete norm, Q + + Q T − = 0, and A few useful relations involving the first-derivative upwind SBP operators are given by By combining relations (6) with flux splitting techniques it is possible to rewrite discretizations using first-derivative upwind operators D + and D − as discretizations using central differences D 1 , with additional terms proportional to H −1 S, where S is negative semi-definite by Definition 2.4. Hence the upwind SBP operators introduce a small dissipative correction to the truncation error, as compared to a corresponding central difference approximation using D 1 . Upwind SBP second-derivative operators are given by the following definition: , are said to be pthorder diagonal-norm second-derivative upwind SBP operators. Here D + and D − are pth-order diagonal-norm upwind SBP operators and A is the diagonal matrix obtained by projecting a smooth function a ≥ 0 onto the diagonal.
The second-derivative upwind operators satisfy the relations D (a) and it therefore follows that D 2± are wide-stencil second-derivative operators, satisfying the SBP properties (4)-(5), replacing D 1 with D ± , respectively. The SBP properties are restated below for sake of completeness.

v, D (a)
where d T 1r ± , d T 1l ± are the first and last rows of D ± respectively.

Framework
For a conserved quantity u with the flux function f(u) ∈ C 1 (R; R d ) defined on the d-dimensional domain , consider the associated scalar conservation law IBVP given by with boundary conditions Lu = g on an inflow boundary , such the resulting problem is well-posed. Discretizing (10) in space using SBP finite difference operators, results in the semi-discrete problem where D f is an SBP finite difference discretization of the advective flux and (ε) is an SBP finite difference discretization of the variable coefficient Laplacian operator, introducing viscous regularization proportional to the RV ε. Boundary conditions are imposed weakly through the SAT, in a stable and consistent fashion. In the present work (11) is discretized with the upwind SBP operators introduced in Section 2.3 in combination with flux-splitting techniques, such that high-order artificial dissipation is introduced. However, the high-order dissipation by itself cannot handle strong discontinuities and shocks, see Section 7.3. For this reason our aim is to extend the RV method presented in Nazarov [32] to the finite difference framework. The RV is constructed as where ε 1 is the first-order viscosity vector and ε r the high-order residual viscosity vector. The purpose of ε 1 is to effectively provide an upper bound for the viscosity, such that first-order convergence is obtained at the shock. The shock location is detected by the error estimator built into ε r . Let i be a d-dimensional multi-index. At grid point x i , ε 1 is computed as where h i is a weighted average of the grid spacings connected to x i , c 1 is an order unity constant and w i is the local wave speed, which should be close to the true wave speed |f (u(x i ))|. Typically, one simply constructs the local wave speed as w i = |f (u) i |. One argument for constructing ε 1 in this way is that when discretizing (11) using second-order central differences and choosing c 1 = 1 2 , one retains the traditional first-order upwind scheme for ε = ε 1 . For problems where f could be zero at a point in the shock, the local wave speed is instead constructed as where the notation max loc is used to denote the maximum taken over the grid points connecting to x i . We refer to (14) as the modified local wave speed. The high-order residual viscosity at grid point x i is computed as where c r is another order unity constant (typically chosen as 1), R(u) is the discrete residual and n(u) i is a normalizing function at grid point x i . The discrete residual R(u) approximates the true residual error R(u) = u t + ∇ · f(u) and is given by where D t (u) and D f (u) are discretizations of u t and ∇ · f(u), respectively. Note that D f (u) may differ from D f (u) in (11).
The discretization of R(u) will be discussed in detail in Section 3.2. The normalization n(u) i is chosen such that the residual is amplified in regions with jumps in u and suppressed elsewhere. It is given by where min loc is defined analogously to max loc . The application of max loc in (15) can be viewed as a post-processing of the normalized residual, which increases the area in which the residual viscosity is active, thus providing additional stabilization.

Remark 3.1.
In the present work, all numerical simulations are performed using c 1 = 1 2 and c r = 1. In addition, we will only consider equidistantly spaced grids, and the weighted grid spacing h i is therefore chosen as |h| where h is the vector holding the grid spacings for each coordinate direction.

Remark 3.2.
Although the RV framework is presented in the context of SBP finite differences, the key elements for proving stability of (11) are the combination of SBP operators and stable boundary treatment. One might therefore consider using other types of SBP operators, for instance the generalized SBP operators presented in e.g., [37,5], or other methods for imposing boundary conditions, such as the projection method presented in Mattsson and Olsson [30].

Discretizing the residual
The key component of the presented method lies in constructing an accurate approximation of the residual error R(u) = u t + ∇ · f(u). We propose two methods of constructing the discrete residual R(u) in (16). For both methods D f (u) is straightforwardly constructed using upwind SBP difference operators, but differ in the construction of the operator D t (u).
The first method involves storing the solution vector over multiple time levels and computing D t (u) using standard highorder backward differentiation formulas (BDFs). The resulting residual is thus constructed as R bdf (u) = BDF(u) +D f (u). (18) BDFs of orders 1-6 can be found in many standard text books on numerics, e.g., Hairer et al. [15,III.1]. In addition, BDFs of orders higher than 6 are listed in Appendix A. The corresponding RV method using R bdf is referred to as RV bdf .
For the second method, note that by discretizing the conservation law (10) in space without adding RV one obtains which will be non-zero if D f (u) and D f (u) are constructed using different operators. This approach is similar to the ones introduced recently by Lu et al. [24] in the spectral element framework, where it was termed advection-based residual estimator, and Upperman and Yamaleev [45] in the SBP spectral collocation framework. More precisely, in [24] different polynomial spaces on a fixed mesh were used to approximate D f (u) and D f (u), whereas in [45] the time derivative in the discrete entropy residual is substituted using the discrete entropy formulation of the PDE. However, in the current work, we propose a slightly different approach. In order to incorporate errors originating both from spatial resolution (i.e., how well the shock is resolved) and operator accuracy into the residual error, D f (u) is computed on a coarser grid using a lower-order accurate operator compared to D f (u). In particular, the discrepancy between D f (u) and D f (u) will be more prominent where u changes rapidly, and smaller where u is smooth. Thus, for the second method the residual is constructed as where D (c) f is the coarse grid advective flux operator, and I C 2F , I F 2C are interpolation operators, where I C 2F interpolates from a coarse grid to the fine solution grid, and vice versa for I F 2C . The RV method using R mg is referred to as RV mg . The subscript mg is shorthand for multi-grid, referring to the fact that the residual is computed on both a coarse and fine grid.

Remark 3.3.
Constructing D f using upwind SBP difference operators will add dissipation targeting high-frequency oscillations of the residual. This is effectively similar to the smoothening of residual oscillations utilized in Guermond et al. [12].

Remark 3.4.
Note that there are no stability requirements to account for when discretizing the residual. One can therefore construct the BDFs and interpolation operators in (18) and (19) of arbitrary accuracy. Of course, there are practical limitations to this in terms of computational cost and memory usage.
Remark 3.5. The interpolation operators in the present work utilize standard high-order polynomial interpolation, with central stencils in the interior and one-sided stencils at the boundaries. In addition, the coarse and fine grids are chosen to be aligned, with a coarse-to-fine ratio of 1:2. Other types of interpolation operators and grid ratios are of course possible, but it is not something that has been considered in the present work.

Computational considerations
In this section we present a short discussion on some practical aspects of using the RV method. Due to the scaling with h 2 i in (15), the formal accuracy of the discrete residual R(u) can be two orders lower compared to the accuracy of the scheme (11), without affecting the order of accuracy of the solution in smooth regions. This opens up for some flexibility in what accuracy to use when discretizing the residual. Unfortunately, it is not clear if there is an optimal choice of the order of accuracy for the operators. To determine what orders to use, the benchmark problems in Section 7 were studied numerically by considering the accuracy and over-and undershoot in the solutions. The qualitative results are summarized below.
To start with, it should be noted that for a pth order accurate scheme (11), the RV method converges with the expected convergence rate in the l 1 -and l 2 -norms as long as the residual is discretized with at least order of accuracy p − 2. The choice of orders does however have an effect on the amount of over-and undershoot present in solutions. For RV bdf we found that choosing the order of accuracy of D f as p − 2, together with a 6th order accurate BDF provided sufficient shock-capturing while maintaining accuracy in smooth regions. This setup has been used for the computations using RV bdf performed in the present work, excluding those in Section 7.3.1. In Section 7.3.1 it is also shown that increasing the order of the BDF has a positive effect on reducing over-and undershoot in solutions containing shocks. Similarly for RV mg , choosing the order of accuracy of D f as p + 1, and of D (c) f and I C 2F as p − 1 was found to provide overall good results. It is also possible to reuse the computed advective flux (i.e., using D f = D f ), but it was found that this approach provided less effective shock-capturing for the problems studied herein.
Note also that R mg in (19) can be constructed for the current time level without requiring information from previous time levels, whereas a pth order BDF requires the solution vector from p + 1 time levels to be stored. This simplifies the use of RV mg in adaptive time-stepping methods. In addition, it is also clear that R mg may be less memory intensive to compute compared to R bdf since the spatial operators can be implemented using a matrix-free approach. However, the additional application of interpolation operators and the coarse grid discrete advective flux may render R mg more computationally costly to evaluate compared to R bdf . The latter is demonstrated at the end of Section 5.

Time integration algorithm
We now go into detail on how the RV is integrated in explicit time-stepping schemes. For simplicity, the temporal domain [0, T ] is discretized using n equidistantly spaced points, with time step size k. In each time step the RV at the current time level is computed, after which the solution vector is updated using a suitable time-integrator and (11). Due to the scaling by h i in (15), the time step may be chosen linearly, i.e., k ∼ h ≡ min i∈1,...,m (h i ), see Section 4.2. For multi-stage methods, the viscosity is held fixed in all substages. Note that RV bdf requires that the solution vectors from previous time levels are stored in order to compute the residual. Denote a BDF of pth order accuracy by BDF (p) . (p) requires the solution vectors from p + 1 time levels. In order to apply the residual viscosity at the first possible time-step, higher-order BDFs may be used as solutions from previous time levels are accumulated. Thus if r and s denote the accuracy of the lowest and highest order accurate BDFs used, the time integration procedure for RV bdf is as follows: Advance to u i+1 using u i and ε 0 end /* Advance in time with increasing BDF orders Compute ε i by RV bdf using u i and BDF (s) (u i ) Advance to u i+1 using u i and ε i end Remark 3.6. Using r < 3 was found to result in a reduced order of accuracy for the smooth problem studied in Section 7.1. In the computations using RV bdf , r = 3 is therefore used.
The time-stepping procedure for RV mg is more straightforward, since it does not require the solution vector to be stored over multiple time levels. It is given by: Algorithm 2 Time-stepping procedure using RV mg .
Compute ε i by RV mg using u i Advance to u i+1 using u i and ε i end As mentioned above, the presented time integration algorithms fix the viscosity prior to updating the solution. When using RV mg and a multi-stage time integrator, such as an explicit Runge-Kutta method, it is possible to compute the RV to design order accuracy in each substage. However, as illustrated by the numerical results in Sections 5 and 7, the RV method is still accurate and robust even if the viscosity is fixed in each substage. Fixing the viscosity is also beneficial in that it reduces the amount of computations needed for multi-stage methods.

Analysis in 1D
In order to clearly demonstrate how the RV method fit into the SBP framework, the problem is first studied in 1D on a Cartesian domain. Consider the 1D scalar conservation law IBVP (10) on the domain = [x l , x r ] given by where now D f (u) is the advective flux. Assuming sufficient smoothness, the advective flux can be expressed in many equivalent forms, e.g., conservative D f (u) = f x (u), quasi-linear D f (u) = f (u)u x and split forms based on the product-or chain rule. For the purpose of the analysis in this section, we will restrict ourselves to cases where the form of D f (u) is such that it admits an energy estimate, which is achieved when D f (u) satisfies the skew-symmetric IBP property Here B T f are boundary terms containing boundary fluxes whose exact form depends on D f (u). For example, such a skewsymmetric form can be obtained for Burgers' equation by splitting the flux by its conservative and quasi-linear parts (see Section 4.1.3). If (21) is satisfied, the non-linear energy estimate is straightforwardly obtained by taking the inner product u, u t , using (20) and (21) (referred to as the energy method), resulting in It is sufficient to consider the case when Lu(x, t) = 0 on inflow boundaries, leading to d dt showing energy stability. In cases where a skew-symmetric form of D f (u) is not explicitly known, linear energy estimates may still be obtained using the method of frozen coefficients. Although such energy estimates are only valid under the assumption of sufficient smoothness of the problem, they might still be practically useful in the development of stable boundary treatment. In Section 6 we present such an analysis.
Remark 4.1. In Tadmor [44] it is shown that a skew-symmetric form of D f (u) exists for the corresponding Cauchy problem to (20) iff an entropy function for the problem exists. In the skew-symmetric form entropy conservation is established straightforwardly using integration by parts. One can therefore view the above analysis to be restricted to problems for which the energy serves as a generalized entropy function. We do however emphasize that the focus of the present work is on viscous regularization and not on energy-or entropy conservative discretizations of the advective flux, although the latter is useful for the purpose of the analysis.

Semi-discrete setting
Now consider the semi-discrete setting, in which time is left continuous while space is discretized. As previously stated, the present work utilizes a combination of high-order upwind dissipation and low-order RV-based diffusion to treat oscillations located both at and away from shocks. To discretize the advective flux using upwind operators, a flux-splitting technique is required. We therefore start by discussing flux-splitting of a general scalar conservation law, before presenting the semi-discrete SBP-SAT discretization and analyze its stability in Section 4.1.2. By the famous theorem of Lax and Wendroff [20], a scheme must preserve conservation of fluxes on the discrete level, in order for the numerical solution to converge to a weak solution of the conservation law. Therefore, a brief analysis of the conservation properties of the discretization is included in addition to the stability analysis. Finally, the analysis is exemplified in Section 4.1.3, by applying it to Burgers' equation in 1D, for which it is also demonstrated how a stable boundary treatment using SAT is derived.

Flux-splitting
Consider the quasi-linear form of a scalar conservation law, Let a(u) = − f (u) and split the flux term into left-and right-going characteristics as u t = a + u x + a − u x , where 2a ± = a ± r for some r ≥ 0 such that a + ≥ 0 and a − ≤ 0. In the present work, a global Lax-Friedrichs type splitting is employed, given by choosing r = a ∞ = u ∞ . In the discrete setting, let F be the matrix obtained by projecting f (u) onto the diagonal and let An upwind discretization of the conservation law in quasi-linear form is then given by Using the relations in (6), the above discretization can be shown to be equivalent to i.e., the scalar conservation law in quasi-linear form discretized using a central difference SBP operator D 1 , but with an added term f (u) ∞ H −1 Su introducing high-order artificial dissipation to the scheme. For further details regarding upwind operators and flux-splitting, see Mattsson [28]. In the present work, a high-order upwind dissipative term f (u) ∞ H −1 Su is applied regardless of the underlying form of the advective flux any form, based on the central difference SBP operator D 1 , the corresponding upwind discretization is formulated as It should be noted that many types of flux splitting techniques are available, and choosing an optimal type of splitting has not been the focus of the present study. However, the numerical studies in Section 5 and Section 7, illustrate the global Lax-Friedrichs type splitting has shown to be effective for the problems studied in the present work.

SBP-SAT discretization -stability & conservation analysis
Let the domain be discretized as in (1). An upwind SBP-SAT discretization of the IBVP (20) regularized with RV is given by where D f (u) is given by (24), D (ε) given by Definition 2.5, and E is the matrix obtained by projecting of the RV vector ε in (12) onto the diagonal. The SAT imposes the boundary conditions Lu = g weakly.
2− is arbitrary. Both are second-derivative upwind operators, and are interchangeable.
In order to obtain a stability estimate, the central difference part D * (24) is assumed to discretize an advective flux in skew-symmetric form, satisfying the SBP property u, D * where S is negative semi-definite, given by Definition 2.4, adding high-order accurate dissipation. We start by analyzing the stability of (25). Stability is shown by obtaining a semi-discrete energy estimate, by applying the energy method. In the semi-discrete setting the energy method consists of taking the inner product u, u t H using (25) and the SBP properties (26) and (9), resulting in For stability it is required that the right-hand-side of (27) is non-positive, assuming boundary data g = 0 is imposed on inflow boundaries. Since S ≤ 0, stability follows if B T ≤ 0. Therefore, the SAT should treat the energy contribution on the boundary arising from u, D 2+ u H , in addition to imposing boundary conditions on inflow boundaries. Thus the above estimate is useful in guiding the construction of stable boundary treatment. Given B T ≤ 0, integrating (27) in time results in the following semi-discrete energy estimate and thus the scheme (25) is stable. Note that (28) is a semi-discrete counterpart of (23), with additional control over the discrete derivative D + u in regions where ε is large, as well as additional upwind dissipation for high-frequency modes.
Next, we analyze the conservation properties of (25). To start with, the addition of the high-order upwind dissipation does not affect the conservation properties of the discrete advective flux D f (u) in (24). To see this, let 1 and 0 respectively be the constant vectors with the values 1, and 0 everywhere. The conservation relation for the discrete advective flux is then given by the inner product 1, D f (u) H . By (24) we have where we used that S is symmetric by definition, and that S1 = 0, following from S consisting of higher-order derivatives that are exact for constant functions (see Mattsson [28] for the derivation of the upwind SBP operators). The conservation properties of D f (u) therefore only depend on the conservation properties of the central difference discretization D * f (u), which is constructed using the standard diagonal-norm D 1 operator given by Definition 2.2. Furthermore, it is shown in Fisher et al. [8] that all split forms of the discrete advective flux that are based on the conservative and product rule forms, and discretized using diagonal-norm skew-symmetric SBP operators are discretely conservative (regardless of whether they admit an energy estimate or not). Thus by the results of Fisher et al. [8], the upwind SBP discretizations of the advective flux D f (u) presented herein are discretely conservative. The SBP conservation relation for the viscous term is established similarly, by the use of (8) where again the fact that D + 1 = 0 is exact for constant functions was used. Thus both the advective and viscous term preserve the conservation property on a discrete level.

Example: Burgers equation
Burgers' equation is the classical scalar example for which a skew-symmetric form of the advective flux is known, and therefore is well-suited to demonstrate how to formulate boundary conditions and a corresponding SAT such that the RV regularized scheme is provably energy stable. The skew-symmetric form of D f (u) is obtained by splitting the flux by one conservative part f x (u) = (u 2 /2) x and one quasi-linear part f (u)u x = uu x via the chain (or product) rule, given by It is easy to verify that D f (u) satisfies the IBP relation It is then clear that imposing Dirichlet boundary conditions on inflow boundaries i.e., results in an energy stable problem. The corresponding SBP upwind discretization (24) of (29) is given by and similarly D f (u) satisfies the SBP relation The estimate (27) is therefore obtained where B T now is given by The SAT should therefore be formulated such that it imposes the boundary conditions (30), and that B T ≤ 0 when considering g l , g r = 0 in order for the scheme to be accurate and stable. The following ansatz is proposed for imposing Dirichlet conditions on inflow boundaries.
where the subscripts l and r denote the left and right boundary respectively. The parameters γ l , γ r , σ l and σ r are penalty parameters which must be tuned for stability. Stability of the scheme is given by the following lemma: Lemma 4.1. The scheme (25) with the advective flux (31) and SATs (33) applied to inflow boundaries is stable for penalty parameters γ l ≤ −1, γ r ≤ −1, σ l = σ r = 1.

Proof.
To determine conditions for stability, by (27) it is enough to consider conditions for when B T ≤ 0, assuming homogeneous boundary data g l = g r = 0. We treat both SAT l and SAT r simultaneously by assuming that there is inflow from both the left and right boundary (i.e., u l ≥ 0, u r ≤ 0), and setting SAT = SAT D l + SAT D r . By (32) we have Since u l ≥ 0, u r ≤ 0, the first two terms are non-positive for γ l ≤ −1, γ r ≤ −1. Furthermore, the last two terms vanish for σ l = σ r = 1. Thus for γ l ≤ −1, γ r ≤ −1, σ l = σ r = 1 we have B T ≤ 0, which completes the proof.
On outflow boundaries, no physical boundary condition should be imposed. However, the SATs must still treat the boundary terms originating from u, D (ε) 2+ u H , and therefore the following ansatz for artificial Neumann boundary conditions is proposed for outflow boundaries The above terms impose the boundary conditions εu x = 0 on the left or right boundary, respectively. The condition is trivially satisfied if ε is zero on the boundary, while if ε = 0 on the boundary, the condition u x = 0 will be enforced. Note however that since lim h→0 ε = 0 by (12), (13) and (15), the boundary term will go to zero in grid refinement. (25) with the advective flux (31) and SATs (34) applied to outflow boundaries is stable for penalty parameters σ l = σ r = 1.

Lemma 4.2. The scheme
Proof. The proof is similar to that of Lemma 4.1. To determine stability conditions for σ l and σ r simultaneously we assume that both the left and right boundaries are outflow boundaries (i.e., u l ≤ 0, u r ≥ 0) and consider (32) The boundary terms vanish for σ l = σ r = 1, such that B T ≤ 0 and stability follows by (27).

A bound on the CFL condition
The aim of the following analysis is to show that the CFL condition is not severely affected by the addition of RV, which is essentially due to the fact that the artificial viscosity ε in (12) is scaled by the grid spacing h j = h for every j = 1, . . . , m (in 1D). This is done by obtaining a coarse, fully discrete energy estimate where, to simplify the analysis, we neglect the treatment of boundaries by considering periodic solutions to the Cauchy problem, and discretize in space using standard central difference SBP operators. The estimate is obtained under the following two assumptions. The first assumption is that the discrete advective flux D f (u) is skew-symmetric, which when discretizing using periodic central difference operators yields the SBP property u, D f (u) H = 0. The second assumption is that D f (u) satisfies the following bound relating to the discrete chain rule, for some constant C independent of h, where W i is the diagonal matrix obtained by projecting the local wave speed w i j , j = 1, . . . , m in (13) onto the diagonal, and the superscript i denoting the time level. Note that if we use the unmodified local wave speed, then w i j = | f (u i ) j |, which clarifies the connection to the discrete chain rule.
Consider the Cauchy problem Discretizing in space using periodic standard central difference SBP operators, regularizing with RV, and discretizing in time using n time levels, results in the following fully discrete initial-value-problem where D τ u i+1 = u i+1 −u i τ is the forward Euler difference operator and τ is the discrete time step. As already stated, the aim of the analysis is to find a bound on the CFL condition. Since the stability region of forward Euler is contained in explicit Runge-Kutta methods the results are relevant also for higher-order Runge-Kutta methods. For the purpose of the analysis we must make the assumption that ε > 0, since otherwise the spatial discretization allows for purely imaginary eigenvalues of the resulting system of equations, for which the forward Euler method is unstable. Practically, ε ≥ 0 is of course still allowed when employing a Runge-Kutta method whose stability region covers purely imaginary eigenvalues. (35). Then, there exists a constant κ independent of h, such that a solution to (36) satisfies the discrete energy estimate

Proposition 4.3. Assume that D f (u) is skew-symmetric and satisfies
Proof. To obtain a discrete energy estimate, the discrete energy method is applied to (36), which consists of taking the since D f (u) is skew-symmetric by assumption and the problem is periodic. Next, by completing the square, the term which, when inserted into (38) yields To bound the right-hand-side term τ 2 D τ u i+1 2 H , observe that by Hölder's inequality. Next, note that by (12) and (13) it follows that holds. Note that (41) also holds when the modified local wave speed in (14) is used. We start by obtaining a bound of the viscous term D where (41) was used, together with the inverse inequality D 1 v H ≤Ch −1 v H =Ch −1 . It then follows that which when combined with (40) and (39) results in What's left to do is show that the term τ 2 D f (u i ) 2 H remains bounded. We show that this is true when the first-order viscosity is active. For this case, the artificial viscosity at grid point x j satisfies ε i j = c 1 hw i j , j = 1, 2, . . . ,m.
Note that the estimate (37) is a fully discrete analog to (28) (neglecting the high-order upwind dissipative term), but is coarse in the sense that it assumes that the first-order viscosity is active everywhere. A sharper estimate would involve a bound of the solution in non-viscous regions (i.e., where ε = ε r ) as well. Obtaining such an estimate could potentially be a first step towards proving convergence of RV FDMs. Extending the theorem to include non-viscous regions is currently ongoing work. The coarse estimate (37) is however instructive since it gives a bound on the CFL condition, showing that the time step restriction is given by k ∼ h. The numerical studies conducted in this paper corroborate the result, and show that a linear time step may be used also for problems on bounded domains where the discrete advective flux D f (u) necessarily is not skew-symmetric (see Section 5 and Section 7).

Remark 4.4.
The question arises when the assumption of the estimate (35) is valid. How to construct or verify that a discrete chain rule is satisfied is not trivial and is ongoing research (see e.g., Ranocha [38]). The estimate can be motivated for the discrete advective flux in conservative form, using second order accurate central difference operators, by considering the mean value theorem on the interval Dividing by 2h, and taking the absolute value, one obtains |(D 1 f (u)) j | = | f (ξ j )(D 1 u) j | (from this expression one can easily form the norm · H ). We now consider three cases: In this case we can always construct a more robust first-order viscosity which is evaluated at some point ξ * j ∈ I j such that f (ξ * j ) = 0. Then we again have =⇒ |(D 1 f (u)) j | ≤ C j |w j (D 1 u) j |. Thus the proof of Proposition 4.3 can be repeated using this modified first-order viscosity. The numerical studies carried out in the present work however illustrate that a non-restrictive CFL can be used for problems where f may be zero, using the modified local wave speed in (14). Finally, if the inequality holds for the conservative form of D f (u) then it also holds for split-forms of D f (u) based on the conservative and quasi-linear form (such as the skew-symmetric form used in Section 4.1.3).

Computations in 1D
Consider Burgers' equation in conservative form, given by (20) with D f (u) = (u 2 /2) x on the domain = [0, 1] and boundary condition u(0, t) = 1. For the initial condition u 0 (x) = 1, the analytic solution is given by The problem is discretized in space using (25)  The solution is computed at t = 0.2 using 401 grid points. The residual is discretized using the orders of accuracy described in Section 3.3. In Figs. 1a-1b the solution is presented for varying orders of accuracy using RV bdf and RV mg . Both RV methods provide significant shock-capturing for all orders of accuracy, with the maximum amount of over-and undershoot being 0.04% and 1.24% respectively for RV bdf . Comparatively the maximum amount of over-and undershoot is 0.65% and 1.31% for RV mg . In Fig. 2a the RV methods are compared to the solution obtained using only high-order upwind dissipation, where it is clear that the RV methods reduce the oscillations at the shock significantly. The comparison is done using a 5th order accurate spatial discretization. Fig. 2b shows the first-order viscosity ε 1 and residual viscosity ε r for the different RV methods. For both methods the residual viscosity is localized to the shock. In Fig. 3 the same problem setup is solved with RV bdf using a 5th order accurate spatial discretization for varying spatial resolutions. In addition, the relative errors in the l 1 -and l 2 -norm and corresponding convergence rates of the solutions, using RV bdf and RV mg , are presented in Table 1. The obtained convergence rates are close to the theoretical rate of 1 and 0.5 measured in the l 1 -and l 2 -norm respectively. Next, we illustrate the effect of the high-order artificial dissipation provided by the upwind SBP operators. This is done by comparing the solution to (25) discretized using a 5th order accurate SBP upwind discretization, and a 6th order central difference SBP discretization. In the central difference discretization the second-derivative D (ε) 2 is constructed using the variable coefficient narrow-stencil SBP second derivatives derived in Mattsson [27]. In both cases RV bdf is used to stabilize the shock. Using the initial condition (45) the solution is computed at t = 0.1, for 201 grid points. The results are presented in Fig. 4. It should be noted that in contrast to a standard wide-stencil central difference second derivative, the narrow-stencil second derivative has better dispersion properties, and thus is less sensitive to numerical oscillations. However, it is clear from Fig. 4 that the upwind discretization is still more effective in suppressing spurious oscillations away from the shock.
To conclude the section, we illustrate that using RV mg may be more computationally costly compared to using RV bdf , due to the additional cost of applying the interpolation and coarse grid discrete advective flux in contrast to applying the BDF stencil. To illustrate this, the wall-clock run time τ (in seconds) needed to evaluate R bdf in (18) and R mg in (19) for   the above problem is measured over 40 samples. The test is implemented in Matlab R2019b and performed on a 2,6 GHz 6-Core Intel Core i7 processor. The minimum and mean of τ is presented in Table 2. The results illustrate that for the present problem R mg is indeed more costly to evaluate compared to R bdf . Table 1 Relative errors in the l 1 -and l 2 -norm and convergence rates q for the solution to Burgers' equation in 1D at t = 0.2.

Analysis in 2D
In this section, the RV method is extended to 2D, and in particular it is demonstrated how to weakly impose boundary conditions for which linear stability estimates may be obtained. The stability estimates in Section 4 were based on the assumption that a skew-symmetric form of the advective flux exists. However, in many cases, such as for non-convex fluxes, a skew-symmetric form is not explicitly known or difficult to obtain. For such problems, linear energy estimates can still be obtained via the method of frozen coefficients. Assuming sufficient smoothness of the solution, the linear energy estimates are also valid for the non-linear problems. Although smoothness of the solution can only be guaranteed up to some finite time (given smooth initial data), the energy estimates are still useful in guiding the discrete boundary treatment. In practice, the boundary treatment may be used for non-smooth problems as well.
Extensions to more complex geometry may be done by considering curvilinear grid transformations and multi-block setups, see e.g., [36,39,2,1], but are outside the scope of the present study. We now assume that no skew-symmetric form of the advective flux is available, (which is the case e.g., for the non-convex conservation law studied in Section 7.3). For this reason, we turn to linear energy estimates, which are obtained by assuming sufficient smoothness of u such that the problem may be analyzed in its quasi-linear form Denoting the components of f by a 1 , and a 2 , i.e., f (x, u) = (a 1 (x, u), a 2 (x, u)), Dirichlet conditions on inflow boundaries are given by Now consider the corresponding frozen coefficient problem i.e., (47) where now a 1 and a 2 are constants. For the case of linear, but spatially variable coefficients, see e.g., Nordström [35], Mishra and Svärd [31]. Assume further that a 1 ≥ 0, a 2 ≤ 0, such that Dirichlet conditions should be imposed on the west and north boundaries. Other combinations of the signs of a 1 , a 2 , are treated similarly. Applying the energy method to the frozen coefficient problem (i.e., taking the inner product u, u t and integrating by parts) results in the estimate d dt where the boundary terms on the right-hand-side are given by Note that B T E ≤ 0 and B T S ≤ 0, contributing to energy dissipation. Assuming zero boundary data and integrating in time it follows that u(t) 2 ≤ u 0 2 and thus the frozen coefficient problem is energy stable. Assuming smoothness of the solution this implies energy stability of the non-linear IBVP (47), see Strang [41].

Semi-discrete setting
We now extend the RV method to 2D and illustrate how to formulate the boundary treatment using SAT leading to linear stability estimates, based on the above continuous analysis. To start with, the SBP operators in Section 2 are extended to 2D.

SBP operators in 2D
The domain is discretized by m x × m y grid points (x i , y j ) where

Any grid function v in 2D is then defined as vector of vectors
The benchmark problems in Section 7 are all defined on square domains, and thus to simplify notation m = m x = m y grid points per coordinate direction will be used from here on. Let I denote the unit matrix of size m × m and let A be the m 2 × m 2 diagonal matrix obtained by projection of a 2D grid function a onto the diagonal. The operators previously defined are extended to 2D as follows: e W = e l ⊗ I, e S = I ⊗ e l , e E = e r ⊗ I, In addition, the following notation for the solution vector and matrices on the boundary is introduced.
Note also that H = H due to the domain being square, and m x = m y . Furthermore the discrete L 2 -inner product on for 2D grid functions v, w is given by v, w H . The SBP properties of the 1D operators in Section 2 are now straightforwardly extended to 2D. Those relevant for the present work are presented below:

SBP-SAT discretization & stability analysis
An upwind SBP-SAT approximation of the 2D IBVP in (10) regularized with RV, is given by where (ε) + is given in (50), and the SAT imposes boundary conditions weakly. As stated in the continuous analysis in the beginning of Section 6, we are concerned with the stability of (52) for cases where no skew-symmetric form of the advective flux is available. Therefore sufficient smoothness of the solution is assumed such that the quasi-linear form of the advective flux, stated in (47) may instead be discretized. The upwind SBP approximation of the advective flux is given by The terms containing S x and S y are high-order upwind dissipative terms, obtained by Lax-Friedrichs flux splitting of f . See Section 4.1.1 for a more detailed presentation of flux-splitting in 1D. Let A 1,2 be the diagonal matrices obtained by projection of components a 1,2 (x, u) of f (u). SATs imposing the Dirichlet conditions in (48) are given by In order to obtain a stability estimate, artificial Neumann conditions must be specified on outflow boundaries (see the discussion in Section 4.1.2). SATs imposing the Neumann conditions εu x= = 0 are constructed as Note that the condition is consistent with setting no boundary condition on the outflow in the limit of grid refinement, since lim h→0 ε = 0 by (12), (13) and (15). Stability of (52) together with the SATs in (54)-(55) is given by the following lemma.
Using (51) it follows that the SBP property of D f (u) is given by For the following analysis, further assume a 1 > 0, and a 2 < 0, i.e., that the west and north boundaries are inflow boundaries. This implies that the SAT should be constructed as The cases involving other combinations of the signs of a 1 , a 2 follow analogously. Stability of the frozen coefficient scheme is shown via the discrete energy method, i.e., taking the inner product u, u t H and using (51) and (56).
where the diffusion from the RV is given by the dissipation from the upwind operators is given by is obtained. Note that S x and S y are negative semi-definite, such that −AD is positive semi-definite. The frozen coefficient scheme is therefore stable, and assuming that a 1,2 (x, u) are smooth the stability of the frozen coefficient scheme implies non-linear stability of the scheme (52) with (53), see Strang [41].
As in the 1D case, the stability estimate (58) is a semi-discrete counterpart of the continuous estimate u(t) 2 ≤ u 0 2 with additional control over the discrete derivatives D x+ u, D y+ u where ε is large via the low-order viscous term D I, as well as additional high-order artificial dissipation from the upwind term A D.

Computations in 2D
In this section, the capabilities of the RV method is demonstrated by applying it to a set of benchmark problems, described in more detail in the following sections. For all problems, time integration is performed using the algorithms presented in Section 3.4, where the solution vector is updated using the standard 4th order Runge-Kutta method and a time step k = 0.1h. Boundary conditions on inflow and outflow boundaries are imposed using (54) and (55) respectively. A compact way to automatically impose the correct type of boundary condition on inflow and outflow boundaries, is achieved by introducing switch functions. First, the components a 1,2 (x, u) of f are split using a Steger-Warming type splitting as Table 3 Relative errors and convergence rates q of the solution for the smooth 2D advection problem at t = 2π , measured in the The above SAT formulation has been used for the computations in the following sections. The residual is discretized using the orders of accuracy described in Section 3.3, unless otherwise stated.

Advection equation in 2D
In order to demonstrate that the RV is localized to regions with discontinuities, and that accuracy in smooth regions is maintained, consider the rotational 2D advection equation, i.e., the IBVP (10) with the advective flux given by on the domain = [−0.5, 0.5] 2 . The problem is discretized in space using (52), with the SBP upwind discrete advective flux D f (u) given by where B 1,2 are the matrices obtained by projecting β 1,2 (x) onto the diagonal, respectively. The SAT is constructed as in (59).
To start with, the advection equation is solved using the smooth initial condition u 0 (x) = 5 15,0). Since the problem is smooth, the accuracy of the solution should not be affected by the addition of RV. Solutions are therefore calculated using; high-order upwind dissipation without RV, RV bdf and RV mg .
The unmodified local wave speed w i = |f (u) i | is used to compute the first-order viscosity in (13). The relative errors in the l 2 -norm and corresponding convergence rates of the solutions are measured at t = 2π , after a full rotation. The results are presented in Table 3. On coarser meshes the RV methods are quite diffusive, but once refined the errors approach the ones obtained using only high-order upwind dissipation. Note that the interior accuracy of the operators is obtained due to the solution not interacting with the boundary. The convergence rates of the RV methods are somewhat higher than expected, but should approach the expected convergence rate once refined sufficiently.
Remark 7.1. For the computations with RV bdf presented in Table 3, BDF (6) is used. For a 9th order spatial discretization, this means that the formal accuracy of the solution is reduced to 8th order (see Section 3.3). The above results however show that no drop in convergence is observed for the utilized grid resolutions.
Next, to demonstrate that the RV is localized to regions with discontinuities, the problem is initialized with a non-smooth solution, consisting of a slotted cylinder, a cone and a smooth profile, as shown in Fig. 6a. For reference, the solution at t = 2π , without any dissipation is presented in Fig. 6b, while solutions using only first-order viscosity, and high-order upwind dissipation without any RV, are presented in Figs. 6c-6d. The solutions using RV bdf and RV mg are presented in Fig. 7. In all of the above cases a 5th order accurate spatial discretization is used on a 401 × 401 grid. The fine features of Fig. 6. Solutions to the advection equation in 2D with non-smooth initial data, at t = 2π without RV, using 5th order accurate spatial discretizations on a 401 × 401 grid. 6a Initial condition; 0 ≤ u ≤ 1, 6b no dissipation; −0.2878 ≤ u ≤ 1.3232, 6c first-order viscosity; 0 ≤ u ≤ 0.6394, 6d upwind; −0.1456 ≤ u ≤ 1.1484. Fig. 7. Solutions to the advection equation in 2D with non-smooth initial data, at t = 2π , using 5th order accurate spatial discretizations on a 401 × 401 grid. 7a RV bdf ; −0.0279 ≤ u ≤ 1.0234, 7b RV mg ; −0.0124 ≤ u ≤ 1.0119. the solutions are preserved using RV bdf and RV mg , while the over-and undershoot in the solutions are significantly smaller compared to the solution obtained using only upwind dissipation. In addition, the relative errors in the l 1 -and l 2 -norms and corresponding convergence rates using RV bdf and RV mg are computed. Similar to the results obtained in Section 5, the expected convergence rates of 1 and 0.5 are obtained for all tested orders of accuracy. The complete results are presented in Table 5 in Appendix B. For comparative studies in the spectral-and finite element settings, see e.g., Lu et al. [24], Guermond et al. [13].

Burgers' equation in 2D
In order to illustrate that the RV method handles problems with shocks of different magnitudes as well as rarefaction waves, consider the scalar Burgers' equation on the unit square = [0, 1] 2 , i.e., the IBVP (10) with the advective flux given For the initial condition x < 0.5, y < 0.5, 0.8, and Dirichlet conditions on inflow boundaries, the analytic solution is given by The problem is discretized in space using (52), with the SBP upwind discrete advective flux D f (u) given by Boundary conditions are imposed using (59), where the Dirichlet boundary data is given by evaluating (63) on the boundary.
Solutions are computed at t = 0.5 using RV bdf and RV mg . The modified local wave speed in (14) is used when computing the first-order viscosity. Fig. 8 shows the solutions u, together with the RV ε using RV bdf , RV mg and global first-order viscosity. The solutions are calculated using 5th order accurate spatial approximations on a 401 × 401 grid. Note that the rarefaction wave in the lower-right region of the solution is significantly less diffused using the RV methods, as are regions with discontinuities in the solution. The relative errors in the l 1 -and l 2 -norm and corresponding convergence rates are measured and, again, the expected rates of 1 and 0.5 are obtained for both RV methods, for all tested orders of accuracy. The complete results are presented in Table 6 in Appendix B. A comparative study where similar convergence rates are obtained is found in Guermond and Nazarov [9] using entropy viscosity with P 1 finite elements.

Remark 7.2.
For the 2D Burgers' equation, non-linear energy and stability estimates can be obtained using the energy method via skew-symmetric splitting of the advective flux in each coordinate direction, as done in the analysis presented in Section 4.1.3. Together with appropriately scaled SATs, this allows Lemmas 4.1-4.2 to be extended to 2D. The analysis is completely analogous to the 1D case and is therefore omitted.

Kurganov-Petrova-Popov rotating wave problem
In the final benchmark, it is demonstrated that the RV method is capable of handling non-convex advective fluxes.
The problem is discretized in space using (52), with the SBP upwind discrete advective flux D f (u) given by where boundary conditions are again imposed using (59). To illustrate the numerical difficulties in obtaining a solution to the KPP problem, the solution is first computed using high-order upwind dissipation without RV using a 5th order accurate spatial discretization and 401 2 grid points. The result is presented in Fig. 9. When solving the non-smooth 2D advection Fig. 9. KPP problem at t = 1 using a 5th order accurate spatial discretization on a 401 × 401 grid. The solution is computed using high-order upwind dissipation without RV.

Table 4
Over-and undershoot for the KPP problem using different order BDFs. 1.28 % 2.04% BDF (7) 0.92 % 1.59 % BDF (8) 0.80 % 0.81 % BDF (9) 0.46 % 0.44 % BDF (10) 0.47 % 0.38 % BDF (11) 0.28 % 0.30 % problem using only high-order upwind dissipation (Fig. 6d), the dominant errors are located to regions with discontinuities, where significant over-and undershoots occur. However, for the KPP problem the errors spread to regions away from the shock, and it is clear that even when disregarding oscillations at discontinuities, high-order upwind dissipation alone is not sufficient to capture the correct wave structure. The solutions obtained using the RV methods are presented in Figs. 10-11. For this problem |f (u)| = 1 and therefore the unmodified wave speed w i = |f (u) i | is used when computing the firstorder viscosity in (13). Some numerical artifacts may be observed, but overall the methods are able to capture the wave structure. Comparative studies are performed in e.g., Kurganov et al. [19] using an adaptive scheme based on Minmod 1 and WENO5, in Guermond and Nazarov [9] using entropy-viscosity with P 1 finite elements, and in Guermond et al. [13] using second-order maximum principle preserving Lagrange finite elements. The maximum overshoot in the solution presented in Guermond and Nazarov [9] is smaller than what is presented in the present work. However, it should be noted that higher-order discretizations are more prone to exhibit small-scale oscillations close to shocks. For RV bdf it is also possible to decrease the amount of over-and undershoot by increasing the order of the BDF, which is illustrated in Section 7.3.1.

Increasing the order of BDFs
The computations performed in the previous sections calculate RV bdf using BDF (6) . In the following, we illustrate that increasing the order of the BDFs may lead to less over-and undershoot in the solution. To illustrate this, the over-and undershoot of the solution compared to its minimum and maximum value, using different order BDFs is presented in Table 4. BDF:s of higher order than 6 are listed in Appendix A. The setup is the same as Section 7.3. From Table 4 it is clear that the over-and undershoot is decreased when the order of the BDF is increased. Fig. 12 shows the solution u together with the RV ε. The solution is computed using RV bdf with BDF (11) . Note that the viscosity at the shock is wider for BDF (11) compared to BDF (6) , see Fig. 12b and Fig. 10b respectively. One can also observe that the tail of the spiral for the BDF (6) case is shorter than the other cases leading to a larger overshoot for BDF (6) in that area, see Fig. 11a. Using BDF (11) would typically be too memory intensive for most large scale applications, but the results nonetheless indicate that there is more to investigate with respect to how the residual can be discretized.

Conclusion and future work
A high-order upwind FDM stabilized using RV has been presented for scalar conservation laws. Two ways of discretizing the residual using SBP operators have been proposed. The method allows for explicit time-stepping under a non-restrictive CFL condition. Numerical results demonstrate that the method is convergent in the l 1 -and l 2 -norm with expected rates, for both smooth and non-smooth problems, and is capable of handling both convex and non-convex fluxes. However, we Fig. 10. KPP problem at t = 1 using a 5th order accurate spatial discretization on a 401 × 401 grid. 10a-10b RV bdf ; solution 0.5607 ≤ u ≤ 11.1368, viscosity 0 ≤ ε ≤ 7.0711 · 10 −3 , 10c-10d RV mg ; solution 0.4110 ≤ u ≤ 11.2200, viscosity 0 ≤ ε ≤ 7.0711 · 10 −3 . Fig. 11. KPP problem at t = 1 using a 5th order accurate spatial discretization on a 401 × 401 grid. 11a RV bdf , 11b RV mg .
do not claim that the solution is overshoot and oscillation free. Additional work is needed to make sure that the method satisfies the maximum principle (Guermond et al. [13,14]) and is ongoing. Currently, the only convergence proof of the RV method present in the literature is for scalar conservation laws and is limited to time-stepping using implicit Euler (see Nazarov [32]). In Section 4.2 a coarse energy estimate in the fully discrete setting is obtained. An extension of the energy estimate to include non-viscous regions could potentially be a first step towards proving convergence of the RV method presented in this paper. Such an extension is currently being investigated. Finally, to make the presented method relevant for physical problems, we plan to extend it to systems of conservation laws and multi-block discretizations in an upcoming paper. Fig. 12. KPP problem at t = 1 using a 5th order accurate spatial discretization on a 401 × 401 grid using RV bdf , with BDF (11) ; solution 0.7520 ≤ u ≤ 11.0262, viscosity 0 ≤ ε ≤ 7.0711 · 10 −3 .

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. − 1524600u i−3 + 762300u i−2 − 304920u i−1 + 83711u i BDF (1) -BDF (6) can be found in many standard numerics textbooks e.g., Hairer et al. [15,III.1]. Table 5 Relative errors the in l 1 -and l 2 -norm and convergence rates q of the solution at t = 2π to the non-smooth 2D advection problem in Section 7.1.  Table 6 Relative errors in the l 1 -and l 2 -norm and convergence rates q for the solution at t = 0.5 to Burgers' equation in 2D in