1 Introduction

We consider the Euler system governing the ideal motion of a compressible gas on a bounded domain \(\Omega \subset {\mathbb {R}}^d \, (d=1,2,3)\), i.e.

$$\begin{aligned} \partial _t {\varvec{U}}+ \mathrm{div}_{ x}{\varvec{F}}({\varvec{U}}) = 0, \quad (t,x) \in (0,T) \times \Omega . \end{aligned}$$
(1.1)

Here \({\varvec{U}}:=(\varrho , {\varvec{m}}, E)^T\) denotes the conservative vector with the fluid density \(\varrho \), momentum \({\varvec{m}}\) and total energy E, while \({\varvec{F}}\) is the flux function given by

$$\begin{aligned} {\varvec{F}}= ({\varvec{m}}, {\varvec{u}}\otimes {\varvec{m}}+ p {\mathbb {I}}, {\varvec{u}}(E+p))^T. \end{aligned}$$
(1.2)

Moreover for positive \(\varrho \), \({\varvec{u}}= \frac{{\varvec{m}}}{\varrho }\) is the velocity of the fluid and p is the pressure satisfying the equation of state

$$\begin{aligned} p = (\gamma -1)\varrho e, \quad \gamma \in (1,2] \end{aligned}$$
(1.3)

with the specific internal energy \(e= \frac{E}{\varrho } -\frac{1}{2}|{\varvec{u}}|^2 \).

We close the system with initial data

$$\begin{aligned} {\varvec{U}}(0,x) ={\varvec{U}}_0 = (\varrho _0, {\varvec{m}}_0:=\varrho _0 {\varvec{u}}_0, E_0)^T \end{aligned}$$
(1.4a)

satisfying

$$\begin{aligned} \varrho _0 > 0 \quad \text{ and } \quad E_0 \in L^1(\Omega ) \end{aligned}$$
(1.4b)

and the impermeability boundary condition

$$\begin{aligned} {\varvec{u}}\cdot {\varvec{n}}|_{\partial \Omega } = 0, \end{aligned}$$
(1.5)

where \({\varvec{n}}\) is the outer normal vector on the boundary \(\partial \Omega \). Taking the second law of thermodynamics into account we further require that the entropy inequality holds, i.e.

$$\begin{aligned} \partial _t \, \eta ({\varvec{U}}) + \mathrm{div}_{ x}\,{\varvec{q}}({\varvec{U}}) \ge 0 \end{aligned}$$
(1.6)

with the physical entropy pair \((\eta , {\varvec{q}})\) defined by

$$\begin{aligned} \eta = C_v \varrho S , \quad {\varvec{q}}= \eta {\varvec{u}}\quad \text{ with } C_v = \frac{1}{\gamma -1} \text{ and } S= \ln \left( \frac{p}{\varrho ^\gamma }\right) . \end{aligned}$$
(1.7)

It is well-known that for the multidimensional Euler system there may exist infinitely many weak entropy solutions, i.e. the solutions satisfying (1.1)–(1.6) in the weak sense, cf. De Lellis and Székelyhidi [6], Chiodaroli et al. [2, 3], and Feireisl et al. [8]. As a consequence, the convergence analysis of standard numerical schemes for the Euler equations and identification of physically reasonable limiting solutions are of fundamental importance.

Over the past few decades there is a rapid development of efficient numerical methods for the Euler equations, cf. Toro [29], Feistauer et al. [14], Li et al. [21], Godunov [15], Shu and Osher [24] and LeVeque [20]. Despite the success in practical simulations, a rigorous convergence analysis of the numerical methods still remains open in general. Most results on error analysis were focused on scalar conservation laws. For example, Kuznetsov [19] showed that the (upper) \(L^1\)-error bound is \({\mathcal {O}}(h^{1/2})\) for multi-dimensional scalar conservation laws under the assumptions on the boundedness of the total variation and continuity in time of numerical solutions, where h is the mesh parameter. Further, Cockburn et al. [4] and Vila [30] extended the result of Kuznetsov and obtained the \(L^1\)-error bounds of \({{{\mathcal {O}}}}(h^{1/4})\) without the assumptions of bounded total variation and continuity in time.

Concerning the linear advection equation, Tang and Teng [26] showed the sharpness of the \({\mathcal {O}}(\sqrt{\Delta x})\) \(L^1\)-error for monotone difference schemes with BV initial data. For the nonlinear scalar equation Teng and Zhang [28] showed the optimal convergence rate of 1 in the \(L^1\)-norm for the viscosity method and monotone schemes if a solution is piecewise constant with finitely many shocks. Moreover, for the piecewise smooth entropy solution with finitely many rarefaction waves, Tang and Teng [27] showed that the error of viscosity solution to the inviscid solution is bounded by \({\mathcal {O}}(\varepsilon |\log \varepsilon | + \varepsilon )\) in the \(L^1\)-norm, where \(\varepsilon \) denotes the viscosity coefficient. Furthermore, Tadmor and Tang [25] studied the pointwise error estimates and showed that the thicknesses of the shock and rarefaction layers are of order \({\mathcal {O}}(\varepsilon )\) and \({\mathcal {O}}(\varepsilon \log ^2 \varepsilon )\), respectively. We point out that the error estimates for scalar hyperbolic conservation laws are typically given in terms of the \(L^1\)-norm in space. We also refer a reader to the recent works of Kröner and Rokyta on the error estimates of finite volume schemes for scalar convection–diffusion equation [17, 18].

When considering the multidimensional nonlinear system of hyperbolic conservation laws, to our best knowledge, the only result was done by Jovanović and Rohde [16]. The authors obtained a convergence rate of 1/2 in terms of the \(L^2\)-errors between the numerical solutions and a classical solution \(({\varvec{U}}\in C^1)\) under the assumptions of uniform boundedness of (numerical) solutions. Moreover, their result required the first-order derivatives to be bounded in the \(L^2\)- and \(L^{\infty }\)-norms.

In our recent work [22] we have proved the convergence of numerical solutions obtained by the Godunov method. In general we have obtained only weak* convergence to a generalised, dissipative measure-valued solution. If the limit is a weak entropy solution then the convergence is also strong. Moreover, if the Euler system admits a strong solution then the numerical solutions converge strongly to the strong solution as long as the latter exists.

Our aim is to extend the previous convergence analysis and show the error estimates for the Godunov method under the assumption that the Euler system (1.1)–(1.5) admits a strong solution. The main tool used in the present paper is the so-called relative energy functional originally introduced by Dafermos [5]. This technique has been largely used in the analysis of the weak–strong uniqueness and singular limit of the compressible fluid flows, see the monograph of Feireisl and Novotný [13], Březina and Feireisl [1], and Feireisl et al. [11, 12]. Recently, this technique has also been successfully applied to the convergence analysis of numerical solutions of compressible viscous fluids, see Feireisl et al. [7] and Mizerová and She [23]. Here we adapt the technique to the Euler system and estimate the corresponding relative energy, which yields the \(L^2\)-error estimates of density, momentum and entropy. We prove a convergence rate of 1/2 for the relative energy in the \(L^1\)-norm and a convergence rate of 1/4 for the \(L^2\)-errors. Furthermore, assuming that the total variation of the numerical solution is bounded, which is weaker than the assumptions in [16], we obtain the same convergence rate.

The rest of the paper is organized as follows. In Sect. 2 we introduce some preliminaries. More precisely, we recall the Godunov method and its consistency formulation proved in Lukáčová and Yuan [22]. We define the strong solution of the Euler system and the relative energy. Then, we prove the relative energy inequality in Sect. 3 and estimate its error in the \(L^1\)-norm. Finally, in Sect. 4 we present some numerical experiments to validate our theoretical results.

2 Preliminaries

In this section we introduce the preliminaries, including the formulation of the Godunov method, its consistency formulation, the definition of strong solution and the concept of the relative energy.

To begin with, we define the following notations for the later use

$$\begin{aligned}&\bullet \quad a \lesssim b \quad \text{ if } \quad a \le c\, b \quad \text{ with } \text{ a } \text{ positive } \text{ constant } \text{ c }, \\&\bullet \quad a \approx b \quad \text{ if } \quad a \lesssim b \quad \text{ and } \quad b \lesssim a. \end{aligned}$$

2.1 Godunov Method

The computational domain \(\Omega \) consists of rectangular meshes \({\overline{\Omega }} := \bigcup _{K} {\overline{K}}\). We denote the set of all mesh cells as \({\mathcal {T}}_h\) and the set of all interior faces of \({\mathcal {T}}_h\) as \(\Sigma _{\mathrm{int}}\). We consider the space of piecewise constant functions

$$\begin{aligned} {\mathcal {Q}}_h(\Omega ) = \{ v :v|_{{\overset{\circ }{K}} } = \text{ constant }, ~ \text{ for } \text{ all }~ K \in {\mathcal {T}}_h \} \end{aligned}$$
(2.1)

and define the projection operator

$$\begin{aligned} \Pi _h { :} L^1(\Omega ) \rightarrow {\mathcal {Q}}_h(\Omega ), \quad \Pi _h[ \phi ]_K = \frac{1}{|K|} \int _K \phi (x) ~dx, \end{aligned}$$
(2.2)

where |K| is the Lebesgue measure of K.

We are looking for \({\varvec{U}}_h(t) \in {\mathcal {Q}}_h(\Omega ;{\mathbb {R}}^{d+2}), \ t\in (0,T)\) satisfying the semi-discrete form of the finite volume method with the Godunov flux, i.e. the Godunov method,

$$\begin{aligned}&\int _{\Omega } {\varvec{\phi }} \frac{d }{d t} {\varvec{U}}_h ~ \,\mathrm{d} { x}- \sum _{\sigma \in \Sigma _{\texttt {int}}} \int _{\sigma } {\varvec{F}}({\varvec{U}}^{RP}_{\sigma }) \cdot {\varvec{n}}[\! [ {\varvec{\phi }} ] \! ] ~dS_{x} = 0, \end{aligned}$$
(2.3a)
$$\begin{aligned}&{\varvec{U}}_{h0} = \Pi _h[ {\varvec{U}}_0 ] . \end{aligned}$$
(2.3b)

Here \({\varvec{\phi }} \in {\mathcal {Q}}_h(\Omega ; {\mathbb {R}}^{d+2})\) is the test function, \({\varvec{U}}^{RP}_{\sigma }\) is the exact solution of the local Riemann problem along the interface \(\sigma \), \({\varvec{F}}\) is the flux function given (1.2), and the notation \([\! [ \cdot ] \! ]\) denotes the jump along the interface.

2.2 Consistency Formulation

This section introduces the necessary results of [22], i.e. the weak BV estimate and the consistency formulation, which will be needed in this paper. We refer to [22] for the details and the proofs. In the following we start with the following assumption.

Assumption 2.1

We assume that the solution to (2.3) satisfies

$$\begin{aligned} 0< {\underline{\varrho }} \le \varrho _h, \quad 0 < E_h \le {\overline{E}} \quad \text{ uniformly } \text{ for } h \rightarrow 0 \end{aligned}$$
(2.4)

for all \(t\in [0,T]\), where \({\underline{\varrho }}, {\overline{E}}\) are some positive constants.

From the above assumption we have the following estimates, see e.g. Feireisl et al. [10, 11].

Lemma 2.2

Under Assumption 2.1 there hold

$$\begin{aligned}&0< {\underline{\varrho }} \le \varrho _h \le {\overline{\varrho }}, \quad |{\varvec{u}}_h| \le {\overline{u}},\quad 0 < {\underline{p}} \le p_h \le {\overline{p}}, ~ \end{aligned}$$
(2.5)
$$\begin{aligned}&|{\varvec{m}}_h| \le {\overline{m}},\quad 0< {\underline{E}} \le E_h \le {\overline{E}}, \quad 0 < {\underline{\vartheta }} \le \vartheta _h \le {\overline{\vartheta }} \end{aligned}$$
(2.6)

uniformly for \(h \rightarrow 0, t\in [0,T]\) with positive constants \( {\overline{\varrho }},\, {\overline{u}},\, {\underline{p}},\, {\overline{p}},\, {\overline{m}},\, {\underline{E}},\,{\underline{\vartheta }},\,{\overline{\vartheta }}\) depending on \({\underline{\varrho }}, {\overline{E}}\), where \(\vartheta := \frac{p}{\varrho }\) is the absolute temperature.

Next, we report the consistency error of the Godunov method, see [22, Theorem 3.1].

Theorem 2.3

(Consistency formulation [22]) Let \((\varrho _h, {\varvec{m}}_h, \eta _h)\) be a numerical solution obtained by the Godunov method (2.3) on the time interval [0, T] satisfying Assumption 2.1. Then for any \(\tau \in (0,T)\) we have:

  • For all \(\phi \in W^{1,\infty }((0,T)\times \Omega )\) it holdsFootnote 1

    $$\begin{aligned} \left[ \int _{\Omega } \varrho _h\phi ~\,\mathrm{d} { x}\right] _{t = 0}^{t = \tau } = \int _{0}^{\tau } \int _{\Omega } \bigg ( \varrho _h\partial _t \phi + {\varvec{m}}_h\cdot \nabla _{ x}\phi \bigg ) \,\mathrm{d} { x}\,\mathrm{d} t + \int _{0}^{\tau } e_{\varrho ,h}(t,\phi ) \,\mathrm{d} t ; \end{aligned}$$
    (2.7)

  • For all \(\varvec{\phi }\in W^{1,\infty }((0,T)\times \Omega ; {\mathbb {R}}^d)\)

    $$\begin{aligned} \begin{aligned} \left[ \int _{\Omega } {\varvec{m}}_h\cdot \varvec{\phi }~\,\mathrm{d} { x}\right] _{t = 0}^{t = \tau } =&\int _{0}^{\tau } \int _{\Omega } \bigg ( {\varvec{m}}_h\cdot \partial _t \varvec{\phi }+ \frac{{\varvec{m}}_h\otimes {\varvec{m}}_h}{\varrho _h} : \nabla _{ x}\varvec{\phi }\\&+ p_h \mathrm{div}_{ x}\varvec{\phi }\bigg )\,\mathrm{d} { x}\,\mathrm{d} t + \int _{0}^{\tau } e_{{\varvec{m}},h}(t,\varvec{\phi }) \,\mathrm{d} t ; \end{aligned} \end{aligned}$$
    (2.8)

  • For all \(\phi \in W^{1,\infty }((0,T)\times \Omega ),\, \phi \ge 0\)

    $$\begin{aligned} \left[ \int _{\Omega } \eta _h \phi ~\,\mathrm{d} { x}\right] _{t = 0}^{t = \tau } \ge \int _{0}^{\tau } \int _{\Omega } \bigg ( \eta _h \partial _t \phi + {\varvec{q}}_h \cdot \nabla _{ x}\phi \bigg )\,\mathrm{d} { x}\,\mathrm{d} t + \int _{0}^{\tau } e_{\eta ,h}(t,\phi ) \,\mathrm{d} t ; \end{aligned}$$
    (2.9)
  • $$\begin{aligned} \int _{\Omega } E_h(\tau ) ~\,\mathrm{d} { x}= \int _{\Omega } E_{0,h} ~\,\mathrm{d} { x}\end{aligned}$$
    (2.10)

with bounded consistency errors \(e_{j,h} \ (j = \varrho , {\varvec{m}}, \eta )\) satisfying

$$\begin{aligned} \begin{aligned} \Vert e_{j,h}\Vert _{L^1(0,\tau )}&\lesssim h \Vert \phi \Vert _{W^{1,\infty }((0,T)\times \Omega )} \int _{0}^{\tau } \sum _{\sigma \in \Sigma _{\texttt {int}}} \int _{\sigma } \left| [\! [ {\varvec{U}}_h ] \! ]\right| ~ \mathrm {d}S_x\,\mathrm{d} t \\&\lesssim h^{1/2} \Vert \phi \Vert _{W^{1,\infty }((0,T)\times \Omega )} \left( \int _{0}^{\tau } \sum _{\sigma \in \Sigma _{\texttt {int}}} \int _{\sigma } \left| [\! [ {\varvec{U}}_h ] \! ]\right| ^2 ~ \mathrm {d}S_x\,\mathrm{d} t \right) ^{1/2}, \end{aligned} \end{aligned}$$
(2.11)

where \({\varvec{U}}_h = (\varrho _h, {\varvec{m}}_h, E_h)\) is the vector of conservative variables.

Theorem 2.3 presents a “weak” form of the Euler system satisfied by the numerical solutions modulus the consistency errors, which depend on the test function as well as on the jumps of numerical solutions on the interfaces. Note that the jumps can be controlled by the weak BV estimate stated in [22, equation (3.10)].

Lemma 2.4

(Weak BV estimate [22]) Under Assumption 2.1 for any \(\tau \in [0, T]\) it holds

$$\begin{aligned} \int _{0}^{\tau }\sum _{\sigma \in \Sigma _{\texttt {int}}} \int _{\sigma } |[\! [ {\varvec{U}}_{h} ] \! ] |^2 ~\mathrm {d}S_x\,\mathrm{d} t \lesssim 1. \end{aligned}$$
(2.12)

2.3 Strong Solution

In this section, the concept of strong solution of the Euler system (1.1)–(1.5) is introduced.

Definition 2.5

(Strong solution) Let \(\Omega \subset {\mathbb {R}}^d\) be a bounded domain with a boundary \(\partial \Omega \) of class \(C^1\). We say that a triple \([{{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}}]\) is the strong solution of the Euler system (1.1)–(1.5) if

$$\begin{aligned}&{{\widetilde{\varrho }}}\in W^{1,\infty }((0,T) \times \Omega ), ~ ~ {{{\widetilde{{\varvec{u}}}}}}\in W^{1,\infty }((0,T) \times \Omega ; {\mathbb {R}}^d), ~ ~ {{\widetilde{\eta }}}\in W^{1,\infty }((0,T) \times \Omega ), \\&{{\widetilde{\varrho }}}>0 \text{ and } \vartheta ({{\widetilde{\varrho }}}, {{\widetilde{\eta }}}) >0 \ \text{ for } \text{ any }\ (t,x) \in [0,T] \times \overline{\Omega } \end{aligned}$$

and the Eqs. (1.1)–(1.5) are satisfied almost everywhere.

Let us point out that we consider \(\varrho \) and \(\eta \) as the independent thermodynamical variables throughout the paper, meaning that all other thermodynamical variables are functions of \((\varrho , \eta )\). Accordingly, we write \({\widetilde{v}} = v({{\widetilde{\varrho }}}, {{\widetilde{\eta }}})\), \(v\in \{p, e, \vartheta , S\}\), for the strong solution. Moreover, we denote \({\widetilde{{\varvec{U}}}}:= ({{\widetilde{\varrho }}}, \widetilde{{\varvec{m}}}, {\widetilde{E}})^T\), where functions \(\widetilde{{\varvec{m}}}, {\widetilde{E}}\) are computed from the strong solution \(({{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}})\).

Since the domain is bounded and \(({{\widetilde{\varrho }}}, {{\widetilde{\vartheta }}})\) is continuous and positive, we have

$$\begin{aligned} 0< {\underline{\varrho }} \le {{\widetilde{\varrho }}},\quad 0 < {\underline{\vartheta }} \le {{\widetilde{\vartheta }}}. \end{aligned}$$
(2.13)

Remark 2.6

According to the definition of strong solution, we know that an entropy solution only containing finitely many rarefaction waves is also a strong solution.

We recall Gibbs’ relation

$$\begin{aligned} \vartheta \mathrm{d}S = \mathrm{d}e + p ~\mathrm{d}\left( 1/\varrho \right) . \end{aligned}$$
(2.14)

Consequently, for any strong solution \(({{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}})\) we obtain the following equations

$$\begin{aligned} \begin{aligned}&\partial _t {{\widetilde{\varrho }}}+ {{{\widetilde{{\varvec{u}}}}}}\cdot \nabla _{ x}{{\widetilde{\varrho }}}+ {{\widetilde{\varrho }}}\,\mathrm{div}_{ x}{{{\widetilde{{\varvec{u}}}}}}= 0, \\&\partial _t {{{\widetilde{{\varvec{u}}}}}}+ {{{\widetilde{{\varvec{u}}}}}}\cdot \nabla _{ x}{{{\widetilde{{\varvec{u}}}}}}+ \frac{1}{{{\widetilde{\varrho }}}}\, \nabla _{ x}{\widetilde{p}} = 0, \\&\partial _t {{\widetilde{\eta }}}+ {{{\widetilde{{\varvec{u}}}}}}\cdot \nabla _{ x}{{\widetilde{\eta }}}+ {{\widetilde{\eta }}}\, \mathrm{div}_{ x}{{{\widetilde{{\varvec{u}}}}}}= 0,\\&\partial _t {\widetilde{p}} + {{{\widetilde{{\varvec{u}}}}}}\cdot \nabla _{ x}{\widetilde{p}} + \gamma {\widetilde{p}} \, \mathrm{div}_{ x}{{{\widetilde{{\varvec{u}}}}}}= 0, \\&\partial _t {\widetilde{S}}+ {{{\widetilde{{\varvec{u}}}}}}\cdot \nabla _{ x}{\widetilde{S}}= 0, \\&\partial _t {{\widetilde{\vartheta }}}+ {{{\widetilde{{\varvec{u}}}}}}\cdot \nabla _{ x}{{\widetilde{\vartheta }}}+ \left( \partial _{{{\widetilde{\eta }}}} {\widetilde{p}} \right) _{{{\widetilde{\varrho }}}} \mathrm{div}_{ x}{{{\widetilde{{\varvec{u}}}}}}= 0, \end{aligned} \end{aligned}$$
(2.15)

which will be used in Sect. 3. See [11, 29] for more details on the existence of such a strong solution.

2.4 Relative Energy

Finally we conclude this section with the relative energy, which can measure the distance between the numerical solution and the strong solution of the Euler system. Moreover, we derive the relationship between the relative energy and the \(L^2\)-error of the numerical solution.

In the context of the compressible Euler system, the relative energy reads

$$\begin{aligned} {{\mathbb {E}}} \left( \varrho , {\varvec{m}}, \eta \mid {{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}}\right) = ~ \frac{1}{2} \varrho \left| \frac{{\varvec{m}}}{\varrho } - {{{\widetilde{{\varvec{u}}}}}}\right| ^2 + \varrho e - \frac{\partial (\varrho e)}{\partial \varrho }\Big |_{({{\widetilde{\varrho }}}, {{\widetilde{\eta }}})} (\varrho - {{\widetilde{\varrho }}}) - \frac{\partial (\varrho e)}{\partial \eta }\Big |_{({{\widetilde{\varrho }}}, {{\widetilde{\eta }}})} (\eta - {{\widetilde{\eta }}}) - {{\widetilde{\varrho }}}{{\widetilde{e}}}\nonumber \\ \end{aligned}$$
(2.16)

for \(\varrho >0\).

Lemma 2.7

Let \(({{\widetilde{\varrho }}},{{{\widetilde{{\varvec{u}}}}}},{{\widetilde{\eta }}})\) be the strong solution of the Euler system in the sense of Definition 2.5 and let \((\varrho _h,{\varvec{m}}_h,\eta _h)\) be the numerical solution of the Euler system obtained by (2.3) and satisfy Assumption 2.1. Then we have the following equivalence

$$\begin{aligned} {\mathbb {E}}\left( \varrho _h, {\varvec{m}}_h, \eta _h \mid {{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}}\right) \approx |{\varvec{m}}_h - \widetilde{{\varvec{m}}}|^2 + |\eta _h - {{\widetilde{\eta }}}|^2 + |\varrho _h- {{\widetilde{\varrho }}}|^2. \end{aligned}$$
(2.17)

Proof

The first step is to prove

$$\begin{aligned} {\mathbb {E}}\left( \varrho _h, {\varvec{m}}_h, \eta _h \mid {{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}}\right)&\approx |{\varvec{u}}_h - {{{\widetilde{{\varvec{u}}}}}}|^2 + |\eta _h - {{\widetilde{\eta }}}|^2 + |\varrho _h- {{\widetilde{\varrho }}}|^2. \end{aligned}$$
(2.18)

To this end we recall the definition of the relative energy (2.16). As \(\varrho _h\) is bounded from above and below by positive constants, we know that the first term on the right hand side of (2.16) can be estimated as

$$\begin{aligned} \frac{1}{2} \varrho _h\left| \frac{{\varvec{m}}_h}{\varrho _h} - {{{\widetilde{{\varvec{u}}}}}}\right| ^2 \approx |{\varvec{u}}_h- {{{\widetilde{{\varvec{u}}}}}}|^2. \end{aligned}$$

Note that the last four terms of (2.16) represent the second order remainder of the Taylor expansion of \(\varrho e\) around \(({{\widetilde{\varrho }}},{{\widetilde{\eta }}})\). Thus they can be estimated with the Hessian matrix \(\nabla _{(\varrho ,\eta )}^2 (\varrho e) \) evaluated at some point between \((\varrho _h, \eta _h)\) and \(({{\widetilde{\varrho }}}, {{\widetilde{\eta }}})\). Hence, the goal is to show the boundedness of the Hessian matrix \(\nabla _{(\varrho ,\eta )}^2 (\varrho e) \).

Taking the derivatives of \( \varrho e \) with respect to \( ( \varrho , \eta )\) we obtain

$$\begin{aligned} \partial _{\varrho } \left( \varrho e \right) =\left( 1+C_v \right) \vartheta - \frac{\eta \vartheta }{\varrho } , \quad \partial _{\eta } \left( \varrho e \right) = \vartheta . \end{aligned}$$
(2.19)

Further, applying the product rule and Gibbs’ relation (2.14) we derive

$$\begin{aligned} \mathrm{d}\vartheta =\frac{\vartheta }{C_v \varrho } \left( 1 - \frac{\eta }{\varrho }\right) \mathrm{d}\varrho + \frac{\vartheta }{C_v \varrho } \mathrm{d}\eta , \end{aligned}$$
(2.20)

and

$$\begin{aligned} \mathrm{d}\left( (1+C_v)\vartheta - \frac{\vartheta \eta }{\varrho } \right) = \left( (1+C_v) - \frac{\eta }{\varrho } \right) \mathrm{d}\vartheta - \frac{\vartheta }{\varrho } \mathrm{d}\eta + \frac{\vartheta \eta }{\varrho ^2} \mathrm{d}\varrho , \end{aligned}$$
(2.21)

which leads to

$$\begin{aligned} \nabla _{(\varrho ,\eta )}^2 (\varrho e) = \frac{\vartheta }{C_v \varrho } \begin{pmatrix} 1 &{} 1-\frac{\eta }{\varrho } \\ 1-\frac{\eta }{\varrho } &{} C_v + \left( 1-\frac{\eta }{\varrho } \right) ^2 \end{pmatrix}. \end{aligned}$$
(2.22)

Since \((\varrho _h, \eta _h)\) and \(({{\widetilde{\varrho }}}, {{\widetilde{\eta }}})\) are bounded, the Hessian matrix \(\nabla _{(\varrho ,\eta )}^2 (\varrho e) \) is symmetric positive definite. Moreover, its eigenvalues are bounded from below and above by positive constants, see “Appendix A”, which implies (2.18).

Next, we recall Assumption 2.1 and the uniform upper bound of \({\varvec{u}}_h\) due to Lemma 2.2 to conclude that

$$\begin{aligned}&|{\varvec{m}}_h - \widetilde{{\varvec{m}}}|^2 { \lesssim } \ |\varrho _h({\varvec{u}}_h- {{{\widetilde{{\varvec{u}}}}}})|^2 + |(\varrho _h - {{\widetilde{\varrho }}}) {{{\widetilde{{\varvec{u}}}}}}|^2 \lesssim |{\varvec{u}}_h- {{{\widetilde{{\varvec{u}}}}}}|^2 + |\varrho _h- {{\widetilde{\varrho }}}|^2,\\&|{\varvec{u}}_h- {{{\widetilde{{\varvec{u}}}}}}|^2 \lesssim |{{\widetilde{\varrho }}}({\varvec{u}}_h- {{{\widetilde{{\varvec{u}}}}}})|^2 \lesssim |{\varvec{m}}_h- \widetilde{{\varvec{m}}}|^2 + |{\varvec{u}}_h({{\widetilde{\varrho }}}- \varrho _h)|^2 \lesssim |{\varvec{m}}_h - \widetilde{{\varvec{m}}}|^2 + |\varrho _h- {{\widetilde{\varrho }}}|^2 . \end{aligned}$$

Substituting the above two inequalities into (2.18) we finish the proof. \(\square \)

Lemma 2.7 shows that the \(L^1\)-norm of \({\mathbb {E}}\left( \varrho _h, {\varvec{m}}_h, \eta _h \mid {{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}}\right) \) is equivalent to the square of the \(L^2\)-error of the numerical solution \((\varrho _h, {\varvec{m}}_h, \eta _h)\) as long as \(( \varrho _h, {\varvec{m}}_h, \eta _h)\) is obtained by a entropy stable scheme and satisfies Assumption 2.1 and \(({{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}})\) is the strong solution of the Euler system in the sense of Definition 2.5.

3 Error Estimates

Equipped with the consistency formulation of the Godunov method we are now ready to estimate the relative energy in the \(L^1\)-norm and the error of the numerical solution in the \(L^2\)-norm.

Theorem 3.1

(Error estimate) Let \(\Omega \subset {\mathbb {R}}^d\), \(d=1,2,3,\) be a bounded domain with a boundary \(\partial \Omega \) of class \(C^1\). Let \(( {{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}})\) be the strong solution of the complete Euler system (1.1) in the sense of Definition 2.5 with initial data satisfying

$$\begin{aligned} \Vert {\varvec{U}}_{h0} -{\varvec{U}}_0\Vert _{L^2(\Omega )} \lesssim h^{1/2}. \end{aligned}$$
(3.1)

Suppose that \(( \varrho _h, {\varvec{m}}_h, \eta _h )\) is the numerical solution obtained by the Godunov method (2.3). Let Assumption 2.1 hold. Then the following estimate of the relative energy holds for any \(\tau \in (0,T]\)

$$\begin{aligned} \int _{\Omega } {{\mathbb {E}}} \left( \varrho _h, {\varvec{m}}_h, \eta _h \mid {{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}}\right) (\tau , \cdot ) \ \,\mathrm{d} { x} \le D \ h^{1/2}, \end{aligned}$$
(3.2)

where D stands for a positive constant which depends only on \(\tau , |\Omega |\) and \(\Vert {\widetilde{{\varvec{U}}}}\Vert _{W^{1, \infty }((0,T) \times \Omega ; R^{d+2})} \).

Proof

The proof can be divided into two steps summarized as follows:

  • Taking suitable functions of the strong solution \(( {{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}})\) as test functions in the consistency formulation, we derive the relative energy inequality between \(( \varrho _h, {\varvec{m}}_h, \eta _h )\) and \(( {{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}})\);

  • Approximating the above inequality such that all terms on the right hand side can be bounded by the discretization parameter h or by the relative energy, we finally estimate the relative energy by Gronwall’s lemma.

Step 1. Rewriting the relative energy (2.16) into a more convenient form we obtain

$$\begin{aligned} \begin{aligned} { {\mathbb {E}}} \left( \varrho _h, {\varvec{m}}_h, \eta _h \mid {{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}}\right)&= ~ \frac{1}{2} \varrho _h \left| \frac{{\varvec{m}}_h}{\varrho _h} - {{{\widetilde{{\varvec{u}}}}}}\right| ^2 + \varrho _h e_h - \left( (1+C_v) {{\widetilde{\vartheta }}}- \frac{{{\widetilde{\vartheta }}}{{\widetilde{\eta }}}}{{{\widetilde{\varrho }}}} \right) (\varrho _h - {{\widetilde{\varrho }}})\\&\quad - {{\widetilde{\vartheta }}}(\eta _h - {{\widetilde{\eta }}}) - {{\widetilde{\varrho }}}{{\widetilde{e}}}\\&= \left[ \frac{1}{2} \frac{|{\varvec{m}}_h|^2}{\varrho _h} + \varrho _h e_h \right] + \varrho _h \left[ \frac{1}{2} |{{{\widetilde{{\varvec{u}}}}}}|^2 - (1+C_v){{\widetilde{\vartheta }}}+ \frac{{{\widetilde{\vartheta }}}{{\widetilde{\eta }}}}{{{\widetilde{\varrho }}}} \right] \\&\quad - {\varvec{m}}_h \cdot {{{\widetilde{{\varvec{u}}}}}}- \eta _h {{\widetilde{\vartheta }}}+ {{\widetilde{p}}}. \end{aligned} \end{aligned}$$
(3.3)

Then, we take \(\frac{1}{2} |{{{\widetilde{{\varvec{u}}}}}}|^2 - (1+ C_v) {{\widetilde{\vartheta }}}+ \frac{{{\widetilde{\vartheta }}}{{\widetilde{\eta }}}}{{{\widetilde{\varrho }}}}\) as the test function in consistency formulation of the density Eq. (2.7) to derive

$$\begin{aligned}&\left[ \int _{\Omega } \varrho _h\left( \frac{1}{2} |{{{\widetilde{{\varvec{u}}}}}}|^2 - (1+ C_v) {{\widetilde{\vartheta }}}+ \frac{{{\widetilde{\vartheta }}}{{\widetilde{\eta }}}}{{{\widetilde{\varrho }}}}\right) ~\,\mathrm{d} { x}\right] _{t = 0}^{t = \tau } = \int _{0}^{\tau } \int _{\Omega } \Bigg ( \varrho _h\partial _t \left( \frac{1}{2} |{{{\widetilde{{\varvec{u}}}}}}|^2 - (1+ C_v) {{\widetilde{\vartheta }}}+ \frac{{{\widetilde{\vartheta }}}{{\widetilde{\eta }}}}{{{\widetilde{\varrho }}}}\right) \\&\quad + {\varvec{m}}_h\cdot \nabla _{ x}\left( \frac{1}{2} |{{{\widetilde{{\varvec{u}}}}}}|^2 - (1+ C_v) {{\widetilde{\vartheta }}}+ \frac{{{\widetilde{\vartheta }}}{{\widetilde{\eta }}}}{{{\widetilde{\varrho }}}}\right) \Bigg ) \,\mathrm{d} { x}\,\mathrm{d} t + \int _{0}^{\tau } e_{\varrho ,h}\left( t, { \frac{1}{2} |{{{\widetilde{{\varvec{u}}}}}}|^2 - (1+ C_v) {{\widetilde{\vartheta }}}+ \frac{{{\widetilde{\vartheta }}}{{\widetilde{\eta }}}}{{{\widetilde{\varrho }}}} } \right) \,\mathrm{d} t . \end{aligned}$$

Analogously, setting \({{{\widetilde{{\varvec{u}}}}}}\) and \({{\widetilde{\vartheta }}}\) respectively as the test functions in consistency formulations of the momentum Eq. (2.8) and entropy inequality (2.9) we obtain

$$\begin{aligned} \left[ \int _{\Omega } {\varvec{m}}_h\cdot {{{\widetilde{{\varvec{u}}}}}}~\,\mathrm{d} { x}\right] _{t = 0}^{t = \tau }&=~\int _{0}^{\tau } \int _{\Omega } \left( {\varvec{m}}_h\cdot \partial _t {{{\widetilde{{\varvec{u}}}}}}+ \frac{{\varvec{m}}_h\otimes {\varvec{m}}_h}{\varrho _h} : \nabla _{ x}{{{\widetilde{{\varvec{u}}}}}}+ p_h \mathrm{div}_{ x}{{{\widetilde{{\varvec{u}}}}}} \right) \ \,\mathrm{d} { x} \,\mathrm{d} t \\&\quad + \int _{0}^{\tau } e_{{\varvec{m}},h}(t, { {{{\widetilde{{\varvec{u}}}}}}}) \,\mathrm{d} t , \end{aligned}$$

and

$$\begin{aligned} \left[ \int _{\Omega } \eta _h {{\widetilde{\vartheta }}}~\,\mathrm{d} { x}\right] _{t = 0}^{t = \tau } \ge ~\int _{0}^{\tau } \int _{\Omega } \left( \eta _h \partial _t {{\widetilde{\vartheta }}}+ \eta _h \frac{{\varvec{m}}_h}{\varrho _h} \cdot \nabla _{ x}{{\widetilde{\vartheta }}} \right) \ \,\mathrm{d} { x}\,\mathrm{d} t + \int _{0}^{\tau } e_{\eta ,h}(t,{ {{\widetilde{\vartheta }}}}) \,\mathrm{d} t . \end{aligned}$$

Thus, substituting the above three formulae together with the energy equality (2.10) into the integral of the relative energy (3.3) we have

$$\begin{aligned} \begin{aligned}&\left[ \int _{\Omega } {\mathbb {E}}\left( \varrho _h, {\varvec{m}}_h, \eta _h \mid {{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}}\right) (t, \cdot ) \ \,\mathrm{d} { x} \right] _{t = 0}^{t = \tau } \\&\quad \le - \int _0^\tau \int _{\Omega } \frac{(\varrho _h{{{\widetilde{{\varvec{u}}}}}}- {\varvec{m}}_h) \otimes (\varrho _h{{{\widetilde{{\varvec{u}}}}}}- {\varvec{m}}_h) }{\varrho _h} : \nabla _{ x}{{{\widetilde{{\varvec{u}}}}}} \ \,\mathrm{d} { x} \,\mathrm{d} t \\&\qquad + \int _0^\tau \int _{\Omega } \big ( ( {{\widetilde{p}}}- p_h ) \mathrm{div}_{ x}{{{\widetilde{{\varvec{u}}}}}}+ (\partial _t {{\widetilde{p}}}+ {{{\widetilde{{\varvec{u}}}}}}\cdot \nabla _{ x}{{\widetilde{p}}}) \big ) \ \,\mathrm{d} { x} \,\mathrm{d} t \\&\qquad + \int _0^\tau \int _{\Omega } (\varrho _h{{{\widetilde{{\varvec{u}}}}}}- {\varvec{m}}_h) \cdot \left[ \partial _t {{{\widetilde{{\varvec{u}}}}}}+ {{{\widetilde{{\varvec{u}}}}}}\cdot \nabla _{ x}{{{\widetilde{{\varvec{u}}}}}}+ \frac{1}{{{\widetilde{\varrho }}}} \nabla _{ x}{{\widetilde{p}}}\right] \ \,\mathrm{d} { x} \,\mathrm{d} t \\&\qquad - \int _0^\tau \int _{\Omega } \left( \varrho _h\partial _t \left( (1+ C_v){{\widetilde{\vartheta }}}- \frac{{{\widetilde{\vartheta }}}{{\widetilde{\eta }}}}{{{\widetilde{\varrho }}}} \right) + {\varvec{m}}_h\cdot \nabla _{ x}\left( (1+ C_v){{\widetilde{\vartheta }}}- \frac{{{\widetilde{\vartheta }}}{{\widetilde{\eta }}}}{{{\widetilde{\varrho }}}} \right) \right) \ \,\mathrm{d} { x} \,\mathrm{d} t \\&\qquad - \int _0^\tau \int _{\Omega } \left( \eta _h \partial _t {{\widetilde{\vartheta }}}+ \eta _h \frac{{\varvec{m}}_h}{\varrho _h} \cdot \nabla _{ x}{{\widetilde{\vartheta }}}+ (\varrho _h{{{\widetilde{{\varvec{u}}}}}}- {\varvec{m}}_h) \frac{1}{{{\widetilde{\varrho }}}} \nabla _{ x}{{\widetilde{p}}} \right) \ \,\mathrm{d} { x} \,\mathrm{d} t \\&\qquad + \int _{0}^{\tau } e_{h}(t, {\widetilde{{\varvec{U}}}}) \,\mathrm{d} t , \end{aligned} \end{aligned}$$
(3.4)

where we have used the following identities

$$\begin{aligned} {{{\widetilde{{\varvec{u}}}}}}\otimes {{{\widetilde{{\varvec{u}}}}}}: \nabla _{ x}{{{\widetilde{{\varvec{u}}}}}}= {{{\widetilde{{\varvec{u}}}}}}\cdot ({{{\widetilde{{\varvec{u}}}}}}\cdot \nabla _{ x}{{{\widetilde{{\varvec{u}}}}}}), \qquad \int _{\Omega } {{{\widetilde{{\varvec{u}}}}}}\cdot \nabla _{ x}{{\widetilde{p}}} \ \,\mathrm{d} { x} = - \int _{\Omega } {{\widetilde{p}}}\,\mathrm{div}_{ x}{{{\widetilde{{\varvec{u}}}}}} \ \,\mathrm{d} { x} , \\ \frac{(\varrho _h{{{\widetilde{{\varvec{u}}}}}}- {\varvec{m}}_h) \otimes (\varrho _h{{{\widetilde{{\varvec{u}}}}}}- {\varvec{m}}_h) }{\varrho _h} : \nabla _{ x}{{{\widetilde{{\varvec{u}}}}}}= ~ \varrho _h{{{\widetilde{{\varvec{u}}}}}}\otimes {{{\widetilde{{\varvec{u}}}}}}: \nabla _{ x}{{{\widetilde{{\varvec{u}}}}}}+ \frac{{\varvec{m}}_h\otimes {\varvec{m}}_h}{\varrho _h} : \nabla _{ x}{{{\widetilde{{\varvec{u}}}}}}\\ - {\varvec{m}}_h\cdot ({{{\widetilde{{\varvec{u}}}}}}\cdot \nabla _{ x}) {{{\widetilde{{\varvec{u}}}}}}- {{{\widetilde{{\varvec{u}}}}}}\cdot \left( {\varvec{m}}_h\cdot \nabla _{ x}\right) {{{\widetilde{{\varvec{u}}}}}}\end{aligned}$$

and the notation

$$\begin{aligned} e_{h}(t, {\widetilde{{\varvec{U}}}}) := e_{\varrho ,h}\left( t, \frac{1}{2} |{{{\widetilde{{\varvec{u}}}}}}|^2 - (1+ C_v) {{\widetilde{\vartheta }}}+ \frac{{{\widetilde{\vartheta }}}{{\widetilde{\eta }}}}{{{\widetilde{\varrho }}}} \right) - e_{{\varvec{m}},h}(t, {{{\widetilde{{\varvec{u}}}}}}) - e_{\eta ,h}(t, {{\widetilde{\vartheta }}}). \end{aligned}$$

Further, employing the relations (2.20) and (2.21) we obtain after lengthy but straightforward calculations from (3.4)

$$\begin{aligned}&\left[ \int _{\Omega } {{\mathbb {E}}} \left( \varrho _h, {\varvec{m}}_h, \eta _h \mid {{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}}\right) (t, \cdot ) \ \,\mathrm{d} { x} \right] _{t = 0}^{t = \tau } \nonumber \\&\quad \le - \int _0^\tau \int _{\Omega } \frac{(\varrho _h{{{\widetilde{{\varvec{u}}}}}}- {\varvec{m}}_h) \otimes (\varrho _h{{{\widetilde{{\varvec{u}}}}}}- {\varvec{m}}_h) }{\varrho _h} : \nabla _{ x}{{{\widetilde{{\varvec{u}}}}}} \ \,\mathrm{d} { x} \,\mathrm{d} t \nonumber \\&\qquad - \int _0^\tau \int _{\Omega } \Big [ p_h - {{\widetilde{p}}}- \partial _{{{\widetilde{\varrho }}}} {{\widetilde{p}}}(\varrho _h- {{\widetilde{\varrho }}})- \partial _{{{\widetilde{\eta }}}} {{\widetilde{p}}}(\eta _h - {{\widetilde{\eta }}}) \Big ] \mathrm{div}_{ x}{{{\widetilde{{\varvec{u}}}}}} \ \,\mathrm{d} { x} \,\mathrm{d} t \nonumber \\&\qquad + \int _0^\tau \int _{\Omega } (\varrho _h{{{\widetilde{{\varvec{u}}}}}}- {\varvec{m}}_h) \cdot \left[ \partial _t {{{\widetilde{{\varvec{u}}}}}}+ {{{\widetilde{{\varvec{u}}}}}}\cdot \nabla _{ x}{{{\widetilde{{\varvec{u}}}}}}+ \frac{1}{{{\widetilde{\varrho }}}} \nabla _{ x}{{\widetilde{p}}}\right] \ \,\mathrm{d} { x} \,\mathrm{d} t \nonumber \\&\qquad + \int _0^\tau \int _{\Omega } \left[ \left( 1- \frac{\varrho _h}{{{\widetilde{\varrho }}}} \right) \partial _{{{\widetilde{\varrho }}}} {{\widetilde{p}}}- \left( \eta _h- \frac{\varrho _h}{{{\widetilde{\varrho }}}} {{\widetilde{\eta }}}\right) \partial _{{{\widetilde{\varrho }}}} {{\widetilde{\vartheta }}}\right] \left[ \partial _t {{\widetilde{\varrho }}}+ {{{\widetilde{{\varvec{u}}}}}}\cdot \nabla _{ x}{{\widetilde{\varrho }}}+ {{\widetilde{\varrho }}}\, \mathrm{div}_{ x}{{{\widetilde{{\varvec{u}}}}}}\right] \ \,\mathrm{d} { x} \,\mathrm{d} t \nonumber \\&\qquad + \int _0^\tau \int _{\Omega } \left[ \left( 1- \frac{\varrho _h}{{{\widetilde{\varrho }}}} \right) \partial _{{{\widetilde{\eta }}}} {{\widetilde{p}}}- \left( \eta _h- \frac{\varrho _h}{{{\widetilde{\varrho }}}} {{\widetilde{\eta }}}\right) \partial _{{{\widetilde{\eta }}}} {{\widetilde{\vartheta }}}\right] \left[ \partial _t {{\widetilde{\eta }}}+ {{{\widetilde{{\varvec{u}}}}}}\cdot \nabla _{ x}{{\widetilde{\eta }}}+ {{\widetilde{\eta }}}\, \mathrm{div}_{ x}{{{\widetilde{{\varvec{u}}}}}}\right] \ \,\mathrm{d} { x} \,\mathrm{d} t \nonumber \\&\qquad - \int _0^\tau \int _{\Omega } \left( \eta _h - \frac{\varrho _h}{{{\widetilde{\varrho }}}} {{\widetilde{\eta }}}\right) \left( \frac{{\varvec{m}}_h}{\varrho _h} - {{{\widetilde{{\varvec{u}}}}}}\right) \cdot \nabla _{ x}{{\widetilde{\vartheta }}} \ \,\mathrm{d} { x} \,\mathrm{d} t + \int _{0}^{\tau } e_{h}(t, {\widetilde{{\varvec{U}}}}) \,\mathrm{d} t \end{aligned}$$
(3.5)

where we have denoted \(\partial _{{{\widetilde{\varrho }}}} {{\widetilde{p}}}:= \frac{\partial p}{\partial \varrho }({{\widetilde{\varrho }}},{{\widetilde{\eta }}})\) and the definitions of \(\partial _{{{\widetilde{\eta }}}} {{\widetilde{p}}}\), \(\partial _{{{\widetilde{\eta }}}} {{\widetilde{\vartheta }}}\) and \(\partial _{{{\widetilde{\eta }}}} {{\widetilde{\vartheta }}}\) are analogous. Then applying the equalities stated in (2.15)–(3.5) we have

$$\begin{aligned} \begin{aligned}&\left[ \int _{\Omega } {\mathbb {E}}\left( \varrho _h, {\varvec{m}}_h, \eta _h \mid {{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}}\right) (t, \cdot ) \ \,\mathrm{d} { x} \right] _{t = 0}^{t = \tau } \\&\quad \le - \int _0^\tau \int _{\Omega } \frac{(\varrho _h{{{\widetilde{{\varvec{u}}}}}}- {\varvec{m}}_h) \otimes (\varrho _h{{{\widetilde{{\varvec{u}}}}}}- {\varvec{m}}_h) }{\varrho _h} : \nabla _{ x}{{{\widetilde{{\varvec{u}}}}}} \ \,\mathrm{d} { x} \,\mathrm{d} t \\&\qquad - \int _0^\tau \int _{\Omega } \Big [ p_h - {{\widetilde{p}}}- \partial _{{{\widetilde{\varrho }}}} {{\widetilde{p}}}(\varrho _h- {{\widetilde{\varrho }}})- \partial _{{{\widetilde{\eta }}}} {{\widetilde{p}}}(\eta _h - {{\widetilde{\eta }}}) \Big ] \mathrm{div}_{ x}{{{\widetilde{{\varvec{u}}}}}} \ \,\mathrm{d} { x} \,\mathrm{d} t \\&\qquad - \int _0^\tau \int _{\Omega } \left( \eta _h - \frac{\varrho _h}{{{\widetilde{\varrho }}}} {{\widetilde{\eta }}}\right) \left( \frac{{\varvec{m}}_h}{\varrho _h} - {{{\widetilde{{\varvec{u}}}}}}\right) \cdot \nabla _{ x}{{\widetilde{\vartheta }}} \ \,\mathrm{d} { x} \,\mathrm{d} t + \int _{0}^{\tau } e_{h}(t, {\widetilde{{\varvec{U}}}}) \,\mathrm{d} t \\ \end{aligned} \end{aligned}$$
(3.6)

Step 2. We begin with the following observation owing to the uniform bounds on \({{\widetilde{\varrho }}}\), \({{\widetilde{\vartheta }}}\) and \({{\widetilde{\eta }}}\), as well as (2.18)

$$\begin{aligned}&\left| \left( \eta _h - \frac{\varrho _h}{{{\widetilde{\varrho }}}} {{\widetilde{\eta }}}\right) \cdot \left( \frac{{\varvec{m}}_h}{\varrho _h} - {{{\widetilde{{\varvec{u}}}}}}\right) \right| \lesssim \left| \ \eta _h - \frac{\varrho _h}{{{\widetilde{\varrho }}}} {{\widetilde{\eta }}}\right| ^2 + \left| \frac{{\varvec{m}}_h}{\varrho _h} - {{{\widetilde{{\varvec{u}}}}}}\right| ^2 \\&\quad \lesssim \left| \ \eta _h -{{\widetilde{\eta }}}\right| ^2+ \left| {{\widetilde{\eta }}}- \frac{\varrho _h}{{{\widetilde{\varrho }}}} {{\widetilde{\eta }}}\right| ^2 + \left| \frac{{\varvec{m}}_h}{\varrho _h} - {{{\widetilde{{\varvec{u}}}}}}\right| ^2 = \left| \ \eta _h -{{\widetilde{\eta }}}\right| ^2+ \left| {{\widetilde{\varrho }}}- \varrho _h\right| ^2 \left( \frac{{{\widetilde{\eta }}}}{{{\widetilde{\varrho }}}} \right) ^2 + \left| \frac{{\varvec{m}}_h}{\varrho _h} - {{{\widetilde{{\varvec{u}}}}}}\right| ^2 \\&\quad \lesssim |\eta _h - {{\widetilde{\eta }}}|^2 + |\varrho _h- {{\widetilde{\varrho }}}|^2+ \left| \frac{{\varvec{m}}_h}{\varrho _h} - {{{\widetilde{{\varvec{u}}}}}}\right| ^2 \lesssim {\mathbb {E}}\left( \varrho _h, {\varvec{m}}_h, \eta _h \mid {{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}}\right) . \end{aligned}$$

Recall that the consistency errors are bounded by \(h^{1/2}\), cf. (2.11) and (2.12). Since the first and the second term on the right hand side of inequality (3.6) are bounded by the relative energy, we obtain

$$\begin{aligned}&\left[ \int _{\Omega } {\mathbb {E}}\left( \varrho _h, {\varvec{m}}_h, \eta _h \mid {{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}}\right) (t, \cdot ) \ \,\mathrm{d} { x} \right] _{t = 0}^{t = \tau } \le c\left( |\Omega |, \Vert {\widetilde{{\varvec{U}}}}\Vert _{W^{1, \infty }((0,T) \times \Omega ; R^{d+2})} \right) h^{1/2} \nonumber \\&\quad + c\left( \Omega , \Vert {\widetilde{{\varvec{U}}}}\Vert _{W^{1, \infty }((0,T) \times \Omega ; R^{d+2})} \right) \int _0^\tau \int _{\Omega } {\mathbb {E}}\left( \varrho _h, {\varvec{m}}_h, \eta _h \mid {{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}}\right) (t, \cdot ) \ \,\mathrm{d} { x} \,\mathrm{d} t . \end{aligned}$$
(3.7)

Recalling the assumption on initial data (3.1) and Lemma 2.7 we obtain the first order estimate of the discrete initial data, i.e.

$$\begin{aligned} \int _{\Omega } {\mathbb {E}}\left( \varrho _{h0}, {\varvec{m}}_{h0}, \eta _{h0} \mid \varrho _0, {\varvec{m}}_0/\varrho _0, \eta _0 \right) \ \,\mathrm{d} { x} \lesssim h.\end{aligned}$$

Now applying Gronwall’s lemma concludes the proof, i.e.

$$\begin{aligned}&\int _{\Omega } {\mathbb {E}}\left( \varrho _h, {\varvec{m}}_h, \eta _h \mid {{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}}\right) (\tau , \cdot ) \ \,\mathrm{d} { x} \le \exp \left( \tau \; c\left( |\Omega |, \Vert {\widetilde{{\varvec{U}}}}\Vert _{W^{1, \infty }((0,T) \times \Omega ; R^d)} \right) \right) \cdot \\&\quad \quad \left( c\left( |\Omega |, \Vert {\widetilde{{\varvec{U}}}}\Vert _{W^{1, \infty }((0,T) \times \Omega ; R^d)} \right) h^{1/2} + \int _{\Omega } {\mathbb {E}}\left( \varrho _{h0}, {\varvec{m}}_{h0}, \eta _{h0} \mid \varrho _0, {\varvec{m}}_0/\varrho _0, \eta _0 \right) \ \,\mathrm{d} { x} \right) \\&\qquad \quad \le { D \ h^{1/2}, } \end{aligned}$$

where D is the constant depending on \(\tau , |\Omega |, \Vert {\widetilde{{\varvec{U}}}}\Vert _{W^{1, \infty }((0,T) \times \Omega ; R^d)}\). \(\square \)

Combining Lemma 2.7 we directly obtained the following a priori error estimates in the \(L^2\) norm.

Proposition 3.2

Under the same condition as Theorem 3.1 it holds for any \(\tau \in (0,T]\)

$$\begin{aligned} \Vert \varrho _h - {{\widetilde{\varrho }}}\Vert _{L^2(\Omega )} \lesssim h^{1/4}, \quad \Vert {\varvec{m}}_h - \widetilde{{\varvec{m}}}\Vert _{L^2(\Omega )} \lesssim h^{1/4}, \quad \Vert \eta _h - {{\widetilde{\eta }}}\Vert _{L^2(\Omega )} \lesssim h^{1/4}. \end{aligned}$$
(3.8)

In what follows we obtain the first order convergence rate in terms of the relative energy under an additional assumption on uniform boundedness of the total variation of numerical solutions.

Theorem 3.3

In addition to the assumptions of Theorem 3.1, we assume that

$$\begin{aligned} \sum _{\sigma \in \Sigma _{\texttt {int}}} \int _{\sigma } |[\! [ {\varvec{U}}_h ] \! ]|_{\sigma } \mathrm {d}S_x\lesssim 1. \end{aligned}$$
(3.9)

Then it holds

$$\begin{aligned} \int _{\Omega } {\mathbb {E}}\left( \varrho _h, {\varvec{m}}_h, \eta _h \mid {{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}}\right) (\tau , \cdot ) \ \,\mathrm{d} { x} \lesssim D \ h, \end{aligned}$$

where D stands for a positive constant which depends only on \(\tau , |\Omega |\) and \(\Vert {\widetilde{{\varvec{U}}}}\Vert _{W^{1, \infty }((0,T) \times \Omega ; R^{d+2})} \).

Proof

With the help of (3.9) and (2.11) the consistency error \(e_{h}(t, {\widetilde{{\varvec{U}}}})\) in (3.6) can be improved as follows

$$\begin{aligned} \Vert e_{h}(t, {\widetilde{{\varvec{U}}}})\Vert _{L^1(0,T)} \lesssim \Vert {\widetilde{{\varvec{U}}}}\Vert _{W^{1,\infty }((0,T)\times {\Omega })} \ h, \end{aligned}$$

which concludes the proof. \(\square \)

Remark 3.4

Here we point out that assumption (3.9) is slightly weaker than the assumption used in the work of Jovanović and Rohde [16]

$$\begin{aligned} \sum _{\sigma \in \Sigma _{\texttt {int}}} \int _{\sigma } \frac{| [\! [ {\varvec{U}}_h ] \! ]| ^2}{h} \mathrm {d}S_x\lesssim 1. \end{aligned}$$
(3.10)

Moreover, for the case of \(d=1\) the assumption (3.9) is exactly the TVB condition, which is a known property for the Godunov method.

Remark 3.5

Let us consider piecewise constant initial data which generate only finitely many rarefaction waves. It is obvious that such kind of initial data fulfills the condition \(\Vert {\varvec{U}}_{h0} - {\varvec{U}}_0 \Vert _{L^2(\Omega )} \lesssim h^{1/2} \) assumed in Theorem 3.1. Moreover, we can expect (3.9) or (3.10) to hold, which consequently implies Theorem 3.3. Thus, in this case it holds \(\Vert {\varvec{U}}_h - {\widetilde{{\varvec{U}}}}\Vert _{L^2(\Omega )} \lesssim h^{1/2} \) for any \(\tau \in (0,T)\).

4 Numerical Experiments

In this section we simulate several one- and two-dimensional Riemann problems. The examples only containing rarefaction waves are used to validate our theoretical results, meanwhile, the other examples containing contact waves or shock waves or both are also tested for comparisons and future interests. We point out that in our simulations there is no projection error of initial data due to these simple Riemann problems and good uniform meshes.

In addition to the Godunov method, we also test the convergence rates of the viscosity finite volume (VFV) method in order to generalize our analysis and for comparison. The VFV method was originally introduced and studied by Feireisl et al. [9]. It can be described as follows

$$\begin{aligned} D_t \varrho _h + \mathrm{div}_h^{\mathrm{{up}}} (\varrho _h, {\varvec{u}}_h)&= 0, \\ D_t {\varvec{m}}_h + \mathrm{div}_h^{\mathrm{{up}}} ({\varvec{m}}_h, {\varvec{u}}_h) + \nabla _h p_h&= { h^{\alpha }} \Delta _h {\varvec{u}}_h,\\ D_t E_h + \mathrm{div}_h^{\mathrm{{up}}} (E_h, {\varvec{u}}_h) +p_h \mathrm{div}_h{\varvec{u}}_h +{\varvec{u}}_h \cdot \nabla _h p_h&= \frac{1}{2} { h^{\alpha }} \Delta _h (|{\varvec{u}}_h|^2). \end{aligned}$$

The discrete operators are defined by

$$\begin{aligned}&\big ( \mathrm{div}_h^{\mathrm{{up}}}(r_h, {\varvec{v}}_h) \big )_K = \sum _{\sigma \in \partial K} \frac{|\sigma |}{|K|} F_h(r_h,{\varvec{v}}_h), \quad (\mathrm{div}_h{\varvec{v}}_h)_K = \sum _{\sigma \in \partial K} \frac{|\sigma |}{|K|} \left\{ \!\left\{ {\varvec{v}}_h\right\} \!\right\} \cdot {\varvec{n}},\\&(\nabla _hr_h)_K = \sum _{\sigma \in \partial K} \frac{|\sigma |}{|K|} \left\{ \!\left\{ r_h\right\} \!\right\} {\varvec{n}}, \quad \Delta _h r_h = \sum _{\sigma \in \partial K} \frac{|\sigma |}{|K|} \frac{[\! [ r_h ] \! ] }{h} , \\&F_h(r_h,{\varvec{v}}_h) = Up[r_h, {\varvec{v}}_h] - h^{\beta } ~ [\! [ r_h ] \! ], \quad 0< \beta <1, \\&Up[r_h, {\varvec{v}}_h] = \left\{ \!\left\{ r_h\right\} \!\right\} \, \left\{ \!\left\{ {\varvec{v}}_h\right\} \!\right\} \cdot {\varvec{n}}- \frac{1}{2} \left| \left\{ \!\left\{ {\varvec{v}}_h\right\} \!\right\} \cdot {\varvec{n}}\right| ~ [\! [ r_h ] \! ], \end{aligned}$$

In the following simulations, we use the forward Euler time discretization for both methods, and take \(\mathrm CFL=0.9\) for the Godunov method, while \(\mathrm CFL=0.3\), \(\alpha =1.8\), \(\beta =0.2\) for the VFV method. Unless otherwise specified, we take \(\gamma =1.4\); the errors of \((\varrho ,{\varvec{m}},\eta ), {\mathbb {E}}\) mean the \(L^2\)-error of \((\varrho _h,{\varvec{m}}_h,\eta _h)\) and the \(L^1\)-norm of the relative energy \({\mathbb {E}}( \varrho _h, {\varvec{m}}_h, \eta _h | {{\widetilde{\varrho }}}, {{{\widetilde{{\varvec{u}}}}}}, {{\widetilde{\eta }}}) \); the convergence rates of \((\varrho ,{\varvec{m}},\eta ), {\mathbb {E}}\) mean the convergence rates of the errors of \((\varrho ,{\varvec{m}},\eta ), {\mathbb {E}}\). In addition, the error diagrams of \(\varrho , {\varvec{m}}, \eta \) and \({\mathbb {E}}\) (also denoted as RE in the plots) will be drawn with symbols “\(\vartriangle \)”, “\(+\)”, “\(\circ \)” and “\(\square \)”, respectively. Moreover, the solid line without a marker represents the reference slope of \(h^{1/2}\).

4.1 One Dimensional Experiments

We start with one dimensional Riemann problems in the computational domain \(\Omega =[0,1]\). Here, the strong solution \({\widetilde{{\varvec{U}}}}\) in the relative energy is taken as the reference (exact) solution computed on the uniform mesh with 20,480 cells.

Example 4.1

(1D single wave) This example is used to measure the convergence rate of three different types of waves—a single contact (C) wave, a single rarefaction (R) wave and a single shock (S) wave.

Given the initial data in Table 1, we compute the contact, rarefaction and shock wave till \(T = 0.2, 0.2\) and 0.25, respectively. Figure 1 (resp. Fig. 2) shows the density \(\varrho \) (resp. the entropy \(\eta \)) obtained on different meshes with \(n =1/h =32, 64, \dots ,1024\) cells. Moreover, we present in Fig. 3 the errors of \((\varrho , {\varvec{m}}, \eta ), {\mathbb {E}}\), see the details in Tables 2 and 3.

The numerical results indicate that :

  • The Godunov method and the VFV method have similar convergence rates. Moreover, the convergence rates of the VFV method are slightly better than those of the Godunov method;

  • For the single rarefaction wave the convergence rate of \((\varrho ,{\varvec{m}},\eta )\) (resp. \({\mathbb {E}}\)) is slightly greater than 1/2 (resp. 1), which is consistent to our theoretical results;

  • For the single contact wave the convergence rate of \((\varrho ,{\varvec{m}},\eta )\) (resp. \({\mathbb {E}}\)) is around 1/4 (resp. 1/2);

  • For the single shock wave the convergence rate of \((\varrho ,{\varvec{m}},\eta )\) (resp. \({\mathbb {E}}\)) is around 1/2 (resp. 1).

In Fig. 2 we see the approximation of the entropy of the order \(10^{-2}\), which is the same as the order of the density error in Fig. 1.

Table 1 Initial data of 1D single wave
Full size table
Fig. 1
figure 1

Example 4.1: density \(\varrho \) obtained by the Godunov method (top) and the VFV method (bottom)

Full size image

Remark 4.2

Here we compare the above observation with the result of Tadmor and Tang [25] for the rarefaction wave and the shock wave.

  • Directly applying the pointwise error estimate for scalar equation in [25], i.e.

    $$\begin{aligned} |(u^{\varepsilon } - u)(x,t)| \approx \text{ dist }(x,R(t))^{-1}\varepsilon \log ^2 \varepsilon \end{aligned}$$

    with rarefaction set R(t), we obtain that the \(L^2\)-error is bounded by \(\varepsilon ^{1/2} \log ^2 \varepsilon \). Setting the vanishing viscosity coefficient \(\varepsilon \approx h\) means that our numerical analysis gives a better upper bound for the convergence rate in the case of a single rarefaction wave.

  • In the case of a shock applying the pointwise error estimate for scalar equation in [25], i.e.

    $$\begin{aligned} |(u^{\varepsilon } - u)(x,t)| \approx \text{ dist }(x,S(t))^{-1} \varepsilon , \end{aligned}$$

    where S(t) is the streamline of shock discontinuities, we obtain that the \(L^2\)-convergence rate is 1/2, which is consistent with our observations.

Fig. 2
figure 2

Example 4.1: entropy \(\eta \) obtained by the Godunov method (top) and the VFV method (bottom)

Full size image
Fig. 3
figure 3

Example 4.1: errors obtained with different stepsizes \(h =1/32, \dots , 1/1024\). The black solid lines without any marker denote the reference slope of \(h^{1/2}\)

Full size image
Table 2 Example 4.1: Errors and convergence rates of \(\varrho ,\eta , {\mathbb {E}}\) obtained by the Godunov method
Full size table
Table 3 Example 4.1: errors and convergence rates of \(\varrho ,\eta , {\mathbb {E}}\) obtained by the VFV method
Full size table

Example 4.3

This experiment is used to further test our theoretical analysis. It describes left-going and right-going rarefaction waves, whose initial data are given by

$$\begin{aligned} (\varrho , u , p)(0,x) \; = \; {\left\{ \begin{array}{ll} (1 ,\, -2 ,\, 0.4 ) , &{} x < 0.5, \\ (1 ,\, 2 ,\, 0.4 ) , &{} x > 0.5. \end{array}\right. } \end{aligned}$$

Figure 4a and b show the density \(\varrho \) obtained at \(T = 0.15\) by the Godunov method and the VFV method, respectively. Moreover, the corresponding \(L^2\)-error of \((\varrho , {\varvec{m}}, \eta )\) as well as the \(L^1\)-norm of \({\mathbb {E}}\) are shown in Fig.  4c and d, see also Table 4.

In our numerical results the convergence rate is approximately 1/2 (resp. 1) for \((\varrho , {\varvec{m}}, \eta )\) (resp. \({\mathbb {E}}\)), which is consistent with our theoretical analysis.

Fig. 4
figure 4

Example 4.3: density \(\varrho \) (top) and the errors of \(\varrho , {\varvec{m}}, \eta , {\mathbb {E}}\) (bottom) obtained by the Godunov method (left) and the VFV method (right). The black solid lines without any marker in the last two subfigures denote the reference slope of \(h^{1/2}\)

Full size image

Example 4.4

This experiment is devoted to the 1D Sod problem. Our aim is to test the convergence rate when the exact solution consists of rarefaction, contact and shock waves. In this example the final time is set to \(T=0.15\) and the initial data are given by

$$\begin{aligned} (\varrho , u , p)(0,x) \; = \; {\left\{ \begin{array}{ll} (1 ,\, 0 ,\, 1 ) , &{} x < 0.5, \\ (0.125 ,\, 0 ,\, 0.1 ) , &{} x > 0.5. \end{array}\right. } \end{aligned}$$

Figure 5a and b show the density \(\varrho \) obtained by the Godunov and VFV methods on different meshes. Moreover, the errors of \((\varrho , {\varvec{m}}, \eta )\) and \({\mathbb {E}}\) are shown in Fig. 5c and d, respectively, see also Table 5 for more details.

These numerical results indicate that the convergence rate of \((\varrho , {\varvec{m}}, \eta )\) (resp. \({\mathbb {E}}\)) on discontinuities is reduced; it is between 1/4 and 1/2 (resp. between 1/2 and 1).

4.2 Two Dimensional Experiments

In this section we present four two-dimensional Riemann problems. The computational domain is taken as \([0,1]^2\). Here the solution \({\widetilde{{\varvec{U}}}}\) used in the relative energy is taken as the reference solution computed by the corresponding numerical method on the uniform mesh with \(4096^2\) cells.

Example 4.5

The first 2D Riemann problem describes the interaction of four rarefaction waves. The initial data are given by

$$\begin{aligned} (\varrho , {\varvec{u}}, p)(0, x) \; = \; {\left\{ \begin{array}{ll} (1 ,\, 0 ,\, 0 ,\, 1 ) , &{} x_1> 0.5,\,x_2>0.5, \\ (0.5197 ,\, -0.7259,\, 0 , \, 0.4) , &{} x_1<0.5,\,x_2>0.5, \\ (1 ,\, -0.7259,\, -0. 7259, \, 1) , &{} x_1<0.5,\,x_2<0.5, \\ (0.5197 ,\, 0,\, -0.7259 , \, 0.4) , &{} x_1>0.5,\,x_2<0.5. \end{array}\right. } \end{aligned}$$

In this example the final time is set to \(T = 0.2\). Figure 6a and b show the density isolines obtained by the Godunov and VFV methods on a mesh with \(4096^2\) cells. Moreover, Fig. 6c and d show the \(L^2\)-errors of \((\varrho , {\varvec{m}}, \eta )\) and \(L^1\)-norm of \({\mathbb {E}}\) on successively refined meshes, see Table 6 for details.

In this case the convergence rates of \((\varrho , {\varvec{m}}, \eta )\) (resp. \({\mathbb {E}}\)) are slightly better than 1/2 (resp. 1). This may indicate that our rigorous error estimates are suboptimal in the case of finitely many rarefaction waves.

Table 4 Example 4.3: Errors and convergence rates of \(\varrho , {\varvec{m}},\eta ,{\mathbb {E}}\) obtained by the Godunov and VFV methods
Full size table
Fig. 5
figure 5

Example 4.4: density \(\varrho \) (top) and the errors of \(\varrho , {\varvec{m}}, \eta , {\mathbb {E}}\) (bottom) obtained by the Godunov method (left) and the VFV method (right). The black solid lines without any marker in the last two subfigures denote the reference slope of \(h^{1/2}\)

Full size image
Table 5 Example 4.4: Errors and convergence rates of \(\varrho , {\varvec{m}},\eta ,{\mathbb {E}}\) obtain the Godunov and VFV methods
Full size table
Fig. 6
figure 6

Example 4.5: density isolines (top) on a mesh with \(4096^2\) cells and the errors of \(\varrho , {\varvec{m}}, \eta , {\mathbb {E}}\) (bottom) obtained by the Godunov method (left) and the VFV method (right). The black solid lines without any marker in the last two subfigures denote the reference slope of \(h^{1/2}\)

Full size image
Table 6 Example 4.5: errors and convergence rates of \(\varrho , {\varvec{m}},\eta ,{\mathbb {E}}\) obtained by the Godunov and VFV methods
Full size table

Example 4.6

The initial data of the second 2D Riemann problem are given by

$$\begin{aligned} (\varrho , {\varvec{u}}, p)(0, x) = {\left\{ \begin{array}{ll} (0.5 ,\, 0.5 ,\, -0.5 ,\, 5 ) , &{} x_1> 0.5,\,x_2>0.5, \\ (1 ,\, 0.5,\, 0.5 , \, 5) , &{} x_1<0.5,\,x_2>0.5, \\ (2 ,\, -0.5,\, 0.5 , \, 5) , &{} x_1<0.5,\,x_2<0.5, \\ (1.5 ,\, -0.5,\, -0.5 , \, 5) , &{} x_1>0.5,\,x_2<0.5. \end{array}\right. } \end{aligned}$$

The exact solution consists of four interacting contact discontinuities yielding vortex sheets with negative signs. We simulate till \(T=0.2\). Figure 7a and b show the density isolines obtained by the Godunov and VFV methods on a mesh with \(4096^2\) cells. The \(L^2\)-errors of \((\varrho , {\varvec{m}}, \eta )\) as well as the \(L^1\)-norm of \({\mathbb {E}}\) are shown in Fig. 7c and d, see also Table 7.

The numerical results indicate that \((\varrho , {\varvec{m}}, \eta )\) converges with a convergence rate about 1/2 and the convergence rate for \({\mathbb {E}}\) is approximately 1.

Fig. 7
figure 7

Example 4.6: density isolines (top) on a mesh with \(4096^2\) cells and the errors of \(\varrho , {\varvec{m}}, \eta , {\mathbb {E}}\) (bottom) obtained by the Godunov method (left) and the VFV method (right). The black solid lines without any marker in the last two subfigures denote the reference slope of \(h^{1/2}\)

Full size image
Table 7 Example 4.6: Errors and convergence rates of \(\varrho , {\varvec{m}},\eta ,{\mathbb {E}}\) obtained by the Godunov and VFV methods
Full size table

Example 4.7

The initial data of the third 2D Riemann problem are given by

$$\begin{aligned} (\varrho , {\varvec{u}}, p)(0, x) \; = \; {\left\{ \begin{array}{ll} (1.5 ,\, 0 ,\, 0 ,\, 1.5 ) , &{} x_1> 0.5,\,x_2>0.5, \\ (0.5323 ,\, 1.206,\, 0 , \, 0.3) , &{} x_1<0.5,\,x_2>0.5,\\ (0.138 ,\, 1.206,\, 1.206, \, 0.029) , &{} x_1<0.5,\,x_2<0.5,\\ (0.5323 ,\, 0,\, 1.206 , \, 0.3) , &{} x_1>0.5,\,x_2<0.5, \end{array}\right. } \end{aligned}$$

which describes the interaction of four shock waves. In this example the final time is set to \(T = 0.35\). Figure 8 shows the density isolines on a mesh with \(4096^2\) cells and the errors of \((\varrho , {\varvec{m}}, \eta )\) and \({\mathbb {E}}\) obtained on successively refined meshes. Table 8 lists the errors and convergence rates.

This example indicates that \((\varrho , {\varvec{m}}, \eta )\) converges with a ratio between 1/4 and 1/2 and \({\mathbb {E}}\) converges with a ratio between 1/2 and 1.

Fig. 8
figure 8

Example 4.7: density isolines (top) on a mesh with \(4096^2\) cells and the errors of \(\varrho , {\varvec{m}}, \eta , {\mathbb {E}}\) (bottom) obtained by the Godunov method (left) and the VFV method (right). The black solid lines without any marker in the last two subfigures denote the reference slope of \(h^{1/2}\)

Full size image
Table 8 Example 4.7: Errors and convergence rates of \(\varrho , {\varvec{m}},\eta ,{\mathbb {E}}\) obtained by the Godunov and VFV methods
Full size table

Example 4.8

The initial data of the fourth 2D Riemann problem are given by

$$\begin{aligned} (\varrho , {\varvec{u}}, p)(0, x) \; = \; {\left\{ \begin{array}{ll} (0.5313 ,\, 0 ,\, 0 ,\, 0.4 ) , &{} x_1> 0.5,\,x_2>0.5, \\ (1 ,\, 0.7276,\, 0 , \, 1) , &{} x_1<0.5,\,x_2>0.5, \\ (0.8 ,\, 0,\, 0, \, 1) , &{} x_1<0.5,\,x_2<0.5, \\ (1 ,\, 0,\, 0.7276 , \, 1) , &{} x_1>0.5,\,x_2<0.5. \end{array}\right. } \end{aligned}$$

This experiment describes the interaction of four discontinuities (the left and bottom discontinuities are two contact discontinuities and the top and right are two shock waves). The final time is set to \(T=0.25\). Figure 9 shows the density isolines obtained by the Godunov and VFV methods on a mesh with \(4096^2\) cells. The \(L^2\)-errors of \(\varrho , {\varvec{m}}, \eta \), and the \(L^1\)-norm of \({\mathbb {E}}\) are presented in Fig. 9 and Table 9.

These numerical results indicate a convergence rate around 1/2 for the \(L^2\)-error of \((\varrho , {\varvec{m}}, \eta )\) and a rate around 1 for the \(L^1\)-norm of the relative energy \({\mathbb {E}}\).

Fig. 9
figure 9

Example 4.8: density isolines (top) on a mesh with \(4096^2\) cells and the errors of \(\varrho , {\varvec{m}}, \eta , {\mathbb {E}}\) (bottom) obtained by the Godunov method (left) and the VFV method (right). The black solid lines without any marker in the last two subfigures denote the reference slope of \(h^{1/2}\)

Full size image
Table 9 Example 4.8: Errors and convergence rates of \(\varrho , {\varvec{m}},\eta ,{\mathbb {E}}\) obtained by the Godunov and VFV methods
Full size table

5 Conclusion

In this paper we have analyzed a priori error estimates between the numerical solution obtained by the Godunov method and the strong exact solution of the multidimensional Euler system via the relative energy. Assuming that there exist a uniform positive lower bound on the density and a positive upper bound on the energy, we showed that the \(L^1\)-norm of the relative energy is equivalent to the square of the \(L^2\)-norm of the error of the numerical solution, see (2.17). Recalling the consistency formulation proved in [22] and applying Gronwall’s lemma, we have derived the estimates for the relative energy in Theorem 3.1. Specifically, the relative energy converges at least at the rate of 1/2 in the \(L^1\)-norm. At the same time, the density, momentum and entropy converge at least at the rate of 1/4 in the \(L^2\)-norm. Being inspired by the fact that the Godunov method for scalar conservation laws has bounded total variations we have formulated an additional hypothesis (3.9). If we assume that (3.9) holds, the convergence rate of density, momentum and entropy (resp. relative energy) can be improved to at least 1/2 (resp. 1), see Theorem 3.3. Finally, we pointed out that our theoretical analysis rigorously holds only for strong solutions, e.g. for a solution that contains only finitely many rarefaction waves.

We have experimentally computed convergence rates for several one- and two-dimensional Riemann problems. From Examples 4.1 and 4.3 containing only rarefaction waves, we observed that the convergence rates of density, momentum and entropy (resp. relative energy) are slightly better than 1/2 (resp. 1), which is consistent with the theoretical results presented in Theorem 3.3. In our experiments for the Riemann problem the Godunov and VFV methods have a convergence rate about 1/4 for the contact wave and about 1/2 for the shock wave, respectively. In future it will be interesting to analyze theoretically the convergence rate towards a unique weak entropy solution containing shock and contact waves.