Convergence of second-order, entropy stable methods for multi-dimensional conservation laws

TL;DR

This paper proves the convergence of a high-order entropy stable numerical method for multi-dimensional hyperbolic conservation laws, establishing theoretical guarantees for its accuracy and stability.

Contribution

It demonstrates how entropy stability leads to convergence in multiple dimensions, filling a gap in rigorous analysis of such methods.

Findings

Entropy stability bounds oscillations

Convergence to weak entropy solutions is proven

Provides theoretical foundation for high-order methods in multiple dimensions

Abstract

High-order accurate, $\textit{entropy stable}$ numerical methods for hyperbolic conservation laws have attracted much interest over the last decade, but only a few rigorous convergence results are available, particularly in multiple space dimensions. In this paper we show how the entropy stability of one such method yields a (weak) bound on oscillations, and using compensated compactness we prove convergence to a weak solution satisfying at least one entropy condition.

Equations

$$\{\begin{array}{ll}{\partial }_{t}u+\nabla \cdot f(u)=0& \forall (x,t)\in {\mathbb{R}}^d\times {\mathbb{R}}_{+}\\ u(x,0)=u_0(x)& \forall x\in {\mathbb{R}}^d\end{array}$$

$${\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{R}}_{+}}{\partial }_{t}u\phi +f(u)\cdot \nabla \phi dxdt+{\int }_{{\mathbb{R}}^d}{u}_0(x)\phi (x,0)dx=0\forall \phi \in C_c^1({\mathbb{R}}^d\times [0,\infty )).$$

$${\eta }_k(u):=|u-k|,q_k(u):=\mathrm{sign}(u-k)(f(u)-f(k)),$$

$${\partial }_t\eta (u)+\nabla \cdot q(u)=0.$$

$${\partial }_t\eta (u)+\nabla \cdot q(u)⩽0$$

$$\{\begin{array}{ll}{\partial }_{t}u+{\partial }_x{f}^x(u)+{\partial }_y{f}^y(u)=0& \forall (x,y,t)\in {\mathbb{R}}^2\times {\mathbb{R}}_{+}\\ u(x,y,0)=u_0(x,y)& \forall (x,y)\in {\mathbb{R}}^2.\end{array}$$

$$\{\{u\}{\}}_{i+\text{nicefrac}12,j}=\frac{u_{i,j}+u_{i+1,j}}{2}$$

$$\{\{u\}{\}}_{i,j+\text{nicefrac}12}=\frac{u_{i,j}+u_{i,j+1}}{2}$$

$$u_{i,j}(t)\approx \frac{1}{\Delta x\Delta y}{\int }_{{\mathcal{C}}_{i,j}}u(x,y,t)d(x,y).$$

$$\frac{d}{dt}{u}_{i,j}(t)+\frac{F_{i+\text{nicefrac}12,j}^x-F_{i-\text{nicefrac}12,j}^x}{\Delta x}+\frac{F_{i,j+\text{nicefrac}12}^y-F_{i,j-\text{nicefrac}12}^y}{\Delta y}=0,$$

$$\frac{d}{dt}\eta (u_{i,j})+\frac{Q_{i+\text{nicefrac}12,j}^x-Q_{i-\text{nicefrac}12,j}^x}{\Delta x}+\frac{Q_{i,j+\text{nicefrac}12}^y-Q_{i,j-\text{nicefrac}12}^y}{\Delta y}=0,$$

$$\frac{d}{dt}\eta (u_{i,j})+\frac{Q_{i+\text{nicefrac}12,j}^x-Q_{i-\text{nicefrac}12,j}^x}{\Delta x}+\frac{Q_{i,j+\text{nicefrac}12}^y-Q_{i,j-\text{nicefrac}12}^y}{\Delta y}⩽0,$$

$${\partial }_{t}u(v)+\nabla \cdot f(u(v))=0.$$

$$\frac{d}{dt}\eta (u_{i,j})+v_{i,j}\frac{F_{i+\text{nicefrac}12,j}^x-F_{i-\text{nicefrac}12,j}^x}{\Delta x}+v_{i,j}\frac{F_{i,j+\text{nicefrac}12}^y-F_{i,j-\text{nicefrac}12}^y}{\Delta y}=0.$$

$$\begin{array}{r}\frac{d}{dt}\eta (u_{i,j})+\frac{Q_{i+\text{nicefrac}12,j}^x-Q_{i-\text{nicefrac}12,j}^x}{\Delta x}+\frac{Q_{i,j+\text{nicefrac}12}^y-Q_{i,j-\text{nicefrac}12}^y}{\Delta y}\\ =\frac{r_{i+\text{nicefrac}12,j}^x+r_{i-\text{nicefrac}12,j}^x}{2\Delta x}+\frac{r_{i,j+\text{nicefrac}12}^y+r_{i,j-\text{nicefrac}12}^y}{2\Delta y}\end{array}$$

$$\begin{array}{rl}{r}_{i+\text{nicefrac}12,j}^x& =[[v]{]}_{i+\text{nicefrac}12,j}{F}_{i+\text{nicefrac}12,j}^x-[[{\psi }^x]{]}_{i+\text{nicefrac}12,j},\\ r_{i,j+\text{nicefrac}12}^y& =[[v]{]}_{i,j+\text{nicefrac}12}{F}_{i,j+\text{nicefrac}12}^y-[[{\psi }^y]{]}_{i,j+\text{nicefrac}12}\\ Q_{i+\text{nicefrac}12,j}^x& =\{\{v\}{\}}_{i+\text{nicefrac}12,j}{F}_{i+\text{nicefrac}12,j}^x-\{\{{\psi }^x\}{\}}_{i+\text{nicefrac}12,j},\\ Q_{i,j+\text{nicefrac}12}^y& =\{\{v\}{\}}_{i,j+\text{nicefrac}12}{F}_{i,j+\text{nicefrac}12}^y-\{\{{\psi }^y\}{\}}_{i,j+\text{nicefrac}12}.\end{array}$$

$$F_{i+\text{nicefrac}12,j}^x={\tilde{F}}_{i+\text{nicefrac}12,j}^x-D_{i+\text{nicefrac}12,j}^x[[v]{]}_{i+\text{nicefrac}12,j},F_{i,j+\text{nicefrac}12}^y={\tilde{F}}_{i,j+\text{nicefrac}12}^y-D_{i,j+\text{nicefrac}12}^y[[v]{]}_{i,j+\text{nicefrac}12}$$

$$[[v]{]}_{i+\text{nicefrac}12,j}{\tilde{F}}_{i+\text{nicefrac}12,j}^x=[[{\psi }^x]{]}_{i+\text{nicefrac}12,j},[[v]{]}_{i,j+\text{nicefrac}12}{\tilde{F}}_{i,j+\text{nicefrac}12}^y=[[{\psi }^y]{]}_{i,j+\text{nicefrac}12}$$

$$\frac{d}{dt}\eta (u_{i,j})+\frac{Q_{i+\text{nicefrac}12,j}^x-Q_{i-\text{nicefrac}12,j}^x}{\Delta x}+\frac{Q_{i,j+\text{nicefrac}12}^y-Q_{i,j-\text{nicefrac}12}^y}{\Delta y}$$

$$=\frac{D_{i+\text{nicefrac}12,j}^x[[v]{]}_{i+\text{nicefrac}12,j}^2+D_{i-\text{nicefrac}12,j}^x[[v]{]}_{i-\text{nicefrac}12,j}^2}{2\Delta x}+\frac{D_{i,j+\text{nicefrac}12}^y[[v]{]}_{i,j+\text{nicefrac}12}^2+D_{i,j-\text{nicefrac}12}^y[[v]{]}_{i,j-\text{nicefrac}12}^2}{2\Delta y}.$$

$$\begin{array}{r}{w}_{i+\text{nicefrac}12,j}^{\pm }=\mathcal{R}{w}_{\cdot ,j}(x_{i+\text{nicefrac}12}\pm 0),w_{i,j+\text{nicefrac}12}^{\pm }=\mathcal{R}{w}_{i,\cdot }(y_{j+\text{nicefrac}12}\pm 0)\end{array}$$

$$⟨⟨w⟩{⟩}_{i+\text{nicefrac}12,j}=w_{i+\text{nicefrac}12,j}^{+}-w_{i+\text{nicefrac}12,j}^{-},⟨⟨w⟩{⟩}_{i,j+\text{nicefrac}12}=w_{i,j+\text{nicefrac}12}^{+}-w_{i,j+\text{nicefrac}12}^{-}.$$

$$\begin{array}{r}{F}_{i+\text{nicefrac}12,j}^x={\tilde{F}}_{i+\text{nicefrac}12,j}-D_{i+\text{nicefrac}12,j}^x⟨⟨\tilde{v}⟩{⟩}_{i+\text{nicefrac}12,j},\\ F_{i,j+\text{nicefrac}12}^y={\tilde{F}}_{i,j+\text{nicefrac}12}-D_{i,j+\text{nicefrac}12}^y⟨⟨\tilde{v}⟩{⟩}_{i,j+\text{nicefrac}12},\end{array}$$

$$[[\tilde{v}]{]}_{i+\text{nicefrac}12,j}⟨⟨\tilde{v}⟩{⟩}_{i+\text{nicefrac}12,j}⩾0,[[\tilde{v}]{]}_{i,j+\text{nicefrac}12}⟨⟨\tilde{v}⟩{⟩}_{i,j+\text{nicefrac}12}⩾0.$$

$$\mathrm{sign}⟨⟨w⟩{⟩}_{i+\text{nicefrac}12}=\mathrm{sign}[[w]{]}_{i+\text{nicefrac}12}.$$

$$\frac{⟨⟨w⟩{⟩}_{i+\text{nicefrac}12}}{[[w]{]}_{i+\text{nicefrac}12}}⩽C_p.$$

$$\begin{array}{r}{F}_{i+\text{nicefrac}12,j}^x={\tilde{F}}_{i+\text{nicefrac}12,j}-D_{i+\text{nicefrac}12,j}^x⟨⟨u⟩{⟩}_{i+\text{nicefrac}12,j},\\ F_{i,j+\text{nicefrac}12}^y={\tilde{F}}_{i,j+\text{nicefrac}12}-D_{i,j+\text{nicefrac}12}^y⟨⟨u⟩{⟩}_{i,j+\text{nicefrac}12},\end{array}$$

$$\underline{D}⩽D_{i+\text{nicefrac}12,j}^x,D_{i,j+\text{nicefrac}12}^y⩽\overline{D}\text{for fixed }\underline{D},\overline{D}>0$$

$$\int_{0}^{T}\sum_{i,j}\Bigl{(}\bigl{|}[[u]]_{{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}},j}\bigr{|}^{3}{\Delta y}+\bigl{|}[[u]]_{i,{j+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}}\bigr{|}^{3}{\Delta x}\Bigr{)}\,dt\leqslant C\|u\|_{L^{\infty}(\mathbb{R}^{2}\times[0,T])},$$

$$\int_{0}^{T}\sum_{i,j}\Big{(}[[u]]_{{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}},j}\langle\langle u\rangle\rangle_{{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}},j}{\Delta y}+[[u]]_{i,{j+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}}\langle\langle u\rangle\rangle_{i,{j+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}}{\Delta x}\Big{)}\,dt\leqslant C.$$

$$\frac{1}{2}\sum _{i,j}(u_{i,j}(T){)}^2\Delta x\Delta y-\frac{1}{2}\sum _{i,j}(u_{i,j}(0){)}^2\Delta x\Delta y=-\mathcal{E},$$

$$\displaystyle\mathcal{E}:=\int_{0}^{T}\sum_{i,j}\Bigl{(}D^{x}_{{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}},j}[[u]]_{{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}},j}\langle\langle u\rangle\rangle_{{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}},j}{\Delta y}+D^{y}_{i,{j+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}}[[u]]_{i,{j+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}}\langle\langle u\rangle\rangle_{i,{j+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}}{\Delta x}\Bigr{)}\,dt.$$

$$\sum_{i\in\mathbb{Z}}\bigl{|}[[w]]_{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}\bigr{|}^{3}\leqslant 2\|w\|_{l^{\infty}}\sum_{i\in\mathbb{Z}}\langle\langle w\rangle\rangle_{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}[[w]]_{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}$$

Full text

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 642768 for NC’s work.

Keywords: Multi-dimensional conservation laws; TECNO scheme; entropy stability; finite volume methods

1. Introduction

Hyperbolic conservation laws appear in a large variety of applications, including gas dynamics, traffic modeling, multi-phase fluid flow problems, and more; see e.g. [2, 7, 10, 13]. We consider a scalar, $d$-dimensional hyperbolic conservation law

[TABLE]

where $u=u(x,t):\mathbb{R}^{d}\times\mathbb{R}_{+}\to\textbf{U}$ is the unknown conserved variable, taking values in some nonempty interval $\mathbf{U}\subset\mathbb{R}$, and the function $f=(f_{1},\dots,f_{d}):\textbf{U}\to\mathbb{R}^{d}$ is the smooth (at least $C^{2}$ on U) and possibly nonlinear flux function. It is well-known that even if the initial datum $u_{0}(x)$ is arbitrarily smooth, the solutions of (1.1) may still be non-smooth [2, 7, 10]. Thus, it is fruitless to look for solutions of (1.1) in the classical sense. Instead these solutions are sought in a weak sense. A function $u\in L^{\infty}(\mathbb{R}^{d}\times\mathbb{R}_{+})$ is said to be a weak solution of (1.1) if it is a distributional solution, i.e.

[TABLE]

It is well known (see e.g. [2, 7, 10]) that weak solutions may be non-unique. Thus to single out a physically relevant solution, the notion of weak solution has to be supplemented with an additional admissibility criterion, namely entropy conditions. A pair of functions $\eta:\mathbb{R}\to\mathbb{R},q:\mathbb{R}\to\mathbb{R}^{d}$ is an entropy pair for (1.1) if $\eta$, the entropy function, is convex and $q$, the entropy flux function, satisfies $q^{\prime}(u)=\eta^{\prime}(u)f^{\prime}(u)$. In particular, for every $k\in\mathbb{R}$ we have the well-known Kružkov entropy pair $(\eta_{k},q_{k})$ given by

[TABLE]

see [11]. Multiplying (1.1) by $\eta^{\prime}(u)$ and using the chain rule we have the following entropy conservation identity for smooth solutions of (1.1),

[TABLE]

Due to the possible non-smoothness of solutions of (1.1), the above derivation cannot be rigorously justified for weak solutions. Instead, the entropy inequality

[TABLE]

is imposed. As was shown by Kružkov [11], this entropy condition guarantees uniqueness and stability of solutions of (1.1).

1.1. Numerical methods for conservation laws

The nonlinear nature of the equation (1.1) and the fact that its solutions are irregular, can make the construction and analysis of numerical methods for (1.1) challenging. We outline here some of the available literature on this subject.

In order to converge to a weak solution satisfying the entropy condition (1.5), the numerical method needs to satisfy some discrete version of the entropy condition. Finite volume methods (to be discussed in subsection 2.1 in detail) with this property are called entropy stable. Harten, Hyman and Lax [8] showed that all monotone schemes for scalar conservation laws are entropy stable with respect to any entropy pair $(\eta,q)$. Osher [14] generalized this to a (presumably) larger class of schemes, the so-called E-schemes. Osher also showed that these E-schemes are at most first-order accurate. In his papers from $1984$ [16] and from $1987$ [17], Tadmor laid a general framework for constructing entropy stable schemes by first constructing entropy conservative methods – schemes satisfying a discrete version of (1.4) – and then adding numerical diffusion to obtain entropy stability. As he proved in [17], entropy conservative schemes are generally second-order accurate; even higher-order entropy conservative schemes were constructed by Lefloch, Mercier and Rohde in [12]. However, the addition of numerical diffusion to any of these entropy conservative schemes, in the way suggested in [16, 17], degrades the order of accuracy to 1.

By combining the high-order accurate entropy conservative schemes in [17, 12] with a judiciously chosen reconstruction method, Fjordholm, Mishra and Tadmor [5, 3] constructed entropy stable methods with an arbitrarily high order of accuracy, the so-called TECNO schemes. By estimating the amount of entropy dissipated by the method (i.e., the amplitude of the left-hand side in (1.5)), the authors could derive a priori weak regularity bounds on the numerical solution, and these bounds were sufficient to prove convergence of the method in the special case of $d=1$ space dimensions.

To the best of our knowledge there is no available proof of convergence of a high-order accurate entropy stable method for a multi-dimensional conservation law. The purpose of the present paper is to prove convergence for a particular case, namely the second-order TECNO scheme.

2. Entropy stable numerical methods

2.1. Finite volume methods

For the sake of notational simplicity we are going to consider the scalar conservation law (1.1) in the particular case of $d=2$ space dimensions, although we emphasize that the results in this paper are equally valid for any number of spatial dimensions $d$.

We write (1.1) in the case $d=2$ as

[TABLE]

Here and in the remainder we will denote the components of all vector-valued functions by $f=(f^{x},f^{y})$.

One of the most popular discretization frameworks is finite volume methods. The spatial domain $\mathbb{R}^{2}$ is partitioned into rectangles of the form $\mathcal{C}_{i,j}=[x_{i-{\nicefrac{\mathrm{1}}{\mathrm{2}}}},x_{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}})\times[y_{{j-{\nicefrac{\mathrm{1}}{\mathrm{2}}}}},y_{j+{\nicefrac{\mathrm{1}}{\mathrm{2}}}})$, where we for the sake of simplicity use uniform grid sizes $x_{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}-x_{i-{\nicefrac{\mathrm{1}}{\mathrm{2}}}}\equiv{\Delta x}$ and $y_{j+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}-y_{j-{\nicefrac{\mathrm{1}}{\mathrm{2}}}}\equiv{\Delta y}$. We denote the midpoint values as $x_{i}=\frac{x_{i-{\nicefrac{\mathrm{1}}{\mathrm{2}}}}+x_{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}}{2}$ and $y_{j}=\frac{y_{j-{\nicefrac{\mathrm{1}}{\mathrm{2}}}}+y_{j+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}}{2}$. For any quantity $(u_{i,j})_{i,j\in\mathbb{Z}}$ defined on this grid, we define the jump and average operators

[TABLE]

We let $u_{i,j}(t)$ be an approximation of the average value of $u$ over the rectangles $\mathcal{C}_{i,j}$, that is,

[TABLE]

The initial data is sampled as $u_{i,j}(0)=\frac{1}{{\Delta x}{\Delta y}}\int_{\mathcal{C}_{i,j}}u_{0}(x,y)\,d(x,y)$. A semi-discrete finite volume method for (2.1) can then be written in the generic form

[TABLE]

where the numerical flux function $F^{x}_{{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}},j}=F^{x}(u_{i,j},u_{i+1,j})$ is computed from an approximate solution of the Riemann problem at the interface $\{(x_{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}},y)\}_{y_{j-{\nicefrac{\mathrm{1}}{\mathrm{2}}}}\leqslant y\leqslant y_{j+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}}$ and $F^{y}_{i,{j+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}}$ is computed analogously [7, 10]. The computed solution generated by the scheme is given by $u^{\mathbf{\Delta}}(x,y,t)=\sum_{i,j}u_{i,j}(t)\chi_{\mathcal{C}_{i,j}}(x,y)$, where ${\mathbf{\Delta}}=({\Delta x},{\Delta y})$ and $\chi_{\mathcal{C}}$ is characteristic function for the rectangle $\mathcal{C}$. We say that the numerical flux function $F$ is consistent with $f$ if $F^{x}(u,u)=f^{x}(u)$ and $F^{y}(u,u)=f^{y}(u)$ for all $u\in\mathbf{U}$. We also say that a numerical flux $F$ is locally Lipschitz continuous if $F^{x},F^{y}$ are locally Lipschitz continuous in each argument.

2.2. Entropy stable numerical methods

In order for any limit $u=\lim_{{\mathbf{\Delta}}\to 0}u^{\mathbf{\Delta}}$ to satisfy the entropy condition (1.5), the numerical method (2.2) must satisfy some discrete form of the entropy condition. In this section we briefly review the theory of so-called entropy conservative and entropy stable schemes, and we define the TECNO schemes, which will be the subject of the rest of the paper.

Definition 2.1** (Entropy conservative methods).**

Let $(\eta,q)$ be an entropy pair. We say that the finite volume method (2.2) is entropy conservative (with respect to $(\eta,q)$) if computed solutions satisfy the discrete entropy equality

[TABLE]

where $Q^{x}_{{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}},j}=Q^{x}(u_{i,j},u_{i+1,j})$ and $Q^{y}_{i,{j+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}}=Q^{y}(u_{i,j},u_{i,j+1})$ are numerical entropy flux functions satisfying $Q^{x}(u,u)=q^{x}(u)$ and $Q^{y}(u,u)=q^{y}(u)$ for all $u\in\mathbf{U}$.

Definition 2.2** (Entropy stable methods).**

Let $(\eta,q)$ be an entropy pair. We say that the finite volume method (2.2) is entropy stable (with respect to $(\eta,q)$) if computed solutions satisfy the discrete entropy equality

[TABLE]

where $Q^{x}_{{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}},j}=Q^{x}(u_{i,j},u_{i+1,j})$ and $Q^{y}_{i,{j+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}}=Q^{y}(u_{i,j},u_{i,j+1})$ are numerical entropy flux functions satisfying $Q^{x}(u,u)=q^{x}(u)$ and $Q^{y}(u,u)=q^{y}(u)$ for all $u\in\mathbf{U}$.

For an entropy pair $(\eta,q)$ the mapping $u\mapsto\eta^{\prime}(u)$ is of particular importance, and we denote this entropy variable by $v=v(u):=\eta^{\prime}(u)$. If $\eta$ is strictly convex, $\eta^{\prime\prime}(u)>0$, then the map $u\mapsto v(u)$ is strictly monotone increasing and hence is invertible. This inverse will be denoted by $u(v)$. Thus, the mapping $u\mapsto v$ induces a change of variables, in terms of which we can pose the conservation law (2.2) as

[TABLE]

We define also the entropy potential $\psi:\mathbf{U}\to\mathbb{R}^{d}$ defined by $\psi(u):=v(u)f(u)-q(u)$, whose name is given by the fact that $\partial_{v}\psi(u(v))=f(u(v))$.

A general approach to designing entropy conservative/stable schemes is as follows. Multiplying both sides of (2.2) by $v_{i,j}:=\eta^{\prime}(u_{i,j})$ and using the chain rule we get

[TABLE]

Adding and subtracting terms yields

[TABLE]

where

[TABLE]

It is straightforward to see that $Q^{x},Q^{y}$ are consistent with $q$ in the sense of Definitions 2.1 and 2.2, as long as $F^{x},F^{y}$ are consistent with $f$. Thus, if $F^{x},F^{y}$ are chosen such that either $r\equiv 0$ or $r\leqslant 0$, then the scheme (2.2) is entropy conservative/stable. In particular, if $F^{x},F^{y}$ are of the form

[TABLE]

for nonnegative coefficients $D^{x},D^{y}\geqslant 0$ and numerical fluxes $\tilde{F}^{x},\tilde{F}^{y}$ satisfying

[TABLE]

then the resulting scheme (2.2) is entropy stable. These observations were first made by Tadmor [16, 17]; see also [18]. For fluxes of the form (2.8) we also get a precise expression for the amount of entropy dissipated in (2.4):

[TABLE]

Under further assumptions on $\eta$ and $D$, this yields explicit bounds on “weak TV” terms of the form $\sum_{i,j}[[v]]_{{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}},j}^{2}{\Delta y}$, which can be used to prove compactness and convergence of the numerical method; see e.g. [3]. We will apply this approach in Section 3.

Remark 2.3**.**

The above observations can be used to design entropy stable schemes, by first finding numerical fluxes $\tilde{F}^{x},\tilde{F}^{y}$ satisfying (2.9), and then adding diffusion in the form (2.8). We note that with this approach, we are only guaranteed that the discrete entropy inequality (2.4) (or (2.3) for entropy conservative schemes) is satisfied for one particular entropy pair $(\eta,q)$.

2.3. The TECNO scheme

The scheme (2.2) with fluxes $\tilde{F}^{x},\tilde{F}^{y}$ satisfying (2.9) is entropy conservative, in the sense of Definition 2.1. It can be shown that two-point entropy conservative schemes are at most second-order accurate, in the sense of truncation error [16, 17]. When adding diffusion in the form (2.8) with $D=O(1)$, the resulting scheme is at most first-order accurate. The TECNO schemes, introduced in [3, 5], represent a systematic approach to designing higher-order accurate entropy stable schemes. Since the convergence proof in Section 3 only applies to the second-order TECNO schemes, we will only describe these methods here, and we refer to [3, 5] for the general construction.

The TECNO scheme has two main ingredients: An entropy conservative flux $\tilde{F}^{x},\tilde{F}^{y}$, and a sign preserving reconstruction method. Since our mesh is a Cartesian grid, we define the reconstruction procedure in a tensorial manner. For a partition $\{\mathcal{C}_{i}\}_{i\in\mathbb{Z}}$ of $\mathbb{R}$ we consider a $p$th order reconstruction operator $\mathcal{R}$, mapping any grid function $(w_{i})_{i\in\mathbb{Z}}$ to a piecewise $(p-1)$th order polynomial $\mathcal{R}w(x)$. Multi-dimensional grid functions $(w_{i,j})_{i,j\in\mathbb{Z}}$ are reconstructed dimension-by-dimension, defining in particular the edge values

[TABLE]

Fix now some entropy pair $(\tilde{\eta},\tilde{q})$. The second-order TECNO scheme [5, 3] is constructed from a flux $\tilde{F}^{x},\tilde{F}^{y}$ which is entropy conservative with respect to $(\tilde{\eta},\tilde{q})$, and applies a second-order reconstruction method to the entropy variables $\tilde{v}=\tilde{\eta}\circ u$. The resulting scheme (2.2) has numerical flux

[TABLE]

where $D^{x},D^{y}\geqslant 0$. As shown in [5], the above scheme is formally second-order accurate, and it satisfies the discrete entropy inequality (2.4) for the entropy pair $(\tilde{\eta},\tilde{q})$, provided the reconstruction operator $\mathcal{R}$ satisfies the sign property

[TABLE]

This is indeed true for the ENO reconstruction method [9]:

Theorem 2.4** (The ENO sign property [6]).**

For some $p\in\mathbb{N}$, let $\mathcal{R}$ denote the $p$-th order ENO reconstruction operator. Then for any grid function $(w_{i})_{i\in\mathbb{Z}}$,

[TABLE]

Moreover, there exists a constant $C_{p}>0$ depending only on $p$ such that

[TABLE]

For the sake of simplicity we henceforth select the entropy for which TECNO is entropy stable as $\tilde{\eta}(u)=u^{2}/2$. The corresponding entropy variable is then simply $\tilde{v}=u$, making the mapping between conserved and entropy variables somewhat easier. A summary of the TECNO scheme that we will analyze in this paper follows:

Definition 2.5**.**

The second-order TECNO scheme for (2.1) is the numerical scheme (2.2) with the numerical flux

[TABLE]

where

• •

$\tilde{F}$* is a consistent and locally Lipschitz continuous numerical flux which is entropy conservative with respect to the entropy $\tilde{\eta}(u)=u^{2}/2$*

• •

the diffusion coefficients $D$ satisfy

[TABLE]

• •

$\langle\langle u\rangle\rangle_{{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}},j}$* and $\langle\langle u\rangle\rangle_{i,{j+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}}$ denote the jumps in the second-order ENO reconstruction of the conserved variable $(u_{i,j}(t))_{i,j\in\mathbb{Z}}$.*

Theorem 2.6**.**

The second-order TECNO scheme (cf. Definition 2.5) has the following properties:

1. (i)

it is entropy stable with respect to the square entropy $\tilde{\eta}(u)=u^{2}/2$ 2. (ii)

the flux $F$ is locally Lipschitz continuous 3. (iii)

there is some $C>0$ independent of ${\Delta x},{\Delta y}$ such that

[TABLE]

Proof.

The entropy stability follows from the calculations in Section 2.2 and the sign property (2.12). Local Lipschitz continuity of $F$ follows from the Lipschitz continuity of $\tilde{F}$ and the upper bound (2.13).

For the “weak TV bounds” (2.15), summing (2.6) over $i,j\in\mathbb{Z}$, integrating over $t\in[0,T]$ and using the specific form of $F$ in (2.14) yields

[TABLE]

Since $D\geqslant 0$ and the reconstruction satisfies the sign property (2.12), we have $\mathcal{E}\geqslant 0$. From the above we also see that $\mathcal{E}\leqslant\frac{1}{2}\sum_{i,j}(u_{i,j}(0))^{2}\,{\Delta x}{\Delta y}\leqslant\frac{1}{2}\|u_{0}\|_{L^{2}(\mathbb{R}^{2})}^{2}$, so we see that the left-hand side of (2.15b) can be upper-bounded by $\underline{D}\mathcal{E}\leqslant\frac{1}{2}\underline{D}\|u_{0}\|_{L^{2}(\mathbb{R}^{2})}^{2}<\infty$.

For the remaining property (2.15a) we use the following fact, proved in [4, Section 4.4]: For every grid function $(w_{i})_{i\in\mathbb{Z}}$ with compact support,

[TABLE]

where $\langle\langle w\rangle\rangle_{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}=w_{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}^{+}-w_{i+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}^{-}$ denotes the jump in the second-order ENO reconstruction of $w$. Thus, the left-hand side of (2.15a) can be bounded by $\|u\|_{L^{\infty}}$ times the left-hand side of (2.15b). ∎

Remark 2.7**.**

For even higher-order TECNO schemes, the results in Theorem 2.6 are still valid, with the exception of (2.15a): The crucial estimate (2.16) has been conjectured but remains unproven; cf.  [4, Section 4.4] or [3, Section 5.5].

3. Convergence of the scheme

Given a numerical solution $(u_{i,j}(t))_{i,j\in\mathbb{Z},t\in\mathbb{R}_{+}}$ computed with the second-order TECNO scheme (cf. Definition 2.5) we define the piecewise constant function

[TABLE]

where ${\mathbf{\Delta}}=({\Delta x},{\Delta y})$. The goal of this section will be to show the following theorem:

Theorem 3.1**.**

Assume that the solution $u^{\mathbf{\Delta}}$ computed by the TECNO scheme (cf. Definition 2.5) is uniformly $L^{\infty}$ bounded,

[TABLE]

for some $M>0$. Then there is some subsequence ${\mathbf{\Delta}}^{\prime}\to 0$ such that $u^{{\mathbf{\Delta}}^{\prime}}\to u$ pointwise a.e. and in $L^{p}(\mathbb{R}^{2}\times[0,T])$ for every $p\in[1,\infty)$. The function $u$ is a weak solution of (2.1) which satisfies the entropy condition (1.5) for the entropy $\eta(u)=u^{2}$.

We will use the method of compensated compactness, and we give the main results required here in Section 3.1. The convergence proof is given in Section 3.2, but we summarize it here:

Proof of Theorem 3.1.

The compactness result, Corollary 3.4 requires the entropy residuals $\{\partial_{t}\eta(u^{\mathbf{\Delta}})+\nabla\cdot q(u^{\mathbf{\Delta}})\}_{{\mathbf{\Delta}}>0}$ to lie in a compact subset of $H^{-1}_{\text{loc}}$. Lemma 3.5 bounds their discrete equivalents $\partial_{t}\eta(u^{\mathbf{\Delta}})+\nabla\cdot Q(u^{\mathbf{\Delta}})$ by terms which, by Theorem 2.6(iii) and (3.1), are bounded in the sense of measures. Lemma 3.6 shows that the remainder $\nabla\cdot(q-Q)$ is small in $H^{-1}_{\text{loc}}$. We then conclude (using Murat’s Lemma) that the entropy residuals are precompact, and hence there is some strongly convergent subsequence.

Lemma 3.7 is a standard “Lax–Wendroff” proof, showing that the limit is a weak solution, and Lemma 3.8 shows consistency with a single entropy condition. ∎

3.1. Compensated compactness

We briefly summarize the technical compactness lemmas here, and refer to [1, 19] for more details.

Lemma 3.2** (Murat’s Lemma).**

Let $\Omega\subset\mathbb{R}^{d}$, $d\geqslant 2$ be an open, bounded subset. Let $(\mu)_{n\in\mathbb{N}}$ be a bounded sequence in $W^{-1,p}(\Omega)$ for some $2<p\leqslant\infty$. Suppose also that $\forall\;n\in\mathbb{N}$

[TABLE]

where $\xi_{n}$ lies in a compact set of $H^{-1}_{\text{loc}}(\Omega)$ and $\pi_{n}$ lies in a bounded set of $\mathcal{M}_{\text{loc}}(\Omega)$. Then $(\mu_{n})_{n\in\mathbb{N}}$ lies in a compact subset of $H^{-1}_{\text{loc}}(\Omega)$.

Theorem 3.3** (Panov, Theorem 5 in [15]).**

Let $(u_{\varepsilon})_{\varepsilon>0}$ be a bounded sequence in $L^{\infty}(\mathbb{R}^{d}\times\mathbb{R}_{+})$ such that for every $k\in\mathbb{R}$, the set

[TABLE]

is precompact in $H^{-1}_{\text{loc}}(\mathbb{R}^{d}\times\mathbb{R}_{+})$. (Here, $(\eta_{k},q_{k})$ denote the Kruzkov entropy pairs (1.3).) Then there is a subsequence $\varepsilon_{n}\to 0$ as $n\to\infty$ and a function $u\in L^{\infty}(\mathbb{R}^{d}\times\mathbb{R}_{+})$ such that

[TABLE]

The following corollary shows that it is enough to consider smooth entropies in the above result.

Corollary 3.4**.**

Let $(u_{\varepsilon})_{\varepsilon>0}$ be a bounded sequence in $L^{\infty}(\mathbb{R}^{d}\times\mathbb{R}_{+})$ such that for every entropy pair $(\eta,q)$ with $\eta\in C^{2}_{b}(\mathbb{R})$, the set

[TABLE]

is precompact in $H^{-1}_{\text{loc}}(\mathbb{R}^{d}\times\mathbb{R}_{+})$. Then there is a subsequence $\varepsilon_{n}\to 0$ as $n\to\infty$ and a function $u\in L^{\infty}(\mathbb{R}^{d}\times\mathbb{R}_{+})$ such that

[TABLE]

Proof.

For $k\in\mathbb{R}$ and $m\in\mathbb{N}$, let $\eta_{k,m}\in C^{2}(\mathbb{R})$ be convex functions satisfying

[TABLE]

Then the corresponding $C^{2}$ entropy flux $q_{k,m}(u)=\int_{k}^{u}f^{\prime}(s)\eta_{k,m}^{\prime}(s)\,ds$ also satisfies $\|q_{k,m}-q_{k}\|_{L^{\infty}(\mathbb{R})}\to 0$ as $m\to\infty$. Fix $k\in\mathbb{R}$ and decompose

[TABLE]

where

[TABLE]

Let $(\varepsilon_{n})_{n\in\mathbb{N}}$ be a sequence converging to 0 as $n\to\infty$ and fix an open, bounded subset $U\subset\mathbb{R}^{d}\times\mathbb{R}_{+}$. It is enough to prove the existence of a further subsequence $E_{\varepsilon^{\prime}_{n}}$ which converges in $H^{-1}(U)$. For every $\varphi\in H_{0}^{1}(U)$ we have

[TABLE]

Let $(\varepsilon_{n}^{\prime})_{n\in\mathbb{N}}$ be a subsequence of $(\varepsilon_{n})_{n\in\mathbb{N}}$, taken such that $(E_{m,\varepsilon_{n}^{\prime}}^{1})_{n\in\mathbb{N}}$ is convergent in $H^{-1}(U)$ for every $m\in\mathbb{N}$. In particular, we may assume that $\|E_{m,\varepsilon_{n}^{\prime}}^{1}-E_{m,\varepsilon_{k}^{\prime}}^{1}\|_{H^{-1}(U)}\leqslant 1/m$ for every $k,n\geqslant m$. Then

[TABLE]

as $n,k\geqslant m$ and $m\to\infty$. ∎

3.2. Convergence of TECNO

The TECNO scheme (Definition 2.5) is guaranteed to dissipate the square entropy $\tilde{\eta}(u)=u^{2}/2$, but the discrete entropy residual (2.6) might have either sign. We can nonetheless show that the entropy residual is not too large, in the following sense:

Lemma 3.5**.**

Assume that the solution computed by the TECNO scheme is $L^{\infty}$ bounded, (3.1). Then for any entropy pair $(\eta,q)$ with $\eta\in C^{2}$, the total discrete entropy production is upper-bounded by

[TABLE]

where $\nabla\cdot Q$ denotes the measure whose integral of any $\varphi\in C_{c}^{0}(\mathbb{R}^{2}\times\mathbb{R}_{+})$ is

[TABLE]

and where $Q^{x},Q^{y}$ are given by (2.7).

Proof.

Let $\tilde{\psi}$ be the entropy potential with respect to the square entropy $\tilde{\eta}=u^{2}/2$, $\tilde{\psi}(u)=uf(u)-\tilde{q}(u)$. We split the entropy residual $r$ in (2.7) as $r=r^{1}+r^{2}$, where

[TABLE]

and similarly for $r^{y}_{i,{j+{\nicefrac{\mathrm{1}}{\mathrm{2}}}}}$. The first part of the entropy residual can be estimated as

[TABLE]

and similarly in the $y$-direction,

[TABLE]

The second part of the entropy residual can be bounded by

[TABLE]

The conclusion now follows. ∎

We can now show precompactness of the sequence of approximations:

Lemma 3.6**.**

Let $\Omega\subset\mathbb{R}^{2}\times[0,T]$ be a bounded subset and assume that the solution computed by the TECNO scheme is $L^{\infty}$ bounded, (3.1). Then there is a subsequence ${\mathbf{\Delta}}^{\prime}\to 0$ such that $u^{{\mathbf{\Delta}}^{\prime}}\to u$ pointwise a.e. and in $L^{p}_{\text{loc}}(\mathbb{R}^{2}\times\mathbb{R}_{+})$ for $1\leqslant p<\infty$, for some $u\in L^{1}\cap L^{\infty}(\mathbb{R}^{2}\times\mathbb{R}_{+})$.

Proof.

Let $(\eta,q)$ be an arbitrary $C^{2}$ entropy pair. By Corollary 3.4 it is sufficient to show that the sequence $\mathcal{E}_{{\Delta x},{\Delta y}}:=\partial_{t}\eta(u^{\mathbf{\Delta}})+\nabla\cdot q(u^{\mathbf{\Delta}})$ is precompact in $H^{-1}_{\text{loc}}$, and to this end we will employ Murat’s lemma. Firstly note that $\mathcal{E}_{{\Delta x},{\Delta y}}$ is bounded in $W^{-1,\infty}(\mathbb{R}^{2}\times\mathbb{R}_{+})$, by the $L^{\infty}$ bound on $u^{\mathbf{\Delta}}$. Decompose

[TABLE]

where $Q$ is given in (2.7). Note that, due to the $L^{\infty}$ bound on $u^{\mathbf{\Delta}}$ and Theorem 2.6, also $Q$ is locally Lipschitz continuous. By Lemma 3.5 and Theorem 2.6(iii), the discrete entropy production $\mathcal{E}^{1}$ is bounded in the space of measures $\mathcal{M}(\mathbb{R}^{d}\times\mathbb{R}_{+})$.

Now to show that $\mathcal{E}^{2}$ is precompact in $H^{-1}_{\text{loc}}(\mathbb{R}^{2}\times[0,T])$, let $\Omega\subset\mathbb{R}^{2}\times[0,T]$ be open and bounded and let $\varphi\in H^{1}_{0}(\Omega)$. Extending $\varphi$ by zero outside $\Omega$, we get

[TABLE]

where the last step follows from 2.15. Thus by invoking Murat’s Lemma 3.2 we can conclude that the sequence $(\mathcal{E}_{{\Delta x},{\Delta y}})_{{\Delta x},{\Delta y}>0}$ is precompact in $H^{-1}_{\text{loc}}(\mathbb{R}^{2}\times[0,T])$. Applying Corollary 3.4 then yields the desired result. ∎

Now we need to show that this limit function $u$ is indeed a weak solution of (2.1). To do so we state and prove the following “Lax–Wendroff result”.

Lemma 3.7**.**

Under the same assumptions as in Lemma 3.6, the limit $u$ is a weak solution of (2.1).

Proof.

Let $\varphi\in C^{1}_{c}(\mathbb{R}^{2}\times(0,T))$ be a test function and select a compact set $K_{x}\times K_{y}\subset\mathbb{R}^{2}$ such that $\operatorname{supp}\varphi\subset K_{x}\times K_{y}\times[0,T]$. Furthermore, denote

[TABLE]

Multiplying the numerical scheme (2.2) by $\varphi_{i,j}(t)$ and integrating/summing in time/space, we get

[TABLE]

After performing integration and summation by parts for temporal and spatial variables respectively we get

[TABLE]

where we can write

[TABLE]

We can write $A^{1}=-\int_{0}^{T}\int_{\mathbb{R}}\int_{\mathbb{R}}u^{\mathbf{\Delta}}\partial_{t}\varphi^{\mathbf{\Delta}}\,dx\,dy\,dt$, and thanks to the convergence of $u^{\mathbf{\Delta}}$ to $u$ from Lemma 3.6 and convergence of $\varphi^{\mathbf{\Delta}}$ to $\varphi$ a.e., we have $\lim_{{\Delta x},{\Delta y}\to 0}A^{1}=-\int_{0}^{T}\int_{\mathbb{R}^{2}}u\partial_{t}\varphi\,dx\,dy\,dt$.

For the second term $A^{2}$, we denote for the sake of simplicity $\Delta_{x}\psi(x,y,t)=\frac{\psi(x+{\Delta x},y,t)-\psi(x,y,t)}{{\Delta x}}$, for any function $\psi$. Since $\tilde{F}^{x}$ is a two-point flux, we can write

[TABLE]

where

[TABLE]

Thanks to the convergence of $u^{\mathbf{\Delta}}$ from Lemma 3.6 and the a.e. convergence of $\varphi^{\mathbf{\Delta}}$ to $\varphi$ we have

[TABLE]

For the term $A^{2,2}$ we apply the Hölder inequality and Lemma 3.6 to get

[TABLE]

as ${\Delta x},{\Delta y}\to 0$ by (2.15a). Analogously, $A^{3}\to-\int_{0}^{T}\int_{\mathbb{R}^{2}}f^{y}(u)\partial_{y}\varphi\,dx\,dy\,dt$ as ${\Delta x},{\Delta y}\to 0$. We conclude that

[TABLE]

It remains to show that $B^{1},B^{2}$ in (3.10) vanish as ${\Delta x},{\Delta y}\to 0$. Indeed,

[TABLE]

as ${\Delta x},{\Delta y}\to 0$, and likewise for $B^{2}$. This completes the proof. ∎

Although we are not able to show that the TECNO scheme converges to the entropy solution, we will show that the weak solution $u$ satisfies at least one of the entropy conditions.

Lemma 3.8**.**

With the same assumptions as in Lemma 3.6, the limit $u$ satisfies

[TABLE]

Proof.

As in (2.6) in Lemma 3.5 we can write

[TABLE]

where, in this particular case, the entropy residuals $r$ on the right-hand side are all nonpositive (see e.g. [16, 5, 3]). Multiplying the above by a nonnegative test function $\varphi\in C^{1}_{c}(\mathbb{R}^{2}\times(0,T))$ and proceeding in the same manner as in Lemma 3.7 we obtain (3.13) in the sense of distribution. ∎

4. Conclusions and outlook

We prove convergence of the second-order TECNO scheme in two space dimensions to a weak solution of the hyperbolic conservation law (1.1); this can easily be generalized to any number of space dimensions. The proof of this result relies on estimating the entropy residual appropriately using the (weak) TV bound obtained from entropy stability with respect to one entropy. Invoking this estimate, precompactness of the sequence of approximate solutions is shown using a corollary derived from a compensated compactness result due to Panov. Finally, to show that the limit function obtained due to the precompactness property is indeed a weak solution of (1.1), a “Lax–Wendroff” type argument is used.

Convergence proofs of even higher-order (i.e. more than second order) TECNO scheme in multiple space dimensions, to a weak solution of the equation (1.1) are still unanswered. In our opinion, this is largely due to the unavailability of weak TV estimates of the type (2.15a), as well as an appropriate version of Lemma 3.5. This should be an object of interest for future research. Last, but not least, one key estimate to prove (2.15b), and consequently (2.15a), is (2.16). For even higher-order ENO reconstruction, this estimate (2.16) (the “ENO-conjecture”) is still not established and is still an open problem.