    # On the TVD property of second order methods for 2D scalar conservation laws

The total variation diminishing (TVD) property is an important tool for ensuring nonlinear stability and convergence of numerical solutions of one-dimensional scalar conservation laws. However, it proved to be challenging to extend this approach to two-dimensional problems. Using the anisotropic definition for discrete total variation (TV), it was shown in <cit.> that TVD solutions of two-dimensional hyperbolic equations are at most first order accurate. We propose to use an alternative definition resulting from a full discretization of the semi-discrete Raviart-Thomas TV. We demonstrate numerically using the second order discontinuous Galerkin method that limited solutions of two-dimensional hyperbolic equations are TVD in means when total variation is computed using the new definition.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

We consider two-dimensional hyperbolic scalar conservation laws

 ut+f(u)x+g(u)y=0, (1)

where and are the flux components in the - and -directions, respectively. Weak solutions of (1) might develop discontinuities in finite time even with smooth initial data. This is one of the challenges in designing robust numerical methods for these equations. Modern methods, which are formally at least second order accurate, develop spurious oscillations near solution discontinuities that need to be controlled in order for the solution to remain stable. At the same time, suppressing oscillations should not substantially reduce solution accuracy in the rest of the domain. For one-dimensional problems, this issue has been largely resolved by enforcing the total variation diminishing (TVD) property . Let be a numerical solution to (1) at . A one-step numerical scheme is called total variation diminishing (TVD) if satisfies the following condition

 TV(Un+1)≤TV(Un),∀n. (2)

Let be a bounded open set and . Then the total variation is defined as 

 TV(u)=supφ∈C1c(Ω){−∫Ωu∇⋅φ dxdy:∥φ(x,y)∥≤1,∀(x,y)∈Ω}, (3)

where and

is a differentiable vector function. For sufficiently smooth functions

, e.g. , where is a Sobolev space, the definition (3) reduces to

 TV(u)=∫Ω√u2x+u2y dxdy. (4)

Several definitions of discrete TV, i.e. TV for functions defined on a computational grid, have been proposed in the literature. Let us assume that is a rectangular domain discretized into a Cartesian grid of elements with centroids at points , and the grid sizes in the and directions and , respectively. Let be a discrete function defined on this grid with being the value associated with element .

In their classical work, Goodman and LeVeque  used discrete TV of the form

 TVa(U)=∑i,j(Δyj|Ui+1,j−Ui,j|+Δxi|Ui,j+1−Ui,j|), (5)

commonly referred to in the literature as anisotropic TV. The definition (5) is a discrete approximation of

 TVa(u)=∫Ω|ux|+|uy| dxdy. (6)

They showed that any conservative, TVD scheme for solving hyperbolic conservation laws in the two-dimensional space is at most first-order accurate for monotone solutions of (1). In order to prove this statement, an associated one-dimensional scheme with the same order of accuracy and special initial data was considered. If the two-dimensional scheme is TVD, then the corresponding one-dimensional scheme should be monotone at least on certain initial data and, therefore, is at most first-order accurate. Thus, the original two-dimensional scheme is also at most first order accurate.

The TVD property is an important tool for proving nonlinear stability and convergence of a numerical scheme for one-dimensional problems. The limiters are commonly used to develop TVD schemes, e.g. [19, 24] for finite volume methods and  for DG methods. In the absence of a suitable TVD property for higher order two-dimensional schemes, most proposed and currently in use limiters are ad hoc. Recently, the local maximum principle (LMP) has been employed in proving stability of numerical solutions in the maximum norm [10, 11, 12, 13, 18]. The LMP property that the solution averages on a computational element at the next time level will not exceed the bounds given by solution averages at the present time in some neighborhood of the element. Thus, the solution is stable in the maximum norm for all times. Unfortunately, stability in the maximum norm is not sufficient to guarantee convergence, though enforcing both linear stability and LMP is property resulted in a convergent scheme . Therefore, it is desirable to develop a TVD definition, that once enforced on solutions of high-order schemes, will guarantee stability and convergence.

While no two-dimensional numerical scheme can be both TVD and second order accurate in the sense of (5), this statement was not proven for all possible definitions of discrete TV, e.g. a widely used

 TVis(U)=∑i,j√Δy2j(Ui+1,j−Ui,j)2+Δx2i(Ui,j+1−Ui,j)2, (7)

which was first introduced in , and is commonly referred to in the literature as isotropic TV. The analytical TV defined by (3) is translation and rotation-invariant. However, the discrete TV definitions (5) and (7) are not isotropic, i.e. if the shape given by is rotated, the corresponding values of , can increase or decrease, see  and Section 3. This makes these definitions not suitable for our purposes, since the TV of the numerical solution of (1) will depend on the orientation of with respect to the grid. Both (5) and (7) are based on the forward difference approximations of , and while (7) improves the accuracy of (5), it remains spatially anisotropic. A comprehensive discussion of this issue can be found in .

A novel approach to computing discrete TV has been recently proposed in . Given a mesh, the authors considered a staggered grid approximation of the divergence operator in (3) and solved an optimization problem to find . This approach, mimicking the standard, by duality, definition of the total variation, is referred to as the discrete dual TV. Since is a convex functional for , the optimization problem resulting from discretization of (3) has a unique global maximum, and can be solved by many standard proximal algorithms for convex optimization. The convergence of discrete dual TV to (3) holds in both weak and strong topologies, e.g. for any , [3, 5]. Moreover, for a broad family of discrete dual TV  any sequence of bounded in discrete functions, with uniformly bounded dual TV contains a convergent subsequence. A fast proximal gradient descent algorithm for solving the optimization problem arising from the implementation of this formulation was proposed in  and further improved in 

. The new dual discrete TV led to improvement in accuracy of TV computation and was used to solve TV-regularization based optimization problems with applications to image denoising, inpainting, motion estimation, and multi-label image segmentation

[20, 7, 1, 21, 4].

The main features distinguishing duality based TV definitions from the classical (5) and (7) are spatial isotropy and asymptotic consistency with (3). Among multiple dual discrete definitions, the one proposed in  is the most “isotropic” according to . However, it is difficult to define the isotropy of a discrete TV, since there is no unique way of defining rotation and translation on a mesh.

In this paper, we consider ,

, and the novel dual definition and apply them to compute TV in means of numerical solutions of two-dimensional hyperbolic equations. In particular, we employ the second order discontinuous Galerkin method with the moment limiter

. We demonstrate that the computed numerical solutions are not TVD in the and sense, as expected. In contrast, the means of limited numerical solutions remain TVD when the dual definition is used. We conjecture that similar to one-dimensional problems, common limiters that are known in practice to stabilize solutions, reduce the slope of the solution enough for it to be TVD in the dual sense. Thus, we propose to use the dual TV definition as a suitable measure of multi-dimensional discrete TV for solutions of PDEs.

The paper is organized as follows. Section 2 introduces the duality-based definitions of discrete total variation and the algorithms used to compute it. Section 3 contains two numerical tests, which are used to compare the isotropy, accuracy, and consistency of the discrete TVs described above. The main numerical results are presented in Section 4 where we use the discontinuous Galerkin method to compute solutions to (1) and their TVs. Finally, Section 5 presents our conclusions.

## 2. Duality-based discrete TV definitions

In this section we introduce duality-based definitions of discrete TV. Let have a bounded distributive derivative , and let a test function be a differentiable vector field. We will further assume that , and on , where is the unit normal vector along the boundary . Then by the divergence theorem the following identity holds

 ∫Ωu∇⋅φ dxdy=−∫ΩφDu dxdy. (8)

Note that can be replaced by if is sufficiently smooth . Using (8) in (3), the total variation of on can be defined as

 TV(u)=supφ∈C1c(Ω){∫ΩφDu dxdy:∥φ(x,y)∥≤1, ∀(x,y)∈Ω}. (9)

We can approximate the left hand side of (8) as

 ∫Ωu∇⋅φ dxdy≈∑i,jui,j(φi+1/2,j−φi−1/2,jΔxi+ψi,j+1/2−ψi,j−1/2Δyj)ΔxiΔyj, (10)

where , , , , and .

Assuming for simplicity that the computational grid is square, i.e. , we apply summation by parts to the right hand side of (10) to obtain a discrete analogue of the right hand side of (8)

 Δx∑i,jui,j((φi+1/2,j−φi−1/2,j)+(ψi,j+1/2−ψi,j−1/2))= (11) −Δx∑i,jφi+1/2,j(ui+1,j−ui,j)+ψi,j+1/2(ui,j+1−ui,j)=−Δx∑i,j⟨Dui,j,φi,j⟩,

where denotes the Euclidean inner product of vectors in , and and are defined below. The vector with the components

 D1ui,j=ui+1,j−ui,j,D2ui,j=ui,j+1−ui,j, (12)

can be viewed as a forward difference approximation of the gradient of at the centroid of , up to division by . Alternatively, can be viewed as a centered approximation of the partial derivative of with respect to at , the midpoint of the right edge of , and as the partial derivative with respect to at , the upper edge’s midpoint, up to division by . Similarly, the values and are combined into a vector, . Note that though the values of and are computed at edge midpoints, we associate with and summation over and in (11).

Replacing the centroid values with , the values of the discrete function on each , we obtain a semi-discrete version of (9), known as Raviart-Thomas TV  or discrete dual total variation

 TVRT(U)=maxφ∈C1c(Ω){Δx∑i,j⟨DUi,j,φi,j⟩:∥φ(x,y)∥≤1, ∀(x,y)∈Ω}. (13)

and the factor before the sum accounts for the size of the cells. Since computation of the inner products in (13) requires the values of the continuous function at edge midpoints only, can be replaced with a discrete function , with defined on the grid of edge midpoints , . To obtain a fully discrete expression for , the constraint on should be replaced with an equivalent constraints on values of . In [3, 7, 16], this idea has been extensively studied and used to construct fully discrete dual TV definitions for TV-regularization based optimization with application to imaging problems.

There are multiple ways to impose the bound on the norm of the discrete test function . An obvious constraint results from the bounds on the values of at edge midpoints, i.e. and . However, by the derivation above and are not included in . Since these values are not available, they need to be defined outside of the definition (13

). We interpolate them by averaging. For example, the value of

on a uniform grid can be obtained by averaging , as

 ψi+1/2,j=ψi,j+1/2+ψi,j−1/2+ψi+1,j+1/2+ψi+1,j−1/24,

see Figure 1 and Figure 2 (right). We can write this in operator notation by setting where

 (P1˜φ)i+1/2,j=((P1˜φ)1i+1/2,j,(P1˜φ)2i+1/2,j)=(φi+1/2,j,ψi,j+1/2+ψi,j−1/2+ψi+1,j+1/2+ψi+1,j−1/24). (14) Figure 1. The stencil of the discrete test function in the dual definition (17) on Ωi,j=[xi−1/2,xi+1/2]×[yj−1/2,yj+1/2]. Components of ˜φ are shown in black. Interpolated values are shown in red and blue. Figure 2. Interpolation stencils for φi,j+1/2 (left), φi,j,ψi,j (center), ψi+1/2,j (right). Components of ˜φ are shown in black, interpolated values are shown in red and blue.

In , the first component is the identity operator and the second component averages the values of on the four horizontal edges around the point and assigns this value to , see Figure 2 (right). Similarly, we define , where the first component is the average of the values of on the four vertical edges around the point and the second component is the identity operator, see Figure 2 (left). Finally, we define the centroid value , as an average of edge values in the horizontal and vertical directions, see Figure 2 (center). The operators and are defined as

 (P2˜φ)i,j+1/2=((P2˜φ)1i,j+1/2,(P2˜φ)2i,j+1/2)=(φi+1/2,j+φi+1/2,j+1+φi−1/2,j+φi−1/2,j+14,ψi,j+1/2), (15) (P3˜φ)i,j=((P3˜φ)1i,j,(P3˜φ)2i,j)=(φi+1/2,j+φi−1/2,j2,ψi,j+1/2+ψi,j−1/22). (16)

Since was assumed to satisfy on , we require , which means that the boundary values of are equal zero, i.e. , . With that, we define the space

 Π(Ω)={˜φ∈R2×(N+1)×(N+1):(˜φ⋅→n)=0 on ∂Ω}.

Then, using the constraints and notations developed above, we arrive at a fully discrete expression for the dual total variation 

 TVd(U)=max˜φ∈Π(Ω){Δx∑i,j⟨DUi,j,˜φi,j⟩: √φ2i+1/2,j+ψ2i+1/2,j≤1, √φ2i,j+1/2+ψ2i,j+1/2≤1, √φ2i,j+ψ2i,j≤1, ∀(i,j)}, (17)

where the subscript stands for “dual”. Using the operator notation (14) - (16), the above can be rewritten as

 TVd(U)=max˜φ∈Π(Ω){Δx∑i,j⟨DUi,j,˜φi,j⟩:∥Pk˜φ∥2≤1, k=1,…,3, ∀(i,j)}, (18)

where denotes the Euclidean norm of a vector in .

Computing the maximizer of the constrained optimization problem (18) is a difficult task. Therefore, an equivalent saddle-point formulation  is derived. For this purpose, we will state the primal-dual formulation of (18). We introduce the gradient field , with approximating at , and , respectively. Let , where

 V={v:∥v∥1,1,2<+∞}, ∥v∥1,1,2=∑k∥vk∥1,2, ∥vk∥1,2=∑i,j∥(vk)i,j∥2=∑i,j√(vxk)2i,j+(vyk)2i,j.

We define the operator , with components

 (F1v)i+1/2,j=((vx1)i+1/2,j+(vx2)i,j+1/2+(vx2)i,j−1/2+(vx2)i+1,j+1/2+(vx2)i+1,j−1/24+(vx3)i,j+(vx3)i+1,j2), (19) (F2v)i,j+1/2=((vy2)i+1/2,j+(vy1)i+1/2,j+(vy1)i+1/2,j+1+(vy1)i−1/2,j+(vy1)i−1/2,j+14+(vy3)i,j+(vy3)i,j+12). (20)

The operator is a projection of onto a coarser grid of edge midpoints. The value of the first component is assigned to the edge midpoint and the value of the second component to the edge midpoint . Then, the primal-dual formulation of (18) can be stated as the following saddle-point problem

 TVd(U)=minv∈Vmax˜φ∈Π(Ω){Δx∑i,j⟨(Fv)i,j,˜φi,j⟩:Fv=DU, ∥Pk˜φ∥2≤1, k=1,2,3, ∀(i,j)}. (21)

Using (21), we can derive a minimization problem for , an approximation of the gradient of ,

 TVd(U)=minv∈V{Δx(∑i,j∥(v1)i,j∥2+∑i,j∥(v2)i,j∥2+∑i,j∥(v3)i,j∥2):Fv=DU}. (22)

The above is called the primal problem, for which (18) is the corresponding dual problem.

There are generally infinitely many ways to define and corresponding . A proper choice of operators allows one to establish the uniqueness of the maximizer in (18) and consistency of the resulting discrete dual TV definition with the analytical TV (3). We now appeal to a general theorem of , which states

###### Theorem 2.1.

Assume the supports and the weights of the convolutions defining operators on the meshes are uniformly bounded. Then discrete dual TV defined by

 TV(U)=minv{Δx∥v∥:Fv=DU}, (23)

where the norm is taken in an appropriate space, -converges to

 TV(U)={|Du|(Ω)if u∈BV(Ω),+∞otherwise . (24)

The convergence holds in both weak and strong topologies, such as for any .

The following result is the direct implication of Theorem 2.1

###### Corollary 2.1.1.

The discrete TV defined by (17) is consistent with (3) in the sense of (24).

Furthermore, the compactness result for follows from Proposition 2.5 of . Let be a sequence of discrete functions, defined on grids. Let us further assume that

 TVd(UN)<+∞, ∀N>0,

and that remains bounded in . Then there exists a subsequence of and a function , such that converges to in .

By Proposition 1 of , a strong duality between (18) and (22) holds. Therefore, the primal and dual problems have the same optimal value. That is, if a maximizer of the primal problem and a minimizer of the dual problem exist, then we have

 TVd(U)=Δx∑i,j⟨DUi,j,˜φ†i,j⟩=Δx∥v†∥1,1,2.

We compute the minimizer using the alternating proximal gradient method (, Algorithm 2), a simplified version of the general alternating direction method of multipliers. The algorithm is given below.

In Algorithm 1, are defined by (14)-(16) and is given by (12). The algorithm converges when and , where is the operator norm of and . Suitable values for the parameter are discussed in Section 3, while is used as suggested in .

## 3. Accuracy of discrete TV definitions

Before applying the three definitions of discrete TV introduced in the previous sections to numerical solutions of hyperbolic PDEs, we investigate how well discrete TVs approximate the analytical TV (3) under mesh refinement and verify convergence of the Algorithm 1.

### 3.1. Asymptotic consistency

We consider a rectangular domain and a function . is discretized into an mesh of square elements . We approximate by a grid-based function , where is set equal to the cell average of on . and are computed using (5) and (7), respectively. is computed according to Algorithm 1 with , which is a sufficient accuracy for the examples considered here and later on in Section 4. The TV values obtained using the three definitions are listed in Table 1. Next, we compare them to the analytical TV of given by (4)

 TV(u)=∫Ω√u2x+u2y dxdy≈0.835249.

It appears that and converge to as increases, while approaches the value of the integral in (6), i.e. at the limit of , for this smooth .

We report the magnitude of the difference between the analytical and discrete TVs, denoted by , in Table 1. We observe that for this smooth function, the and converge linearly to , and converges linearly to .

Remark. Convergence of Algorithm 1 depends on a proper choice for the value of parameter . In particular, we have observed in our numerical experiments that dependents on the size of a mesh and for a carefully tuned value of the convergence of can be superlinear. To find a suitable value for a mesh of a particular resolution, we ran a direct search in the interval with a step size . We choose the value of that provides the smallest . The number of steps to convergence also depends on mesh size. For example, it took iterations to achieve the accuracy with . The values of found for each mesh are listed in the last column of Table 1. These values were also used in the numerical tests in Section 4.

### 3.2. Spatial isotropy

Next, we investigate how discrete TV of a discontinuous shape varies when the shape is rotated. We consider a rectangular domain and a square pulse defined on by

 u(x,y)={1,if x∈[−1/√2,1/√2],y∈[−1/√2,1/√2],0,otherwise. (25)

Similar to the previous example, is set equal to the cell average of on , an element in a uniform, mesh of with . To simplify calculations, we choose to be a multiple of . The resulting discrete function with is illustrated in Figure 3 (A). takes three distinct values: in the interior of the square, on the elements containing the edges of and on the elements containing the corners. Here, denotes the fractional part of a real number. It follows from (5) that

 TVa(U)=4(2⌊N/(4√2)⌋+2{N/(4√2)})Δx,

where is the floor function of a real number. Under mesh refinement, i.e. as , tends to . Using (7), it is straightforward to show that also converges to . Finally, we evaluate numerically using Algorithm 1. We tabulate the values of TV according to the three definitions in Table 2. We observe that , , and converge to under mesh refinement.

Next, we consider the function , which corresponds to counterclockwise rotation of by the angle . Since rotation does not change analytical TV of according to (3), we have . Choosing , we obtain a square shape whose diagonals are aligned with the coordinate axes. Using the same meshes as above, we construct discrete functions from (Figure 3 (B)). takes three values: in the exterior of the square, in the interior, and on the elements containing edges of . The explicitly calculated value of the is

 TVa(V)=(2N−4)Δx=8−16/N,

which tends to . Then, is greater than by a factor of . The value for the

 TVis(V)=((3√2+2)N/4−5√2/2+2)Δx,

converges to , , which is greater than the by a factor of .

We observe that , , on the same mesh. This is expected because and are different projections of on a finite mesh. However, we would like to see convergence between these values under mesh refinement. TV values according to the three definitions are listed in Table 2. We notice that is the only one such that , i.e. the difference in TV values between the original and rotated functions, diminishes as increases.

We conclude that both and values change when the shape is rotated, while does not. Thus, neither , nor are isotropic.

## 4. Numerical Examples

In this section, we present a number of numerical examples, whose aim is to demonstrate that the second order DG method is TVD in the sense and is not TVD in the and

sense. In discretizing the spatial variable, we used the tensor-product, orthogonal basis functions of degree one. The system of ODEs resulted from the DG spatial discretization was solved with the Heun’s time integrator. We employed the moment limiter, introduced in

. All problems were solved on a square domain discretized into square cell, meshes.

We computed the total variation of the obtained solutions using the conventional , and dual (17) definitions, where is the solution average on . was calculated with Algorithm 1 and . The value of parameter is given in Table 1.

### 4.1. Rotation of a hill

We begin by studying rotation of a cosine hill around the origin described by

 ut+2πxux−2πyuy=0, (26)

with the initial condition

 (27)

and suitable boundary conditions. We solve the problem with and without the limiter until the final time , which corresponds to the counter-clockwise rotation of the hill by . The initial condition and solution at on the mesh are shown in Figure 4. Figure 4. Limited solution for (26) with initial condition (27) at t = 0.125 (left) and its contour plot (right).

We compute the errors and convergence rates at and list them in Table 3. We observe that both limited and unlimited solutions exhibit the second rate of convergence. Figure 5. TV for solutions of (26) with the initial condition (27) for the first 10 (N=80) and 20 (N=160) time steps, as a function of time t.

We compute the total variation of the solutions and report the results in Table 4 and Figure 5. In Figure 5, we plot TV values of solutions computed with and without the limiter for the first 10 and 20 time steps on and meshes, respectively. We observe that the total variation of the limited solutions monotonically decreases with each time step for all definitions of TV. The solution on the coarser mesh is more diffusive and is more affected by the limiter. Consequently, its TV decreases faster.

The TV of the unlimited solutions behaves differently for , and . initially oscillates about the exact value of one and then decays. This is expected of unlimited high-order methods, with TV changes on the order of discretization error from one time step to another. Moreover, a similar behavior was reported for one-dimensional problems in . The and initially increase and then decay, but the initial increase is a magnitude larger that in .

We further report TV for the unlimited solutions at ten time instances between and in Table 4. For this problem with a radial symmetry of the initial profile, we observe that all three TV share similar long-term behaviour. We also notice that is not a good approximation of value, as was discussed in Section 3. Figure 6. TV for the limited solutions of (26) with the initial condition (28) for the first 10 (N=80) and 20 (