The gradient discretisation method for linear advection problems

TL;DR

This paper extends the Gradient Discretisation Method (GDM) to linear hyperbolic equations, unifying and analyzing various numerical schemes like finite elements and finite volumes for scalar advection problems.

Contribution

It adapts GDM to hyperbolic equations, enabling unified design and convergence analysis of multiple numerical schemes for scalar advection.

Findings

Convergence of the adapted GDM scheme is established.

Numerical tests confirm the effectiveness of the method.

The scheme accommodates various discretisation techniques.

Abstract

We adapt the Gradient Discretisation Method (GDM), originally designed for elliptic and parabolic partial differential equations, to the case of a linear scalar hyperbolic equations. This enables the simultaneous design and convergence analysis of various numerical schemes, corresponding to the methods known to be GDMs, such as finite elements (conforming or non-conforming, standard or mass-lumped), finite volumes on rectangular or simplicial grids, and other recent methods developed for general polytopal meshes. The scheme is of centred type, with added linear or non-linear numerical diffusion. We complement the convergence analysis with numerical tests based on the mass-lumped P1 conforming and non conforming finite element and on the hybrid finite volume method.

Tables

Table 1. Table 1: Case 1, results with the centred scheme, using the CVFE method, δ ​ t = 0.4 ​ h 𝛿 𝑡 0.4 ℎ {\delta\!t}=0.4h

$h$	errl2	rate	errl1	rate	umin	umax
0.250	2.95E-01	-	1.95E-01	-	0.108	0.137
0.125	2.55E-01	0.212	1.37E-01	0.504	0.014	0.174
0.062	2.32E-01	0.136	1.23E-01	0.158	0.000	0.344
0.031	1.77E-01	0.394	8.55E-02	0.525	-0.001	0.734
0.016	1.23E-01	0.524	4.73E-02	0.853	-0.013	1.003

Table 2. Table 2: Case 1, results with the centred scheme, using the MLNC– ℙ 1 subscript ℙ 1 \mathbb{P}_{1} method, δ ​ t = 0.4 ​ h 𝛿 𝑡 0.4 ℎ {\delta\!t}=0.4h

$h$	errl2	rate	errl1	rate	umin	umax
0.250	2.52E-1	-	1.10E-1	-	0.043	0.054
0.125	2.65E-1	-0.076	1.51E-1	-0.457	0.016	0.194
0.062	2.37E-1	0.165	1.31E-1	0.208	0.000	0.361
0.031	1.82E-1	0.381	8.64E-2	0.597	0.000	0.687
0.016	1.33E-1	0.456	5.34E-2	0.694	0.000	0.960

Table 3. Table 3: Case 1, results with the centred scheme, using the HFV method, δ ​ t = 0.4 ​ h 𝛿 𝑡 0.4 ℎ {\delta\!t}=0.4h

h	errl2	rate	errl1	rate	umin	umax
0.35	2.80E-1	-	2.08E-1	-	0.152	0.155
0.18	2.79E-1	0.001	1.54E-1	0.436	0.044	0.124
0.09	2.59E-1	0.111	1.30E-1	0.236	0.001	0.220
0.04	2.10E-1	0.300	1.08E-1	0.276	0.000	0.499
0.02	1.47E-1	0.520	6.57E-2	0.713	0.000	0.906

Table 4. Table 4: Case 1, results with the upstream ℙ 1 subscript ℙ 1 \mathbb{P}_{1} scheme, δ ​ t = 0.4 ​ h 𝛿 𝑡 0.4 ℎ {\delta\!t}=0.4h

$h$	errl2	rate	errl1	rate	umin	umax
0.250	2.59E-01	-	1.65E-01	-	0.005	0.313
0.125	2.32E-01	0.159	1.19E-01	0.462	0.000	0.286
0.062	2.13E-01	0.122	1.10E-01	0.122	0.000	0.454
0.031	1.85E-01	0.205	9.13E-02	0.266	0.000	0.672
0.016	1.53E-01	0.270	6.93E-02	0.398	0.000	0.868

Table 5. Table 5: Case 2, results with Scheme ( 21 ), using the CVFE method, δ ​ t = 0.4 ​ h 𝛿 𝑡 0.4 ℎ {\delta\!t}=0.4h

$h$	errl2	rate	errl1	rate	errl$\infty$	rate
0.250	4.96E-02	-	4.34E-02	-	0.138	-
0.125	1.82E-02	1.44	1.42E-02	1.61	7.17E-02	0.94
0.062	5.89E-03	1.62	4.26E-03	1.73	3.59E-02	0.99
0.031	1.81E-03	1.70	1.16E-03	1.87	1.78E-02	1.01
0.016	5.51E-04	1.71	3.06E-04	1.92	8.86E-03	1.01

Table 6. Table 6: Case 2, results with the upstream scheme, using the CVFE method, δ ​ t = 0.4 ​ h 𝛿 𝑡 0.4 ℎ {\delta\!t}=0.4h

$h$	errl2	rate	errl1	rate	errl$\infty$	rate
0.250	5.37E-02	-	5.01E-02	-	9.25E-02	-
0.125	2.88E-02	0.90	2.59E-02	0.95	5.85E-02	0.66
0.062	1.57E-02	0.87	1.35E-02	0.94	3.36E-02	0.80
0.031	8.55E-03	0.88	6.94E-03	0.96	2.21E-02	0.60
0.016	4.57E-03	0.90	3.53E-03	0.98	1.50E-02	0.55

Equations

$${\partial }_t\bar{u}+\mathrm{div}(\bar{u}\vec{\mathbf{v}})+\bar{u}{q}^P=f{q}^I,in\; \Omega \times (0,T),$$

$$\bar{u}(\mathbf{x},0)=u_{\mathrm{ini}}(\mathbf{x}),for\; a.e.x\in \Omega ,$$

$$\Omega \text{mbox}isanopenboundedconnectedpolyhedralsubsetof{\mathbb{R}}^d,d\in {\mathbb{N}}^{\star }\text{mbox}andT>0,$$

$$u_{\mathrm{ini}}\in L^2(\Omega )\text{mbox}andf\in L^2(\Omega \times (0,T)),$$

$$q^I,q^P\in L^{\infty }(\Omega \times (0,T))\text{mbox}with{q}^I\ge 0\text{mbox}and{q}^P\ge 0\text{mbox}a.e.in\Omega \times (0,T),$$

$$\vec{\mathbf{v}}\in W^{1,\infty }(\Omega \times (0,T){)}^d\text{mbox}satisfies\vec{\mathbf{v}}\cdot \mathbf{n}=0\; on\partial \Omega \times (0,T)\text{mbox}and\mathrm{div}\vec{\mathbf{v}}=q^I-q^P\text{mbox}a.e.in\Omega \times (0,T),$$

$$\displaystyle\bar{u}\in L^{2}(\Omega\times(0,T))\mbox{ and, for all $\varphi\in C^{\infty}_{c}(\mathbb{R}^{d}\times[0,T))$,}$$

$$-{\int }_0^T{\int }_{\Omega }\bar{u}{\partial }_t\phi \mathrm{d}\mathbf{x}\mathrm{d}t-{\int }_{\Omega }{u}_{\mathrm{ini}}(\mathbf{x})\phi (\mathbf{x},0)\mathrm{d}\mathbf{x}+{\int }_0^T{\int }_{\Omega }(-\bar{u}\vec{\mathbf{v}}\cdot \nabla \phi +\bar{u}{q}^P\phi )\mathrm{d}\mathbf{x}\mathrm{d}t={\int }_0^T{\int }_{\Omega }f{q}^I\phi \mathrm{d}\mathbf{x}\mathrm{d}t,$$

$$\displaystyle u^{(n+1)}\in V,~{}u^{(n+\theta)}=\theta u^{(n+1)}+(1-\theta)u^{(n)}\mbox{ and, for all $v\in V$,}$$

$${\int }_{\Omega }\frac{u^{(n+1)}-u^{(n)}}{{\delta t}^{(n+\frac{1}{2})}}v\mathrm{d}\mathbf{x}+{\int }_{\Omega }(-u^{(n+\theta )}{\vec{\mathbf{v}}}^{(n+\frac{1}{2})}\cdot \nabla v+u^{(n+\theta )}(q^P{)}^{(n+\frac{1}{2})}v)\mathrm{d}\mathbf{x}={\int }_{\Omega }{f}^{(n+\frac{1}{2})}(q^I{)}^{(n+\frac{1}{2})}v\mathrm{d}\mathbf{x},$$

$${\int }_{\Omega }{u}^{(n+\theta )}\vec{\mathbf{v}}\cdot \nabla u^{(n+\theta )}\mathrm{d}\mathbf{x}={\int }_{\Omega }\vec{\mathbf{v}}\cdot \nabla \frac{(u^{(n+\theta )}{)}^2}{2}\mathrm{d}\mathbf{x}=-{\int }_{\Omega }\frac{(u^{(n+\theta )}{)}^2}{2}\mathrm{div}\vec{\mathbf{v}}\mathrm{d}\mathbf{x}.$$

$$\displaystyle u^{(n+1)}\in X_{{\cal D}},~{}u^{(n+\theta)}=\theta u^{(n+1)}+(1-\theta)u^{(n)}\mbox{ and, for all $v\in X_{{\cal D}}$,}$$

$${\int }_{\Omega }{\Pi }_{\mathcal{D}}\frac{u^{(n+1)}-u^{(n)}}{{\delta t}^{(n+\frac{1}{2})}}{\Pi }_{\mathcal{D}}v\mathrm{d}\mathbf{x}$$

$$+{\int }_{\Omega }(-{\Pi }_{\mathcal{D}}{u}^{(n+\theta )}{\vec{\mathbf{v}}}^{(n+\frac{1}{2})}\cdot {\nabla }_{\mathcal{D}}v+{\Pi }_{\mathcal{D}}{u}^{(n+\theta )}(q^P{)}^{(n+\frac{1}{2})}{\Pi }_{\mathcal{D}}v)\mathrm{d}\mathbf{x}={\int }_{\Omega }{f}^{(n+\frac{1}{2})}(q^I{)}^{(n+\frac{1}{2})}{\Pi }_{\mathcal{D}}v\mathrm{d}\mathbf{x}.$$

$$\frac{1}{2}\mathrm{div}\vec{\mathbf{v}}+q^P=\frac{1}{2}\left(q^I-q^P\right)+q^P=\frac{1}{2}(q^I+q^P),$$

$$\forall v\in C_c^{\infty }({\mathbb{R}}^d\times [0,T)),$$

$$-{\int }_0^T{\int }_{\Omega }\bar{u}{\partial }_{t}v\mathrm{d}\mathbf{x}\mathrm{d}t-{\int }_{\Omega }{u}_{\mathrm{ini}}(\mathbf{x})v(\mathbf{x},0)\mathrm{d}\mathbf{x}$$

$$+{\int }_0^T{\int }_{\Omega }\left(\frac{1}{2}\nabla \bar{u}\cdot \vec{\mathbf{v}}v-\frac{1}{2}\bar{u}\vec{\mathbf{v}}\cdot \nabla v+\frac{1}{2}\bar{u}(q^I+q^P)v\right)\mathrm{d}\mathbf{x}\mathrm{d}t={\int }_0^T{\int }_{\Omega }f{q}^{I}v\mathrm{d}\mathbf{x}\mathrm{d}t.$$

$$\displaystyle u^{(n+1)}\in X_{{\cal D}},~{}u^{(n+\theta)}=\theta u^{(n+1)}+(1-\theta)u^{(n)}\mbox{ and, for all $v\in X_{{\cal D}}$,}$$

$$\displaystyle\int_{\Omega}\Pi_{\cal D}\dfrac{u^{(n+1)}-u^{(n)}}{{\delta\!t}^{(n+{\frac{1}{2}})}}\Pi_{\cal D}v\ {\rm d}{\bm{x}}+\int_{\Omega}\Big{(}{\frac{1}{2}}\nabla_{\cal D}u^{(n+\theta)}\cdot\vec{{\bm{v}}}^{(n+{\frac{1}{2}})}\Pi_{\cal D}v-{\frac{1}{2}}\Pi_{\cal D}u^{(n+\theta)}\vec{{\bm{v}}}^{(n+{\frac{1}{2}})}\cdot\nabla_{\cal D}v$$

$$\displaystyle\qquad\qquad\qquad+{\frac{1}{2}}\Pi_{\cal D}u^{(n+\theta)}\left[(q^{I})^{(n+{\frac{1}{2}})}+(q^{P})^{(n+{\frac{1}{2}})}\right]\Pi_{\cal D}v\Big{)}{\rm d}{\bm{x}}=\int_{\Omega}f^{(n+{\frac{1}{2}})}\ (q^{I})^{(n+{\frac{1}{2}})}\Pi_{\cal D}v\ {\rm d}{\bm{x}}.$$

$$\mathrm{div}\left(\left[\bar{u}-h\frac{\vec{\mathbf{v}}}{|\vec{\mathbf{v}}|}\cdot \nabla \bar{u}\right]\vec{\mathbf{v}}\right)=\mathrm{div}(\bar{u}\vec{\mathbf{v}}-h\Lambda \nabla \bar{u}),with\; \Lambda (\mathbf{x},t)=\frac{\vec{\mathbf{v}}(\mathbf{x},t)}{|\vec{\mathbf{v}}(\mathbf{x},t)|}\otimes \vec{\mathbf{v}}(\mathbf{x},t).$$

$$\Lambda (\mathbf{x},t)=\frac{\vec{\mathbf{v}}(\mathbf{x},t)}{|\vec{\mathbf{v}}(\mathbf{x},t)|}\otimes \vec{\mathbf{v}}(\mathbf{x},t)+\mu \mathrm{Id},$$

$$-h^{\alpha }\mathrm{div}(\Lambda |\nabla \bar{u}{|}_{\Lambda }^{p-2}\nabla \bar{u})\text{mbox}where|\nabla \bar{u}{|}_{\Lambda }=\sqrt{\Lambda \nabla \bar{u}\cdot \nabla \bar{u}},$$

$$\displaystyle\bar{u}\in L^{2}(\mathbb{R}^{d}\times(0,T))\mbox{ and, for all $\varphi\in C^{\infty}_{c}(\mathbb{R}^{d}\times[0,T))$,}$$

$$-{\int }_0^T{\int }_{{\mathbb{R}}^d}\bar{u}{\partial }_t\phi \mathrm{d}\mathbf{x}\mathrm{d}t-{\int }_{{\mathbb{R}}^d}{u}_{\mathrm{ini}}(\mathbf{x})\phi (\mathbf{x},0)\mathrm{d}\mathbf{x}$$

$$+{\int }_0^T{\int }_{{\mathbb{R}}^d}(-\bar{u}\vec{\mathbf{v}}\cdot \nabla \phi +\bar{u}{q}^P\phi )\mathrm{d}\mathbf{x}\mathrm{d}t={\int }_0^T{\int }_{{\mathbb{R}}^d}f{q}^I\phi \mathrm{d}\mathbf{x}\mathrm{d}t.$$

$$-{\int }_0^T{\Theta }^{\prime }(t){\bar{U}}_{\psi }(t)\mathrm{d}t-\Theta (0){\int }_{{\mathbb{R}}^d}{u}_{\mathrm{ini}}(\mathbf{x})\psi (\mathbf{x})\mathrm{d}\mathbf{x}$$

$$+{\int }_0^T\Theta (t){\int }_{{\mathbb{R}}^d}(-\bar{u}(\mathbf{x},t)\vec{\mathbf{v}}(\mathbf{x},t)\cdot \nabla \psi (\mathbf{x})+\bar{u}(\mathbf{x},t)q^P(\mathbf{x},t)\psi (\mathbf{x}))\mathrm{d}\mathbf{x}\mathrm{d}t$$

$$={\int }_0^T\Theta (t){\int }_{{\mathbb{R}}^d}f(\mathbf{x},t)q^I(\mathbf{x},t)\psi (\mathbf{x})\mathrm{d}\mathbf{x}\mathrm{d}t.$$

$$\bar{U}_{\psi}^{\prime}(t)=\int_{\mathbb{R}^{d}}\Big{(}(f({\bm{x}},t)\ q^{I}({\bm{x}},t)-\bar{u}({\bm{x}},t)\ q^{P}({\bm{x}},t))\psi({\bm{x}})\ +\bar{u}({\bm{x}},t)\ \vec{{\bm{v}}}({\bm{x}},t)\cdot\nabla\psi({\bm{x}})\Big{)}{\rm d}{\bm{x}}.$$

$$\frac{1}{2}{\int }_0^T{\int }_{\Omega }(\bar{u}(\mathbf{x},t){)}^2\mathrm{d}\mathbf{x}\mathrm{d}t+\frac{1}{2}{\int }_0^T(T-t){\int }_{\Omega }(q^I(\mathbf{x},t)+q^P(\mathbf{x},t))\bar{u}(\mathbf{x},t{)}^2\mathrm{d}\mathbf{x}\mathrm{d}t$$

$$=\frac{T}{2}{\int }_{\Omega }{u}_{\mathrm{ini}}(\mathbf{x}{)}^2\mathrm{d}\mathbf{x}+{\int }_0^T(T-t){\int }_{\Omega }f(\mathbf{x},t)q^I(\mathbf{x},t)\bar{u}(\mathbf{x},t)\mathrm{d}\mathbf{x}\mathrm{d}t.$$

$$\forall \mathbf{x},t\in {\mathbb{R}}^d\times [0,T],\phi (\mathbf{x},t)=(T-t){\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{R}}^d}\bar{u}(\mathbf{z},t){\rho }_n(\mathbf{y}-\mathbf{z}){\rho }_n(\mathbf{y}-\mathbf{x})\mathrm{d}\mathbf{z}\mathrm{d}\mathbf{y},$$

$$\forall \mathbf{x},t\in {\mathbb{R}}^d\times [0,T],\nabla \phi (\mathbf{x},t)=(T-t){\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{R}}^d}\bar{u}(\mathbf{z},t)\nabla {\rho }_n(\mathbf{y}-\mathbf{z}){\rho }_n(\mathbf{y}-\mathbf{x})\mathrm{d}\mathbf{z}\mathrm{d}\mathbf{y},$$

$$\text{Cr}{t1}^{(n)}={\int }_0^T{\int }_{{\mathbb{R}}^d}\bar{u}(\mathbf{x},t){\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{R}}^d}\bar{u}(\mathbf{z},t){\rho }_n(\mathbf{y}-\mathbf{z}){\rho }_n(\mathbf{y}-\mathbf{x})\mathrm{d}\mathbf{z}\mathrm{d}\mathbf{y}\mathrm{d}\mathbf{x}\mathrm{d}t,$$

$$\displaystyle\Cr{t2}^{(n)}=-\int_{0}^{T}(T-t)\int_{\mathbb{R}^{d}}\bar{u}({\bm{x}},t)\int_{\mathbb{R}^{d}}\partial_{t}\Big{(}\int_{\mathbb{R}^{d}}\bar{u}({\bm{z}},t)\rho_{n}({\bm{y}}-{\bm{z}})\rho_{n}({\bm{y}}-{\bm{x}}){\rm d}{\bm{z}}\Big{)}{\rm d}{\bm{y}}{\rm d}{\bm{x}}{\rm d}t,$$

$$\text{Cr}{t3}^{(n)}=-T{\int }_{{\mathbb{R}}^d}{u}_{\mathrm{ini}}(\mathbf{x}){\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{R}}^d}{u}_{\mathrm{ini}}(\mathbf{z}){\rho }_n(\mathbf{y}-\mathbf{z}){\rho }_n(\mathbf{y}-\mathbf{x})\mathrm{d}\mathbf{z}\mathrm{d}\mathbf{y}\mathrm{d}\mathbf{x},$$

$$\text{Cr}{t4}^{(n)}=-{\int }_0^T(T-t){\int }_{{\mathbb{R}}^d}\bar{u}(\mathbf{x},t)\vec{\mathbf{v}}(\mathbf{y},t)\cdot {\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{R}}^d}\bar{u}(\mathbf{z},t)\nabla {\rho }_n(\mathbf{y}-\mathbf{z}){\rho }_n(\mathbf{y}-\mathbf{x})\mathrm{d}\mathbf{z}\mathrm{d}\mathbf{y}\mathrm{d}\mathbf{x}\mathrm{d}t,$$

$$\text{Cr}{t4a}^{(n)}={\int }_0^T(T-t){\int }_{{\mathbb{R}}^d}(-\bar{u}(\mathbf{x},t)(\vec{\mathbf{v}}(\mathbf{x},t)-\vec{\mathbf{v}}(\mathbf{y},t))\cdot {\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{R}}^d}\bar{u}(\mathbf{z},t)\nabla {\rho }_n(\mathbf{y}-\mathbf{z}){\rho }_n(\mathbf{y}-\mathbf{x})\mathrm{d}\mathbf{z}\mathrm{d}\mathbf{y}\mathrm{d}\mathbf{x}\mathrm{d}t,$$

$$\text{Cr}{t5}^{(n)}={\int }_0^T(T-t){\int }_{{\mathbb{R}}^d}\bar{u}(\mathbf{x},t)q^P(\mathbf{x},t){\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{R}}^d}\bar{u}(\mathbf{z},t){\rho }_n(\mathbf{y}-\mathbf{z}){\rho }_n(\mathbf{y}-\mathbf{x})\mathrm{d}\mathbf{z}\mathrm{d}\mathbf{y}\mathrm{d}\mathbf{x}\mathrm{d}t,$$

Full text

Keywords: linear scalar hyperbolic equation, Gradient Discretisation Method, convergence analysis, numerical tests.

AMS subject classification: 65N12, 65N30

1 Introduction

We are interested here in designing and analysing an approximation of $\bar{u}$, solution to the linear advection problem stated in its strong form as

[TABLE]

with the following assumptions on the data:

[TABLE]

where $\bm{n}$ is the outer normal to $\partial\Omega$. Since the normal boundary value of $\vec{{\bm{v}}}$ vanishes, there is no need for a boundary condition on (1a).

The model (1) typically arises in oil recovery from underground reservoirs [1, 15] or in underground water resources management [24], in which case $q^{I}$ and $q^{P}$ may represent the injection and production wells and $\bar{u}$ is the concentration of injected solvent or pollutant. The problem (1) is often discretised by the upstream weighting finite volume scheme (see, for example, [16, Chapters 5 and 6] and references therein), which is easy to implement even on unstructured meshes since the problem is first order. There are also numerous papers studying Galerkin methods for this type of problems, which are based on the following weak formulation: a function $\bar{u}$ is said to be a weak solution of Problem (1) if:

[TABLE]

where $C^{\infty}_{c}(\mathbb{R}^{d}\times[0,T))$ is the set of the restrictions of functions of $C^{\infty}_{c}(\mathbb{R}^{d}\times(-\infty,T))$ to $\mathbb{R}^{d}\times[0,T)$.

Let $0=t^{(0)}<t^{(1)}<\cdots<t^{(N)}=T$ be a discretisation of the time interval, and let ${\delta\!t}^{(n+{\frac{1}{2}})}=t^{(n+1)}-t^{(n)}$. We recall that, for $V\subset H^{1}(\Omega)$ a finite dimensional space and $\theta\in[0,1]$, the $\theta$-scheme takes the following form: $u^{(0)}\in V$ being a chosen interpolate of $u_{\rm ini}$, the scheme consists in finding, for all $n=0,\ldots,N-1$,

[TABLE]

with suitable time approximations of the data indexed by $(n+{\frac{1}{2}})$. This scheme is $L^{2}$ stable provided that $\theta\geq{\frac{1}{2}}$, which is proved letting $v=u^{(n+\theta)}$ and following the calculus formula

[TABLE]

Weak convergence properties are then obtained for the approximate solution, which generally displays oscillations. See [14] for a complete study of the particular case of Finite element methods, and [7] for a comparison of different Galerkin schemes. A convergence result is proved in [13] under strong regularity hypotheses on the solution and with a constant velocity field.

This paper is focused on the case where the approximation of $u$ is no longer done in a subspace of $H^{1}(\Omega)$. In a number of situations, coupled problems including terms of different nature (e.g. diffusive, advective…) must be solved in an industrial context where the discretisation method, imposed by the use of an existing code, is based on non conforming finite element, discontinuous Galerkin or hybrid methods (with face and cell unknowns), for example.

In order to handle such a situation, we use the Gradient Discretisation Method (GDM) framework, which gives a unified formulation of a large class of conforming and nonconforming methods; we refer the reader to the monograph [12] for details. The idea of the GDM is to replace, in a weak formulation of the continuous problem, the continuous space by the vector space of the degrees of freedom of the method $X_{\cal D}$, the functions $u$ and $v$ by their reconstruction $\Pi_{\cal D}u$ and $\Pi_{\cal D}v$, and the gradient $\nabla v$ by the reconstruction of a discrete gradient $\nabla_{\cal D}v$. For conforming methods, $\Pi_{\cal D}(X_{\cal D})$ is a subspace of $H^{1}(\Omega)$ and, for $v\in X_{\cal D}$, $\nabla_{\cal D}v=\nabla(\Pi_{\cal D}v)$; for non-conforming finite element methods, $\Pi_{\cal D}(X_{\cal D})$ is a space of piecewise polynomial functions and, for all $v\in X_{\cal D}$, $\nabla_{\cal D}v$ is the broken gradient of $\Pi_{\cal D}v$. Discontinuous Galerkin methods, which are popular in the framework of hyperbolic problems, can also be embedded in the GDM; for these methods, $\Pi_{\cal D}(X_{\cal D})$ is again a space of piecewise polynomial functions, the expression of $\nabla_{\cal D}v$ takes into account both the broken gradient of $\Pi_{\cal D}v$ and the jump terms, and no additional stabilisation term has to be introduced in the formulation of the scheme (see [12, Chapter 11]). Note that for fully discrete methods or mass-lumped versions of the previous schemes, $\Pi_{\cal D}$ is a genuine function reconstruction (see the schemes used in Section 5).

A natural scheme would then be: given an interpolate $u^{(0)}\in X_{\cal D}$ of $u_{\rm ini}$, solve for $n=0,\ldots,N-1$,

[TABLE]

Unfortunately, it does not seem possible to establish the stability (and thus the convergence) of (6) due to the absence of the equivalent of the calculus chain (5) in this fully discrete setting involving function and gradient reconstructions $\Pi_{\cal D}$ and $\nabla_{\cal D}$ instead of the classical differential operators. To obtain a scheme amenable to a convergence analysis, we thus consider an alternative formulation, using a skew-symmetric reformulation of the advective term.

If $\nabla\bar{u}\in L^{2}(\Omega\times(0,T))$, owing to the relation

[TABLE]

a function $\bar{u}\in L^{2}(\Omega\times(0,T))$ is a solution to (3) if and only if it satisfies

[TABLE]

The idea to discretise (1a) is then to mimick the formulation (7) instead of (3) in the discrete setting (this idea is in the same line as the weak formulation chosen in [4, Hypothesis (A1)]). Indeed, similarly to the standard skew-symmetric formulation of the convective term in the Navier-Stokes equations, the advection component in (7) vanishes when the solution is taken as a test function. The GDM scheme based on (7) is thus: take $u^{(0)}\in X_{\cal D}$ and interpolant of $u_{\rm ini}$ and, for all $n=0,\ldots,N-1$,

[TABLE]

Letting $v=u^{(n+\theta)}$ in (8) leads to an estimate on $\Pi_{\cal D}u^{(n+\theta)}$, which entails a weak convergence property for the reconstruction of the function. However a new difficulty arises: the scheme (8) does not yield any estimate on $\nabla_{\cal D}u^{(n+\theta)}$; this prevents us from obtaining any limit (even weak) for this term, and thus from passing to the limit to recover the continuous problem.

This issue is solved by introducing a stabilisation term that yields a weak bound on $\nabla_{\cal D}u^{(n+\theta)}$. Several versions of such a stabilisation term can be found [22, 19], such as the symmetric linear stabilisation of [4], or the Streamline-Upwind/Petrov-Galerkin (SUPG) stabilisation [3, 23, 21, 10]. The latter is equivalent to replacing, in the term $\bar{u}\vec{{\bm{v}}}$ of (1a), $\bar{u}$ by $\bar{u}-h\frac{\vec{{\bm{v}}}}{|\vec{{\bm{v}}}|}\cdot\nabla\bar{u}$ (this is a kind of continuous upstream weighting for a mesh with size $h$). This leads to the term

[TABLE]

It is then numerically more stable to complete the SUPG scheme by modifying $\Lambda$ into

[TABLE]

for a small value $\mu>0$. This choice of stabilisation term ${\rm div}(-h\Lambda\nabla\bar{u})$ can be generalised into

[TABLE]

for some $p\in(1,+\infty)$ and $\alpha>0$, and $\Lambda({\bm{x}},t)$ symmetric positive definite with uniformly bounded eigenvalues. An obvious and easy choice is $p=2$ and $\Lambda=\rm{Id}$, which leads to the classical Laplace operator. However, using $p\neq 2$ may lead to a smaller numerical diffusion, see Section 5; let us note that in this case, the linear model (1) is approximated by a non-linear problem, which is not in general much of a problem, since the complete coupled physical model usually involves other non-linear terms. In this paper, we stabilise the scheme (8) by introducing the discrete version of the stabilisation term (9), which leads to Scheme (21). Since the GDM method also includes meshless schemes, the stabilisation term depends on a parameter $h_{\cal D}$ which is an adaptation to the hyperbolic setting of the space size of gradient discretisation for elliptic problems, see Definition 3.4 below.

In addition to providing a generic formulation that applies to a large variety of schemes, this paper presents the following original features:

The analysis applies to mesh-based as well as meshless schemes, owing to Definition 3.4 of the size of gradient discretisation which gives us a way to introduce an intrinsic vanishing viscosity without referring to any mesh size. 2. 2.

We study and compare, for different values of $p$, the effect of the stabilisation of a hyperbolic scheme by $p$-Laplace vanishing diffusion. Numerical examples show that in some cases, values of $p$ different from 2 lead to more accurate solutions. 3. 3.

The strong convergence of the stabilised scheme is obtained through an energy estimate, proved in the $L^{2}$ framework by regularisation as in [9]; this energy estimate is also used for the proof of uniqueness of the solution. 4. 4.

Convergence is established without assuming additional regularity on the solution or the velocity field, and a uniform-in-time weak convergence is proved.

This paper is organised as follows. The continuous problem and the energy estimate are studied in Section 2. We then apply in Section 3 the gradient discretisation tools to Problem (3), and derive some estimates which are used in Section 4 to establish the convergence of the scheme; as a by-product of this convergence, we also obtain an existence result for the solution to (3). In Section 5 some numerical results are provided, using three different schemes that fit into the GDM framework.

2 The continuous problem

Since the flux is null on the boundary $\partial\Omega$, the problem (3) may be reformulated on the whole space $\mathbb{R}^{d}$ by extending $\vec{{\bm{v}}}$, $q^{I}$ and $q^{P}$ to $\mathbb{R}^{d}\times\mathbb{R}$: we first choose an extension $\vec{{\bm{v}}}\in W^{1,\infty}(\mathbb{R}^{d}\times\mathbb{R})^{d}$, and then set $q^{I}=\max({\rm div}\vec{{\bm{v}}},0)$ and $q^{P}=\max(-{\rm div}\vec{{\bm{v}}},0)$ outside $\Omega\times(0,T)$. We also extend $\bar{u}$, $f$ and $u_{\rm ini}$ by the value [math] outside $\Omega\times(0,T)$ and $\Omega$ respectively. With these extensions and assuming (2), the problem (3) is equivalent to the following problem, posed on the whole space:

[TABLE]

Lemma 2.1** (Weak continuity with respect to time).**

Assuming (2), let $\bar{u}$ be a solution of (3), or to (10) after extending $\bar{u}$ by 0 outside of $\Omega$. Let $\psi\in C^{\infty}_{c}(\mathbb{R}^{d})$. Then the function $\bar{U}_{\psi}~{}:~{}t\mapsto\int_{\mathbb{R}^{d}}\bar{u}({\bm{x}},t)\psi({\bm{x}})\ {\rm d}{\bm{x}}$ satisfies $\bar{U}_{\psi}\in H^{1}(0,T)\subset C^{0}([0,T])$ and $\bar{U}_{\psi}(0)=\int_{\mathbb{R}^{d}}u_{\rm ini}({\bm{x}})\psi({\bm{x}})\ {\rm d}{\bm{x}}$. Hence, $\bar{u}\in C_{\rm w}([0,T],L^{2}(\Omega))$, where $C_{\rm w}([a,b],L^{2}(\Omega))$ stands for the space of functions $[a,b]\to L^{2}(\Omega)$ that are continuous weakly in $L^{2}(\Omega)$.

Proof.

Let $\Theta\in C^{\infty}_{c}([0,T))$. Taking $\varphi({\bm{x}},t)=\Theta(t)\psi({\bm{x}})$ in (10) yields

[TABLE]

Restricting to $\Theta\in C^{\infty}_{c}(0,T)$ this shows that, in the weak derivative sense,

[TABLE]

Since the right hand side of the above equation belongs to $L^{2}(0,T)$, this concludes the proof that $\bar{U}_{\psi}\in H^{1}(0,T)\subset C^{0}([0,T])$. The relation $\bar{U}_{\psi}(0)=\int_{\mathbb{R}^{d}}u_{\rm ini}({\bm{x}})\psi({\bm{x}})\ {\rm d}{\bm{x}}$ is proved taking $\Theta\in C^{\infty}_{c}([0,T))$ such that $\Theta(0)=1$ in (11), integrating-by-parts in time and using (12). ∎

Lemma 2.2** (Energy estimate).**

Assuming (2), any solution $\bar{u}$ of (3) satisfies:

[TABLE]

Proof.

By density of $C^{\infty}_{c}(\mathbb{R}^{d}\times[0,T))$ in $H^{1}(\mathbb{R}^{d}\times(0,T))$, we can consider functions $\varphi\in C^{0}_{c}(\mathbb{R}^{d}\times[0,T))\cap H^{1}(\mathbb{R}^{d}\times(0,T))$ in (10). Letting $\rho$ be a mollifier on $\mathbb{R}^{d}$, and $\rho_{n}({\bm{x}})=n^{d}\rho(n{\bm{x}})$ for all ${\bm{x}}\in\mathbb{R}^{d}$ and $n\in\mathbb{N}^{\star}$, we choose the function $\varphi$ defined by

[TABLE]

which satisfies $\varphi\in C^{0}_{c}(\mathbb{R}^{d}\times[0,T))\cap H^{1}(\mathbb{R}^{d}\times[0,T))$ owing to Lemma 2.1. Using an integration by parts with respect to ${\bm{y}}$, we notice that

[TABLE]

With this choice of $\varphi$ in (10) leads to $\Cl[ctrterm]{t1}^{(n)}+\Cl[ctrterm]{t2}^{(n)}+\Cl[ctrterm]{t3}^{(n)}+\Cl[ctrterm]{t4}^{(n)}+\Cl[ctrterm]{t4a}^{(n)}+\Cl[ctrterm]{t5}^{(n)}=\Cl[ctrterm]{t6}^{(n)}$, with

[TABLE]

Introducing the function $\bar{u}_{n}({\bm{y}},t)=\int_{\mathbb{R}^{d}}\bar{u}({\bm{z}},t)\rho_{n}({\bm{y}}-{\bm{z}}){\rm d}{\bm{z}}$, which converges to $\bar{u}$ in $L^{2}(\mathbb{R}^{d}\times(0,T))$ as $n\to\infty$ and satisfies $\bar{u}_{n}\in H^{1}(\mathbb{R}\times(0,T))$ and $\bar{u}_{n}({\bm{y}},0)=\int_{\mathbb{R}^{d}}u_{\rm ini}({\bm{z}})\rho_{n}({\bm{y}}-{\bm{z}}){\rm d}{\bm{z}}$ (see Lemma 2.1), we have

[TABLE]

and, using an integration-by-parts,

[TABLE]

Gathering these results leads to

[TABLE]

and therefore

[TABLE]

Turning to $\Cr{t4}^{(n)}$ we write, using the divergence formula and ${\rm div}\vec{{\bm{v}}}=q^{I}-q^{P}$,

[TABLE]

Hence,

[TABLE]

We then easily see that, as $n\to\infty$,

[TABLE]

and

[TABLE]

The proof is completed by gathering all the above convergence results and by proving that

[TABLE]

In order to do so, we follow the technique of [9, Lemma II.1] and [17, Lemma B.4]. An integration-by-parts gives

[TABLE]

with

[TABLE]

Since the function $(T-t)\bar{u}_{n}$ converges to $(T-t)\bar{u}$ in $L^{2}(\mathbb{R}^{d}\times(0,T))$ as $n\to\infty$, the proof of (15) is complete if we can show that $a_{n}\to 0$ weakly in $L^{2}(\mathbb{R}^{d}\times(0,T))$. We have

[TABLE]

By Lipschitz continuity of $\vec{{\bm{v}}}$, there exists $C_{\vec{{\bm{v}}}}>0$ depending only on $\vec{{\bm{v}}}$ such that $|\nabla\rho_{n}({\bm{y}}-{\bm{x}})\cdot(\vec{{\bm{v}}}({\bm{x}},t)-\vec{{\bm{v}}}({\bm{y}},t))|\leq C_{\vec{{\bm{v}}}}|{\bm{y}}-{\bm{x}}|\,|\nabla\rho_{n}({\bm{y}}-{\bm{x}})|$. Noting that the sequence of functions ${\bm{z}}\mapsto|{\bm{z}}|\,|\nabla\rho_{n}({\bm{z}})|$ is bounded in $L^{1}(\mathbb{R}^{d})$, Young’s inequality for convolution shows that the first term in the right-hand side of (16) is bounded in $L^{2}(\mathbb{R}^{d}\times(0,T))$. The same Young inequality also easily shows that the second term in this right-hand side is also bounded in the same space, which proves that $a_{n}$ itself remains bounded in $L^{2}(\mathbb{R}^{d}\times(0,T))$. The weak convergence of $a_{n}$ therefore only needs to be assessed for smooth functions. Taking $\psi\in C^{\infty}_{c}(\mathbb{R}^{d}\times(0,T))$, we have

[TABLE]

Hence, using the Lipschitz continuity of $\vec{{\bm{v}}}$ and the fact that $\rho_{n}$ is supported in the ball centred at [math] and of radius $1/n$, there exists $C>0$ depending only on $\bar{u},\vec{{\bm{v}}}$ and $\psi$ such that

[TABLE]

Hence $a_{n}$ converges to [math] weakly in $L^{2}(\mathbb{R}^{d}\times(0,T))$, which concludes the proof of (15) and of the lemma. ∎

Corollary 2.3** (Uniqueness).**

Assuming (2), there exists at most one solution $\bar{u}$ to (3).

Proof.

The difference of two solutions to (3) is a solution for the same problem with right-hand-side $f=0$ and initial condition $u_{\rm ini}=0$. The energy estimate (13) shows that this difference is a.e. equal to 0. ∎

3 The gradient discretisation method for the linear advection equation

The gradient discretisation method (GDM) is a general framework for nonconforming approximations of elliptic or parabolic problems, see [12] for a general presentation of the method and of some models and schemes it applies to.

The principle of the GDM is to design a set of discrete elements (space, operators) called a gradient discretisation (GD), which is substituted in the weak formulation of the PDE in lieu of the related continuous elements leading to a discretisation scheme.

Definition 3.1** (Gradient discretisation).**

Let $p\in(1,+\infty)$ be given and let $p^{\prime}\in(1,+\infty)$ with $\frac{1}{p}+\frac{1}{p^{\prime}}=1$. A gradient discretisation ${\cal D}$ is defined by ${\cal D}=(X_{{\cal D}},\Pi_{\cal D},\nabla_{\cal D})$ where:

the set of discrete unknowns $X_{{\cal D}}$ is a finite dimensional vector space on $\mathbb{R}$, 2. 2.

the linear mapping $\Pi_{\cal D}~{}:~{}X_{{\cal D}}\to L^{\max(2,p^{\prime})}(\Omega)$ reconstructs functions, 3. 3.

the linear mapping $\nabla_{\cal D}~{}:~{}X_{{\cal D}}\to L^{\max(2,p)}(\Omega)^{d}$ reconstructs approximations of their gradients, 4. 4.

the quantity $\|\cdot\|_{{\cal D}}:=\|\Pi_{\cal D}\cdot\|_{L^{2}(\Omega)}+\|\nabla_{\cal D}\cdot\|_{L^{p}(\Omega)^{d}}$ defines a norm on $X_{{\cal D}}$.

Remark 3.2**.**

In the above definition, the definition of the norm is not standard in the GDM setting (in the sense of [12, Definition 2.1]), because of the simultaneous use of the $L^{p}$, $L^{p^{\prime}}$ and $L^{2}$ norms.

This notion is extended to evolution problems in the following definition.

Definition 3.3** (Space-time gradient discretisation).**

A family ${\cal D}^{T}=(X_{{\cal D}},\Pi_{\cal D},\nabla_{\cal D},{\mathcal{I}}_{\cal D},(t^{(n)})_{n=0,\ldots,N})$ is a space-time gradient discretisation if

• •

${\cal D}=(X_{{\cal D}},\Pi_{\cal D},\nabla_{\cal D})$* is a gradient discretisation of $\Omega$, in the sense of Definition 3.1,*

• •

${\mathcal{I}}_{\cal D}~{}:~{}L^{2}(\Omega)\to X_{{\cal D}}$* is an interpolation operator,*

• •

$t^{(0)}=0<t^{(1)}\ldots<t^{(N)}=T$.

We then set ${\delta\!t}^{(n+{\frac{1}{2}})}=t^{(n+1)}-t^{(n)}$, for $n=0,\ldots,N-1$, and ${\delta\!t}_{\cal D}=\max_{n=0,\ldots,N-1}{\delta\!t}^{(n+{\frac{1}{2}})}$.

The properties of GDs are assessed through the two following functions $S_{{\cal D}}$ and $W_{{\cal D}}$. The first one measures an interpolation error:

[TABLE]

whilst the second one is a measure of a conformity defect (i.e. the defect in a discrete integration-by-parts formula): letting $W^{1,\infty}_{\bm{n},0}(\Omega)^{d}$ be the set of elements of $W^{1,\infty}(\Omega)^{d}$ with zero normal trace on $\partial\Omega$,

[TABLE]

Let us now define the space size of a GD relative to some regularity spaces. This definition, which holds for both mesh-based and meshless methods, is a measure of the approximation properties of a given GD (this notion is defined in the framework of elliptic problems with homogeneous boundary conditions in [12, Definition 2.22]).

Definition 3.4** (Space size of a GD).**

Let ${\cal D}$ be a gradient discretisation. The space-size of ${\cal D}$ is $h_{\cal D}$ defined by

[TABLE]

Remark 3.5** (Link between $h_{\cal D}$ and the size of the mesh for mesh-based GDs).**

In the case of the mesh-based GDs detailed in [12, Chapters 8-14], $h_{{\cal D}}$ is related to the size of the mesh $h_{\mathcal{M}}$ by $h_{\cal D}\leq Ch_{\mathcal{M}}$ (see, e.g., [12, Remark 2.24]).

Definition 3.6** (Consistent and limit-conforming sequence of space-time gradient discretisation).**

A sequence $({\cal D}^{T}_{m})_{m\in\mathbb{N}}$ of space-time gradient discretisations is said to be consistent and limit-conforming if $h_{{\cal D}_{m}}$, ${\delta\!t}_{{\cal D}_{m}}$ and, for all $w\in L^{2}(\Omega)$, $\|w-\Pi_{{\cal D}_{m}}{\mathcal{I}}_{{\cal D}_{m}}w\|_{L^{2}(\Omega)}$ tend to 0 as $m\to\infty$.

Remark 3.7** (Link with the core properties of a GD in the framework of elliptic or parabolic problems).**

An adaptation of [12, Lemma 2.25] to elliptic problems with homogeneous Neumann boundary conditions yields an equivalence between Definition 3.6 and [12, Definitions 3.4 and 3.5] of consistent and limit-conforming sequences of gradient discretisations, assuming that the sequence $({\cal D}^{T}_{m})_{m\in\mathbb{N}}$ is compact (this holds true for the GDs detailed in [12, Chapters 8-14]).

Given a space–time gradient discretisation ${\cal D}^{T}=(X_{{\cal D}},\Pi_{\cal D},\nabla_{\cal D},{\mathcal{I}}_{\cal D},(t^{(n)})_{n=0,\ldots,N})$ (in the sense of Definition 3.3), we now describe the gradient scheme defined from this GD. For $n=0,\ldots,N-1$ and a given space-time function $g\in L^{1}(\Omega\times(0,T))^{\ell}$ with $\ell=1$, $\ell=d$ or $\ell=d\times d$ ($g$ could be $\Lambda$, $f$, $\vec{{\bm{v}}}$, $q^{I}$ or $q^{P}$), set, for a.e. ${\bm{x}}\in\Omega$ and for all $n=0,\ldots,N-1$,

[TABLE]

Let $\theta\in[\frac{1}{2},1]$ and $\alpha\in(0,p)$. The ($\theta$-implicit) scheme for Problem (3) is defined by replacing the continuous space and operators in (7) with their discrete counterparts given by ${\cal D}$, as follows: find $u=(u^{(n)})_{n=0,\ldots,N}$ such that

[TABLE]

denoting for short

[TABLE]

We introduce the following notations $\Pi^{\theta}_{\cal D}$ and ${\nabla^{\theta}_{\!{\cal D}}}$ for reconstructed space-time functions: given $v=(v_{n})_{n=0,\ldots,N}$ in $X_{\cal D}^{N+1}$, we set

[TABLE]

We extend these definitions to $t=0$ by setting $\Pi_{\cal D}v({\bm{x}},0)=\Pi_{\cal D}v_{0}$ and ${\nabla_{\!{\cal D}}}v({\bm{x}},0)={\nabla_{\!{\cal D}}}v_{0}$.

4 Convergence analysis

Our main convergence result is stated in the following theorem. We recall (see [11, Definition 2.11]) that a sequence $(v_{n})_{n\in\mathbb{N}}$ of bounded functions $[0,T]\to L^{2}(\Omega)$ is said to converge uniformly on $[0,T]$ weakly in $L^{2}(\Omega)$ towards a function $v$ if for all $\varphi\in L^{2}(\Omega)$, the sequence of functions $t\in[0,T]\mapsto\langle v_{n}(t),\varphi\rangle_{L^{2}(\Omega)}\in\mathbb{R}$ converges uniformly on $[0,T]$ towards the function $t\mapsto\langle v(t),\varphi\rangle_{L^{2}(\Omega)}$.

Theorem 4.1** (Convergence of the GDM).**

Assuming (2), let $({\cal D}^{T}_{m})_{m\in\mathbb{N}}$ be a consistent and limit-conforming sequence of space-time gradient discretisations in the sense of Definition 3.6. Let $\theta\in[{\frac{1}{2}},1]$, $p\in(1,+\infty)$ and $\alpha\in(0,p)$ be given. Then, for any $m\in\mathbb{N}$, there exists a unique $u_{m}$ solution to Scheme (21) with ${\cal D}^{T}={\cal D}_{m}^{T}$.

Moreover, as $m\to\infty$, $\Pi_{{\cal D}_{m}}^{\theta}u_{m}$ converges in $L^{2}(\Omega\times(0,T))$, and uniformly on $[0,T]$ weakly in $L^{2}(\Omega)$, to the unique solution $\bar{u}$ of Problem (3).

Remark 4.2** (Theoretical order of convergence).**

Assuming sufficient smoothness of the continuous solution and the velocity field, and letting $p=2$, it seems possible to derive a theoretical error estimate in $L^{\infty}(0,T;L^{2}(\Omega))$ norm with order $h_{\cal D}^{\min(\alpha,2-\alpha)}$. This provides a maximal order 1 if $\alpha=1$. However, in the numerical tests with a regular solution (see Section 5.3), much better numerical orders of convergence are obtained, even letting $\alpha=2$. The question of the theoretical derivation of these better rates remains open.

The uniqueness component of this theorem is the most straightforward part, and the purpose of the following lemma.

Lemma 4.3** (Uniqueness of a discrete solution).**

Assuming (2), let ${\cal D}^{T}=(X_{{\cal D}},\Pi_{\cal D},\nabla_{\cal D},{\mathcal{I}}_{\cal D},(t^{(n)})_{n=0,\ldots,N})$ be a space-time gradient discretisation in the sense of Definition 3.3. Let $\theta\in[{\frac{1}{2}},1]$, $p\in(1,+\infty)$ and $\alpha\in(0,p)$ be given. Then there exists at most one solution to Scheme (21).

Proof.

The scheme defines exactly one approximation $u^{(0)}$. Let us assume that, for a given $n\in\mathbb{N}$ and for a given $u^{(n)}$, there exist two solutions $u^{(n+1)}$ and $\widehat{u}^{(n+1)}$ to Scheme (21). Let us create the difference of the two equations (21), and let us choose $v=u^{(n+\theta)}-\widehat{u}^{(n+\theta)}=\theta(u^{(n+1)}-\widehat{u}^{(n+1)})$ in the resulting equation. We obtain

[TABLE]

It is classical (see for instance [12, Lemma 2.40] or [2, Lemma 2.1]), that

[TABLE]

Applying this inequality in (24) with $\chi=(\Lambda^{(n+{\frac{1}{2}})})^{1/2}{\nabla_{\!{\cal D}}}u^{(n+\theta)}$ and $\xi=(\Lambda^{(n+{\frac{1}{2}})})^{1/2}{\nabla_{\!{\cal D}}}\widehat{u}^{(n+\theta)}$ (in which the left hand side is therefore the sum of non-negative terms), we get that ${\nabla_{\!{\cal D}}}u^{(n+\theta)}={\nabla_{\!{\cal D}}}\widehat{u}^{(n+\theta)}$ a.e., and therefore ${\nabla_{\!{\cal D}}}u^{(n+1)}={\nabla_{\!{\cal D}}}\widehat{u}^{(n+1)}$ as well as $\Pi_{\cal D}u^{(n+1)}=\Pi_{\cal D}\widehat{u}^{(n+1)}$. Hence, thanks to the property of the norm assumed in Definition 3.1, $u^{(n+1)}=\widehat{u}^{(n+1)}$, which concludes the proof of uniqueness by induction.

∎

The proof of Theorem 4.1 hinges on a priori estimates stated in the following lemma.

Lemma 4.4** ($L^{\infty}(0,T;L^{2}(\Omega))$ and discrete $L^{2}(0,T;H^{1}_{0}(\Omega))$ estimates, existence of a discrete solution).**

Assuming (2), let ${\cal D}^{T}=(X_{{\cal D}},\Pi_{\cal D},\nabla_{\cal D},{\mathcal{I}}_{\cal D},(t^{(n)})_{n=0,\ldots,N})$ be a space-time gradient discretisation in the sense of Definition 3.3. Let $\theta\in[{\frac{1}{2}},1]$, $p\in(1,+\infty)$ and $\alpha\in(0,p)$ be given. Then there exists one and only one solution to Scheme (21). Moreover, this solution satisfies, for all $k=1,\ldots,N$,

[TABLE]

and there exists $\Cl[ctrcst]{estimld}>0$, depending only on $C_{\rm ini}\geq\|u_{\rm ini}-\Pi_{\cal D}{\mathcal{I}}_{\cal D}u_{\rm ini}\|_{L^{2}(\Omega)}$, $\underline{\lambda}$, $f$ and $q^{I}$ such that

[TABLE]

and

[TABLE]

Remark 4.5** (Weak $BV$ estimate).**

The estimate (27) is the adaptation in the GDM framework of the classical weak $BV$ estimate used for finite volumes see [5] for the seminal paper and [16, chapters 5 & 6] for more general results. This estimate is used in two occasions: first to pass to the limit in the skew-symmetric term, and second to show that the stabilisation term vanishes at the limit.

Proof.

Before establishing the existence of at least one discrete solution to Scheme (21), let us first prove that any solution to this scheme satisfies (25)–(27). We first notice that for all $a,b\in\mathbb{R}$,

[TABLE]

Hence, letting $v={\delta\!t}^{(n+{\frac{1}{2}})}u^{(n+\theta)}$ in (21) and applying the estimate above with $a=\Pi_{\cal D}u^{(n+1)}$ and $b=\Pi_{\cal D}u^{(n)}$, we obtain

[TABLE]

Taking $k=1,\ldots,N$ and summing this inequality over $n=0,\ldots,k-1$ proves (25).

The Young inequality and the property $0\leq q^{I}_{\cal D}\leq q^{I}_{\cal D}+q^{P}_{\cal D}$ yield

[TABLE]

Plugging this into (25) leads to

[TABLE]

where we have used the Jensen inequality to bound the $L^{2}$-norm of $f_{\cal D}$ by the $L^{2}$-norm of $f$. Estimate (27) directly follows from (28) with $k=N$. Estimate (26) is also a consequence of (28), once we notice that $\Pi_{\cal D}^{\theta}u({\bm{x}},t)=\theta\Pi_{\cal D}u^{(n+1)}({\bm{x}})+(1-\theta)\Pi_{\cal D}u^{(n)}({\bm{x}})$ for a.e. ${\bm{x}}\in\Omega$, all $t\in(t^{(n)},t^{(n+1)})$ and all $n=0,\ldots,N-1$.

We can now prove the existence of a solution to Scheme (21) (the uniqueness is proved in Lemma 4.3). If $p=2$ then, at each time step, (21) describes a linear square system on $u^{(n+\theta)}$ (after substituting $u^{(n+1)}=u^{(n)}+\theta^{-1}(u^{(n+\theta)}-u^{(n)})$). The estimates (26) and (27) show that any solution $u^{(n+\theta)}$ to this system satisfies a priori bounds. The kernel of the matrix of this linear system is therefore reduced to $\{0\}$, and the matrix is invertible, which establishes the existence of a unique solution $u^{(n+\theta)}$ (and thus of $u^{(n+1)}$) to the system at time step $n+1$.

If $p\neq 2$ we use the topological degree [8]. Let us assume the existence of $u^{(n)}$. Let us substitute the term $|{\nabla_{\!{\cal D}}}u^{(n+\theta)}|_{\Lambda}^{p-2}{\nabla_{\!{\cal D}}}u^{(n+\theta)}$ of the scheme by $\nu|{\nabla_{\!{\cal D}}}u^{(n+\theta)}|_{\Lambda}^{p-2}{\nabla_{\!{\cal D}}}u^{(n+\theta)}+(1-\nu){\nabla_{\!{\cal D}}}u^{(n+\theta)}$ for $\nu\in[0,1]$. It is clear that the above estimates still hold (again after substituting $u^{(n+1)}=u^{(n)}+\theta^{-1}(u^{(n+\theta)}-u^{(n)})$) so that $\|\Pi_{\cal D}u^{(n+1)}\|_{L^{2}(\Omega)}$ and $\nu\|{\nabla_{\!{\cal D}}}u^{(n+1)}\|_{L^{p}(\Omega)}^{p}+(1-\nu)h_{\cal D}^{\alpha}\|{\nabla_{\!{\cal D}}}u^{(n+1)}\|_{L^{2}(\Omega)}^{2}$ remain bounded independently of $\nu$. We infer from this latter estimate a bound on $\|u^{(n+1)}\|_{{\cal D}}$ that is uniform with respect to $\nu$. Hence, all solutions to the scheme with the above substitution remain bounded independently of $\nu$. This shows that, on a large enough ball, the topological degree of the non-linear mapping defining the scheme is independent of $\nu$. For $\nu=0$ this mapping is linear and the arguments developed in the case $p=2$ show that its topological degree is non-zero. The degree for the original scheme (corresponding to $\nu=1$) is therefore also non-zero, proving that this scheme has at least one solution. ∎

We can now prove our convergence results, starting with the uniform-in-time weak-in-space convergence.

Proof of Theorem 4.1: uniform-in-time weak-in-space convergence.

Owing to (26) there is $\bar{u}\in L^{2}(\Omega\times(0,T))$ and a subsequence of $({\cal D}_{m})_{m\in\mathbb{N}}$ such that $\Pi^{\theta}_{{\cal D}_{m}}u_{m}$ converges to $\bar{u}$ in $L^{\infty}(0,T;L^{2}(\Omega))$ weak-$\star$ as $m\to\infty$. Let $m\in\mathbb{N}$, and let us denote ${\cal D}={\cal D}_{m}$ (belonging to the above subsequence); we drop some indices $m$ to simplify the notations.

Let $\varphi\in C^{\infty}_{c}([0,T))$ and $w\in C^{\infty}_{c}(\mathbb{R}^{d})$, and let $P_{m}w\in X_{\cal D}$ that realises the minimum in $S_{{\cal D}}(w)$. We denote by $P_{m}\varphi:(0,T)\to\mathbb{R}$ the function equal to $\varphi^{(n+1-\theta)}:=\theta\varphi(t^{(n)})+(1-\theta)\varphi(t^{(n+1)})$, on $(t^{(n)},t^{(n+1)})$, for all $n=0,\ldots,N-1$.

For $n=0,\ldots,N-1$ and $t\in(t^{(n)},t^{(n+1)})$, let $u^{(\theta)}_{m}(t)=\theta u_{m}^{(n+1)}+(1-\theta)u_{m}^{(n)}\in X_{\cal D}$, and notice that ${\Pi^{\theta}_{\!{\cal D}}}u_{m}(t)=\Pi_{\cal D}u_{m}^{(\theta)}(t)$ and ${\nabla^{\theta}_{\!{\cal D}}}u_{m}(t)={\nabla_{\!{\cal D}}}u_{m}^{(\theta)}(t)$. By definition (18) of $W_{\cal D}$ and (19) of $h_{\cal D}$, since $w\vec{{\bm{v}}}(t)\in W^{1,\infty}_{\bm{n},0}(\Omega)^{d}$ we have, for a.e. $t\in(0,T)$, recalling the definition (19) of $h_{\cal D}$,

[TABLE]

Thanks to (26)–(27), there is $\Cr{estimldd}$ depending only on $\Cr{estimld}$ and $T$ such that

[TABLE]

This right-hand side tends to [math] as $m\to\infty$ (remember that $\alpha<p$) and thus, since $P_{m}\varphi$ is bounded in $L^{\infty}(0,T)$,

[TABLE]

By strong convergence of $P_{m}\varphi$ to $\varphi$ in $L^{2}(0,T)$ and weak convergence of ${\Pi^{\theta}_{\!{\cal D}}}u_{m}$ to $\bar{u}$ in $L^{2}(\Omega\times(0,T))$ we infer

[TABLE]

A Cauchy–Schwarz inequality yields

[TABLE]

and, by definition of $P_{m}\varphi$, $h_{\cal D}$ and $P_{m}w$,

[TABLE]

Therefore, using (27) again,

[TABLE]

We take ${\delta\!t}^{(n+{\frac{1}{2}})}\varphi(t^{(n)})P_{m}w$ as test function in (21) and sum the resulting equation over $n=0,\ldots,N-1$. This gives

[TABLE]

with

[TABLE]

The summation-by-parts formula [12, Eq. (D.17)] reads

[TABLE]

Using this relation to transform, in the sum appearing in $\Cr{1m}^{(m)}$, the term ${\delta\!t}^{(n+{\frac{1}{2}})}\varphi^{(n+1-\theta)}\delta_{\cal D}^{(n+{\frac{1}{2}})}u_{m}$ into $(\varphi(t^{(n)})-\varphi(t^{(n+1)}))\Pi_{\cal D}u_{m}^{(n+\theta)}$, we see that

[TABLE]

and so, since ${\Pi^{\theta}_{\!{\cal D}}}u_{m}\to\bar{u}$ weakly in $L^{2}(\Omega\times(0,T))$, ${\Pi_{\!{\cal D}}}P_{m}w\to w$ strongly in $L^{2}(\Omega)$, and $\Pi_{\cal D}u_{m}^{(0)}=\Pi_{{\cal D}_{m}}{\mathcal{I}}_{{\cal D}_{m}}u_{\rm ini}\to u_{\rm ini}$ in $L^{2}(\Omega)$,

[TABLE]

Noticing that

[TABLE]

the relation (29) yields

[TABLE]

Moreover, since

[TABLE]

${\Pi^{\theta}_{\!{\cal D}}}u$ weakly converges to $u$ in $L^{2}(\Omega\times(0,T))$, ${\Pi_{\!{\cal D}}}(P_{m}w)$ strongly converges in $L^{2}(\Omega)$, and ${\nabla_{\!{\cal D}}}(P_{m}w)$ strongly converges to $\nabla w$ in $L^{2}(\Omega)^{d}$, we have

[TABLE]

The Hölder inequality and (27) show that

[TABLE]

The boundedness of ${\nabla_{\!{\cal D}}}(P_{m}w)$ in $L^{p}(\Omega)$ (since this sequence converges in this space) and $\alpha>0$ then yield $\lim_{m\to\infty}\Cr{grad}^{(m)}=0$.

Finally, using (31) again,

[TABLE]

Passing to the limit $m\to\infty$ in (30) shows that $\bar{u}$ satisfies (3) for any test function of the form $\varphi(t)w({\bm{x}})$, and thus for sums of such test functions. Since the set $\mathcal{T}=\{\sum_{i=1}^{q}\varphi_{i}(t)w_{i}({\bm{x}})\,:\,q\in\mathbb{N},\varphi_{i}\in C^{\infty}_{c}[0,T),w_{i}\in C^{\infty}_{c}(\mathbb{R}^{d})\}$ is dense in the set of the restrictions to $\overline{\Omega}\times[0,T)$ of the elements of $C^{\infty}_{c}(\mathbb{R}^{d}\times[0,T))$, we conclude that $\bar{u}$ is a solution of (3).

It now suffices to prove the uniform-in-time weak-$L^{2}(\Omega)$ convergence of $\Pi_{{\cal D}_{m}}^{\theta}u_{m}$ to $\bar{u}$. Let $w\in C^{\infty}_{c}(\mathbb{R}^{d})$ and define $P_{m}w$ as before. For $0\leq s\leq t$, writing $\Pi_{{\cal D}_{m}}^{\theta}u_{m}({\bm{x}},t)-\Pi_{{\cal D}_{m}}^{\theta}u_{m}({\bm{x}},s)$ as the sum of its jumps at each $t^{(n)}\in(s,t)$ (see [12, Proof of Theorem 4.19] for details), Scheme (21) and the estimates in Lemma 4.4 give the existence of $\Cl[ctrcst]{cquatre}$, depending only on the data introduced in 2, such that

[TABLE]

Hence, introducing $\pm w$, and using (26) again,

[TABLE]

Using $\sqrt{t-s+2{\delta\!t}_{m}}\leq\sqrt{t-s}+\sqrt{2{\delta\!t}_{m}}$, we get

[TABLE]

with $g(a,b)=\sqrt{a}\Cl[ctrcst]{ccinq}+b$, $\Cr{ccinq}=\Cr{cquatre}\sup_{m}\|P_{m}w\|_{{\cal D}_{m}}$ and $h^{\varphi}_{m}=(2{\delta\!t}_{m})^{1/2}\Cr{ccinq}+2C_{1}\|w-\Pi_{{\cal D}_{m}}(P_{m}w)\|_{L^{2}(\Omega)}$. We then may apply [12, Theorem C.11] or [11, Theorem 6.2] to deduce that $\Pi_{{\cal D}_{m}}^{\theta}u_{m}$ weakly tends to $u$ in $L^{2}(\Omega)$ uniformly on $[0,T]$. ∎

Proof of Theorem 4.1: strong convergence.

The proof makes use of the continuous energy estimate (13) and a discrete version thereof, in a similar way as in the proof of [11, Theorem 2.16]. Let us first establish this discrete energy estimate. We remark that for all $a,b,c,d\in\mathbb{R}$,

[TABLE]

Setting $t^{(n+1-\theta)}=\theta t^{(n)}+(1-\theta)t^{(n+1)}$, letting $v={\delta\!t}^{(n+{\frac{1}{2}})}u^{(n+\theta)}(T-t^{(n+1-\theta)})$ in (21), applying the above relation with $a=\Pi_{\cal D}u^{(n+1)}$, $b=\Pi_{\cal D}u^{(n)}$, $\alpha=2\theta-1$, $c=t^{(n)}$, $d=t^{(n+1)}$, and dropping the last addend (which is positive), we obtain

[TABLE]

Summing the obtained inequality on $n=0,\ldots,N-1$, and denoting by $t_{\cal D}$ the function equal to $t^{(n+1-\theta)}$ for all $t\in(t^{(n)},t^{(n+1)})$, we get

[TABLE]

Taking the superior limit as $m\to\infty$ of the above inequality for ${\cal D}={\cal D}_{m}$, we get

[TABLE]

We then use (13) to substitute the right-hand side of this inequality and find

[TABLE]

Developing the square $(\Pi_{{\cal D}_{m}}^{\theta}u_{m}-\bar{u})^{2}$ we have

[TABLE]

The limit of the second (resp. third) term in the right-hand side is obtained by weak/strong (resp. strong) convergence:

[TABLE]

and

[TABLE]

Hence, using (35) to deal with the first term in the right-hand side of (36) we find

[TABLE]

and therefore, since $1+(T-t_{{\cal D}_{m}})\left[q^{I}_{{\cal D}_{m}}+q^{P}_{{\cal D}_{m}}\right]\geq 1$,

[TABLE]

which concludes the proof of the convergence of $\Pi_{{\cal D}_{m}}^{\theta}u_{m}$ to $\bar{u}$ in $L^{2}(\Omega\times(0,T))$.

∎

5 Numerical results

Let $\Omega=(0,1)^{2}$ and consider meshes $\mathfrak{T}=({\mathcal{M}},{\mathcal{F}},{\cal P},\mathcal{V})$ as per [12, Definition 7.2]: ${\mathcal{M}}$ is the set of polygonal/polyhedral cells $K$, ${\mathcal{F}}$ is the set of faces $\sigma$, ${\cal P}$ is a set of points $({\bm{x}}_{K})_{K\in{\mathcal{M}}}$ with $K$ star-shaped with respect to ${\bm{x}}_{K}$ for all $K\in{\mathcal{M}}$, and $\mathcal{V}$ is the set of vertices $\mathsf{v}$. Let us define two test cases.

Case 1.

This test case is divergence free (it corresponds to the pure transport of a tracer). We choose $T=5$, $u_{\rm ini}({\bm{x}})=1$ if ${\bm{x}}=(x_{1},x_{2})\in(0.1,0.4)\times(0.1,0.4)$ and $u_{\rm ini}(x)=0$ elsewhere, $q^{I}=q^{P}=0$ and $\vec{{\bm{v}}}$ is given by

[TABLE]

Case 2.

This test case includes source terms. We choose $T=1$, $u_{\rm ini}\equiv 0$, $\vec{{\bm{v}}}$ is given by

[TABLE]

$q^{I}({\bm{x}},t)=\max({\rm div}(\vec{{\bm{v}}})({\bm{x}}),0)=\max(2-2(x_{1}+x_{2}),0)$ and $f({\bm{x}},t)=1$, $q^{P}({\bm{x}},t)=\max(-{\rm div}(\vec{{\bm{v}}})({\bm{x}}),0)=\max(2(x_{1}+x_{2})-2,0)$. Then the solution $\bar{u}({\bm{x}},t)$ to (1) can be analytically calculated. Find first ${\bm{X}}(s)=(X_{1}(s),X_{2}(s))$ by solving the differential equation ${\bm{X}}^{\prime}(s)=\vec{{\bm{v}}}({\bm{X}}(s))$, letting ${\bm{X}}(0)=\widehat{{\bm{x}}}=(\widehat{x}_{1},\widehat{x}_{2})$ for any $\widehat{{\bm{x}}}\in\Omega$ (this is easy owing to the expression of $\vec{{\bm{v}}}$); set then $\bar{v}(s)=\bar{u}({\bm{X}}(s),s)$, which leads to $\bar{v}^{\prime}(s)=(f({\bm{X}}(s))-\bar{v}(s))q^{I}({\bm{X}}(s))$ with $\bar{v}({\bm{X}}(0))=0$. This requires to compute $\widehat{s}$ such that $X_{1}(\widehat{s})+X_{2}(\widehat{s})=1$, since $q^{I}({\bm{X}}(s))=0$ for $s\geq\widehat{s}$. Finally, we get that $\bar{u}({\bm{x}},t)=\bar{v}(t)$, when $\widehat{{\bm{x}}}$ is chosen such that ${\bm{X}}(t)={\bm{x}}$, and we have $\bar{u}({\bm{x}},t)=0$ if $\widehat{x}_{1}+\widehat{x}_{2}>1$. Denoting by $\alpha:=\min(1,e^{\widehat{s}-t})$, this leads to the following expressions.

[TABLE]

5.1 Case 1, different schemes with $p=2$

We apply Scheme (21) with three different gradient discretisations, corresponding respectively to the mass-lumped conforming $\mathbb{P}_{1}$ finite element method (or CVFE method, see [18] for the seminal paper and [12, Chapter 8] for the study in the GDM framework), to the mass-lumped non-conforming $\mathbb{P}_{1}$ (MLNC–$\mathbb{P}_{1}$ for short) finite element method [12, Chapter 9], and to (a variant of) the Hybrid Finite Volume method (HFV), a member of the family of Hybrid Mimetic Mixed methods [12, chapter 13]. For the sake of completeness we briefly recall the definition of these gradient discretisations.