# A fully local hybridised second-order accurate scheme for advection-
diffusion equations

Hanz Martin Cheng  Department of Mathematics and Computer Science, Eindhoven
University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands.
<EMAIL_ADDRESS>
###### Abstract.

In this paper, we present a fully local second-order upwind scheme, applicable
on generic meshes. This is done by hybridisation, which is achieved by
introducing unknowns on each edge of the mesh. By doing so, fluxes only depend
on values associated to a single cell, and thus, this scheme can easily be
applied even on cells near the boundary of the domain. Another advantage of
hybridised schemes is that static condensation can be employed, leading to a
very efficient implementation. A convergence analysis, which also covers a
flux-limited TVD variant of the scheme, is then presented. Numerical results
are also given in order to compare this with a hybridised first-order upwind
scheme and a classical cell-centered second-order upwind type scheme.

##  1\. Introduction

In this paper, we study a family of finite volume methods for stationary
advection-diffusion equations. We start by presenting the choice of
discretisation for the diffusive fluxes. This will be done via the hybrid
mimetic mixed (HMM) method [6], which is equivalent to the SUSHI method [9].
For the advective fluxes, we propose a fully local second-order scheme. To
motivate the problem, we start with a revision of the first-order upwind
scheme, which is the easiest to implement. One of the main disadvantages of
the first-order upwind scheme is that it may easily introduce too much
numerical diffusion into the solution of the system [11]. High-order schemes
have been proposed in order to help mitigate the introduction of too much
numerical diffusion. However, in order to achieve a high-order discretisation,
more degrees of freedom (DOFs) are needed, which translates into a higher
computational cost. In this work, we focus on second-order schemes. Classical
second-order upwind schemes on Cartesian meshes [16] involve a 9-point
stencil; hence, some interpolation techniques or introduction of ghost cells
are needed in order to apply them on cells near the boundary of the domain.

The novelty of this work is the introduction of a hybridised fully local
second-order scheme, inspired by the ideas in [3]. To do so, in addition to
unknowns at the cell center, we introduce one additional unknown to each cell
face. This results to more unknowns in the system: now having $\sharp$ of
cells $+$ $\sharp$ of faces (edges) unknowns, as compared to $\sharp$ of cells
for cell-centered schemes. However, this allows us to reduce the dependence of
the fluxes on values from neighboring cells. In the case of square cells (see
Figure 1), we only need one cell value, and four interface values. Since we
remove the direct dependence of the fluxes on the values from neighboring
cells, this allows us to directly apply this hybridised second-order scheme
even near the boundaries of the domain. Moreover, static condensation can be
employed to make the implementation more efficient. Following this, we then
present some convergence results, which cover, for the advective component,
the classical second-order upwind type schemes, and also include some
nonlinear total variation diminishing (TVD) methods [4, 13].

Figure 1. Values needed for computing a hybridised second-order flux for edge
$\sigma$.

The paper is organised into the following sections: We start by considering a
stationary advection-diffusion equation, and present it in its finite volume
form. Afterwards, we provide a discretisation of the diffusive fluxes, which
will be done via the HMM method. We then give a short review of the cell-
centered first and second-order upwind schemes on Cartesian meshes and give
their extension onto generic polygonal meshes. Following this, we propose a
hybridised fully local second-order scheme, applicable on generic meshes.
Here, we show that by introducing additional unknowns along the faces (edges)
of each cell, we can create a stencil that does not depend on values from
neighboring cells. Convergence results will then be presented. Numerical tests
will then be performed to illustrate the accuracy of the scheme, and also to
compare this with a hybridised first-order upwind scheme and a classical cell-
centered second-order upwind type scheme.

##  2\. The model problem

Consider the following stationary, linear scalar advection-diffusion equation
with Dirichlet boundary conditions on a polytopal domain $\Omega$: Find $c\in
H^{1}(\Omega)$ such that

$\displaystyle\nabla\cdot(-\Lambda\nabla c+c\mathbf{V})$
$\displaystyle=f\qquad\mathrm{in}\quad\Omega,$ (1) $\displaystyle c$
$\displaystyle=g\qquad\mathrm{in}\quad\partial\Omega.$

Here, $\Lambda$ is a diffusion tensor, $\mathbf{V}$ is a velocity field, and
$f$ is a source term. We begin by stating the assumptions on the data:

  1. (A1)

$\Lambda$ is a measurable function from $\Omega$ to the set of $d\times d$
symmetric positive definite matrices, and there exists
$\underline{\lambda},\overline{\lambda}>0$ such that, for a.e. $x\in\Omega$,
the eigenvalues of $\Lambda(x)$ are in
$[\underline{\lambda},\overline{\lambda}]$;

  2. (A2)

$\mathbf{V}\in C^{1}(\overline{\Omega})^{d}$ with ${\mathop{\rm
div}}(\mathbf{V})\geq 0$;

  3. (A3)

$f\in L^{2}(\Omega)$.

###  2.1. Finite volume discretisation of the advection-diffusion equation

We now write a finite volume discretisation of the model (1). First, we define
a mesh in the simplest intuitive way: a partition of $\Omega$ into polygonal
(in 2D) or polyhedral (in 3D) sets. Following the notations in [5, Definition
7.2], we denote $\mathcal{T}=({\mathcal{M}};{\mathcal{E}})$ to be the set of
cells $K$ and faces (edges in 2D) $\sigma$ of our mesh, respectively. For each
cell $K\in{\mathcal{M}}$, we denote by $|K|$ its $d$ dimensional measure,
$\mathrm{diam}(K)$ its diameter, ${\mathcal{E}}_{K}\subset{\mathcal{E}}$ the
set of faces (edges) of cell $K$, and ${\boldsymbol{x}}_{K}$ the cell center
of gravity. The collection of faces is a disjoint union of two sets,
${\mathcal{E}}={\mathcal{E}}_{\mathrm{int}}\cup{\mathcal{E}}_{\mathrm{ext}}$,
where ${\mathcal{E}}_{\mathrm{int}}$ and ${\mathcal{E}}_{\mathrm{ext}}$ denote
the set of interior and exterior faces, respectively. For each interior face
$\sigma\in{\mathcal{E}}_{K}$, we denote by $K^{\prime}$ the cell that shares
the face $\sigma$ with $K$. Also, for $\sigma\in{\mathcal{E}}$, we denote by
$|\sigma|$ its $d-1$ dimensional measure, ${\boldsymbol{x}}_{\sigma}$ its
center, and $\mathbf{n}_{K,\sigma}$ its normal direction pointing out of $K$.
We also denote by $d_{K,\sigma}$ the orthogonal distance from
${\boldsymbol{x}}_{K}$ to $\sigma$ (see Figure 2).

Denoting by ${\mathcal{M}}_{h}$ the mesh such that,
$\max_{K\in{\mathcal{M}}_{h}}(\mathrm{diam}(K))=h$, the analysis performed in
Section 5 will require the following assumptions on the mesh.

  1. (MR1)

Every cell $K\in{\mathcal{M}}_{h}$ is star-shaped with respect to
${\boldsymbol{x}}_{K}$.

  2. (MR2)

The mesh regularity parameter, which is defined as

$\mathrm{regul}({\mathcal{M}}_{h}):=\max\bigg{(}\max_{\sigma\in{\mathcal{E}}_{h,\mathrm{int}}}\frac{d_{K,\sigma}}{d_{K^{\prime},\sigma}},\max_{\sigma\in{\mathcal{E}}_{h,\mathrm{ext}}}\frac{\mathrm{diam}(K)}{d_{K,\sigma}},\max_{K\in{\mathcal{M}}_{h}}\mathrm{card}({\mathcal{E}}_{K})\bigg{)},$

where $\mathrm{card}({\mathcal{E}}_{K})$ denotes the number of edges of a cell
$K$, is uniformly bounded as $h\rightarrow 0$.

Figure 2. Notations in a mesh cell.

We now denote the unknowns of the system, one on each cell and one on each
edge by

$X_{\mathcal{D}}:=\\{\big{(}(w_{K})_{K\in{\mathcal{M}}},(w_{\sigma})_{\sigma\in{\mathcal{E}}},w_{K},w_{\sigma}\in\mathbb{R}\big{)}\\}.$

Now, we take the integral of (1) over a cell $K$, and use Gauss’ Theorem to
obtain

$\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}\int_{\sigma}(-\Lambda\nabla
c+c\mathbf{V})\cdot\mathbf{n}_{K,\sigma}\,\,ds=\int_{K}f\,\,d{\boldsymbol{x}}.$

Following the ideas in [3], we write
$F_{K,\sigma}^{D}\approx\frac{1}{|\sigma|}\int_{\sigma}\Lambda\nabla
c\cdot\mathbf{n}_{K,\sigma}\,\,ds$ and
$F_{K,\sigma}^{A}\approx\frac{1}{|\sigma|}\int_{\sigma}c\mathbf{V}\cdot\mathbf{n}_{K,\sigma}\,\,ds$
approximations of the average diffusive and advective fluxes along $\sigma$,
respectively. The discrete equation for flux balance in each cell $K$ reads

$\displaystyle\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|(F_{K,\sigma}^{D}+F_{K,\sigma}^{A})=|K|f_{K},$

where $f_{K}$ is the average value of the source term at cell $K$. Here, we
see that the main ingredients for obtaining an approximation for the solution
$c$ is the definition of the diffusive and advective fluxes $F_{K,\sigma}^{D}$
and $F_{K,\sigma}^{A}$.

##  3\. Diffusive fluxes

In this section, we discuss the discretisation of the diffusive fluxes. This
will be done via the hybrid-mimetic-mixed (HMM) method [6], for which the
diffusive fluxes are based on the bilinear form
$a(u,w):=\int_{\Omega}\Lambda\nabla u\cdot\nabla w\,\,d{\boldsymbol{x}}$,
which stems from the weak formulation of the advection-diffusion equation (1).
In particular, we define the diffusive fluxes such that for $c\in
X_{\mathcal{D}}$, we have for all $K\in{\mathcal{M}},v\in X_{\mathcal{D}}$,

$\sum_{\sigma\in{\mathcal{E}}_{K}}|\sigma|F_{K,\sigma}^{D}(v_{K}-v_{\sigma})=\int_{K}\Lambda_{K}\nabla_{{\mathcal{D}}}c\cdot\nabla_{{\mathcal{D}}}v\,\,d{\boldsymbol{x}},$
(2)

where $\Lambda_{K}$ is an approximation of $\Lambda$ at cell $K$, and
$\nabla_{{\mathcal{D}}}$ is a _stabilised discrete gradient_. The discrete
gradient is stabilised in the sense that it consists of a consistent term,
added to a stabilisation term. The consistent term is linearly exact, and is
defined via a cellwise discrete gradient $\overline{\nabla}_{\mathcal{D}}$,
such that for all $q\in X_{\mathcal{D}}$ and for all $K\in{\mathcal{M}}$,
$(\overline{\nabla}_{\mathcal{D}}q)_{|K}=\overline{\nabla}_{K}q$, where

$\overline{\nabla}_{K}q:=\frac{1}{|K|}\sum_{\sigma\in{\mathcal{E}}_{K}}|\sigma|(q_{\sigma}-q_{K})\mathbf{n}_{K,\sigma},\
$ (3)

and $\mathbf{n}_{K,\sigma}$ is the unit outward normal vector of $\sigma$. In
order to ensure the coercivity of the diffusive flux, a stabilisation term
needs to be added to the discrete gradient (3). For the HMM method, the
stabilisation term is defined such that for all $q\in X_{\mathcal{D}}$ and for
all $K\in{\mathcal{M}}$, $S_{K}:X_{\mathcal{D}}\rightarrow L^{2}(K)$ is given
by

$S_{K}(q):=\sum_{\sigma\in{\mathcal{E}}_{K}}S_{K,\sigma}\mathbbm{1}_{D_{K,\sigma}},$
(4)

where $(D_{K,\sigma})_{\sigma\in{\mathcal{E}}_{K}}$ are convex hulls of
$\sigma$ and ${\boldsymbol{x}}_{K}$ (see Figure 3), and

$S_{K,\sigma}:=\frac{\sqrt{d}}{d_{K,\sigma}}[q_{\sigma}-q_{K}-\overline{\nabla}_{K}q\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K})]\mathbf{n}_{K,\sigma}.$
(5)

Here, $d_{K,\sigma}$ is the orthogonal distance from ${\boldsymbol{x}}_{K}$ to
$\sigma$, and ${\boldsymbol{x}}_{\sigma}$ is the center of the edge $\sigma$.

Figure 3. Notations on cell $K$.

The stabilised discrete gradient is then defined such that for all $q\in
X_{\mathcal{D}}$ and for all $K\in{\mathcal{M}}$

$(\nabla_{{\mathcal{D}}}q)_{|K}:=\overline{\nabla}_{K}q+S_{K}(q).$ (6)

One important property of the stabilisation term, which will be needed to
establish the coercivity of the scheme, is the following orthogonality
condition.

######  Lemma 3.1 (Orthogonality of the stabilisation term.).

Let $q\in X_{\mathcal{D}}$ and let $K\in{\mathcal{M}}$. Then, the
stabilisation $S_{K}:X_{\mathcal{D}}\rightarrow L^{2}(K)$ defined as in (4)
satisfies the following orthogonality condition.

$\int_{K}S_{K}q\cdot\phi\,\,d{\boldsymbol{x}}=0\quad\forall\phi\in\mathbb{R}^{d}.$
(7)

A proof of Lemma 3.1 can be found in [5, Chapter 13], but for completeness, we
present an alternative proof below.

Proof.  Since $S_{K}q$ is piecewise constant with value $S_{K,\sigma}$ on each
convex hull $D_{K,\sigma}$, we can write

$\displaystyle\int_{K}S_{K}q\cdot\phi\,\,d{\boldsymbol{x}}$
$\displaystyle=\sum_{\sigma\in{\mathcal{E}}_{K}}\int_{D_{K,\sigma}}S_{K,\sigma}\cdot\phi\,\,d{\boldsymbol{x}}$
$\displaystyle=\sum_{\sigma\in{\mathcal{E}}_{K}}|D_{K,\sigma}|S_{K,\sigma}\cdot\phi.$

Now, using the definition of $S_{K,\sigma}$ in (5) and the geometric relation

$|D_{K,\sigma}|=\frac{1}{d}|\sigma|d_{K,\sigma},$ (8)

we get

$\displaystyle\sum_{\sigma\in{\mathcal{E}}_{K}}|D_{K,\sigma}|S_{K,\sigma}\cdot\phi$
$\displaystyle=\frac{\sqrt{d}}{2}\sum_{\sigma\in{\mathcal{E}}_{K}}|\sigma|[q_{\sigma}-q_{K}-\overline{\nabla}_{K}q\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K})]\mathbf{n}_{K,\sigma}\cdot\phi$
$\displaystyle=\frac{\sqrt{d}}{2}\sum_{\sigma\in{\mathcal{E}}_{K}}|\sigma|(q_{\sigma}-q_{K})\mathbf{n}_{K,\sigma}\cdot\phi$
$\displaystyle\,\,-\frac{\sqrt{d}}{2}\sum_{\sigma\in{\mathcal{E}}_{K}}|\sigma|[\overline{\nabla}_{K}q\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K})]\mathbf{n}_{K,\sigma}\cdot\phi.$

Using (3), we obtain

$\sum_{\sigma\in{\mathcal{E}}_{K}}|\sigma|(q_{\sigma}-q_{K})\mathbf{n}_{K,\sigma}\cdot\phi=|K|\overline{\nabla}_{K}q\cdot\phi.$
(9)

We then note that $\overline{\nabla}_{K}q$ is a constant vector, and hence,
$\overline{\nabla}_{K}q\cdot({\boldsymbol{x}}-{\boldsymbol{x}}_{K})$ is a
polynomial of degree one. This means that if ${\boldsymbol{x}}_{\sigma}$ is
the edge midpoint of $\sigma$, then

$|\sigma|[\overline{\nabla}_{K}q\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K})]\mathbf{n}_{K,\sigma}\cdot\phi=\int_{\sigma}\overline{\nabla}_{K}q\cdot({\boldsymbol{x}}-{\boldsymbol{x}}_{K})\mathbf{n}_{K,\sigma}\cdot\phi\,\,ds.$

Taking the sum over $\sigma\in{\mathcal{E}}_{K}$, using Green’s Theorem and
the fact that $\phi$ is constant, we then get

$\displaystyle\sum_{\sigma\in{\mathcal{E}}_{K}}|\sigma|[\overline{\nabla}_{K}q\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K})]\mathbf{n}_{K,\sigma}\cdot\phi$
$\displaystyle=\sum_{\sigma\in{\mathcal{E}}_{K}}\int_{\sigma}\overline{\nabla}_{K}q\cdot({\boldsymbol{x}}-{\boldsymbol{x}}_{K})\mathbf{n}_{K,\sigma}\cdot\phi\,\,ds$
$\displaystyle=\int_{K}{\mathop{\rm
div}}(\phi\overline{\nabla}_{K}q\cdot({\boldsymbol{x}}-{\boldsymbol{x}}_{K}))\,\,d{\boldsymbol{x}}$
$\displaystyle=\int_{K}\overline{\nabla}_{K}q\cdot\phi\,\,d{\boldsymbol{x}},$

or equivalently,

$\sum_{\sigma\in{\mathcal{E}}_{K}}|\sigma|[\overline{\nabla}_{K}q\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K})]\mathbf{n}_{K,\sigma}\cdot\phi=|K|\overline{\nabla}_{K}q\cdot\phi.$
(10)

The proof is then concluded by combining the expressions (9) and (10).

##  4\. Advective fluxes

###  4.1. Cell-centered first and second-order upwind fluxes

In this section, we discuss the advective fluxes, starting with the standard
first order upwind scheme for Cartesian meshes. For each control volume $K$,
we assign one discrete unknown $c_{K}$, which approximates the average value
of $c$ at the center ${\boldsymbol{x}}_{K}$ of $K$. Consider now a cell $K$
with eastern edge $\sigma$. Adapting the compass notation, the cells to the
west and east of $K$ are denoted by $W,E$ respectively. The cell to the east
of $E$ is then denoted by $EE$ (see Figure 4).

Figure 4. Cells involved in computing a second-order upwind flux for edge
$\sigma$.

Denoting by ${\boldsymbol{x}}_{\sigma}$ the midpoint of the edge $\sigma$ and
using the midpoint rule for computing the integral
$\int_{\sigma}c\mathbf{V}\cdot\mathbf{n}_{K,\sigma}\,\,ds$, the upwind fluxes
are then given by

$F_{K,\sigma}^{A}=c_{K}\mathbf{V}_{K,\sigma}^{+}-c_{E}\mathbf{V}_{K,\sigma}^{-},$
(11)

where
$\mathbf{V}_{K,\sigma}=\frac{1}{|\sigma|}\int_{\sigma}\mathbf{V}\cdot\mathbf{n}_{K,\sigma}\,\,ds$,$\mathbf{V}_{K,\sigma}^{+}=\max(\mathbf{V}_{K,\sigma},0)$,$\mathbf{V}_{K,\sigma}^{-}=\max(-\mathbf{V}_{K,\sigma},0)$,
and the value $c_{\sigma}=c({\boldsymbol{x}}_{\sigma})$ is approximated from
the upwind direction.

Essentially, this tells us that if material is flowing out (into) cell $K$
through the edge $\sigma$, corresponding to
$\mathbf{V}\cdot\mathbf{n}_{K,\sigma}$ being positive (negative), then the
value of $c_{\sigma}$ is approximated by $c_{K}$ ($c_{E}$), which makes sense
since the material traveling through $\sigma$ comes from cell $K$ ($E$). From
this, we see that the first-order upwind scheme has a very natural physical
interpretation, and is fairly simple to implement. Moreover, the formulation
(11) can straightforwardly be extended onto generic meshes, by treating $E$ as
a generic neighboring cell that shares $\sigma$ with $K$. However, the upwind
scheme may introduce too much numerical diffusion. Several methods can be
employed in order to mitigate this, but for this work, we focus on second-
order methods.

We now discuss the second-order upwind scheme. Again, we look at the eastern
edge $\sigma$ of cell $K$. The second-order upwind scheme approximates the
advective flux $F_{K,\sigma}$ in the following manner:

$F_{K,\sigma}^{A}=\bigg{(}c_{K}+\frac{1}{2}(c_{K}-c_{W})\bigg{)}\mathbf{V}_{K,\sigma}^{+}-\bigg{(}c_{E}+\frac{1}{2}(c_{E}-c_{EE})\bigg{)}\mathbf{V}_{K,\sigma}^{-}.$
(12)

Here, we see that the main difference between the second-order upwind flux
(12) and the first-order upwind flux is that (11) approximates the value of
$c_{\sigma}$ by a first-order Taylor expansion centered at
${\boldsymbol{x}}_{K}$ (${\boldsymbol{x}}_{E}$) when
$\mathbf{V}\cdot\mathbf{n}_{K,\sigma}$ is positive (negative), whereas (12)
gives a second-order Taylor expansion in the same direction. Another thing we
notice is that due to the presence of $-c_{W}$ and $-c_{E}$ in the second-
order upwind flux (12), it is no longer guaranteed that the discrete maximum
principle is satisfied. This is typically the case for linear schemes which
are second-order or higher. We now note that

$\frac{1}{2}(c_{K}-c_{W})\approx\frac{\partial c}{\partial
x}({\boldsymbol{x}}_{K})|{\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K}|=\nabla
c\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K}),$

and similarly

$\frac{1}{2}(c_{EE}-c_{E})\approx\frac{\partial c}{\partial
x}({\boldsymbol{x}}_{E})|{\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{E}|=\nabla
c\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{E}).$

Here, the gradient of $c$ is approximated from the upwind direction. To be
specific, we see that if $\mathbf{V}_{K,\sigma}>0$, then $\frac{\partial
c}{\partial x}({\boldsymbol{x}}_{K})$ is approximated by values from cells $K$
and $W$, which are located along the upwind direction.

To generalise, second-order fluxes are obtained from first-order upwind fluxes
by adding a correction term, where the correction term is related to a
discrete gradient. In particular, the fluxes $F_{K,\sigma}^{A}$ for second-
order schemes can be written in the following manner:

$F_{K,\sigma}^{A}=\bigg{(}c_{K}+\widetilde{\nabla}_{{\mathcal{D}}}c_{K}\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K})\bigg{)}\mathbf{V}_{K,\sigma}^{+}-\bigg{(}c_{E}+\widetilde{\nabla}_{{\mathcal{D}}}c_{E}\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{E})\bigg{)}\mathbf{V}_{K,\sigma}^{-},$
(13)

where $\widetilde{\nabla}_{{\mathcal{D}}}c_{K}$ is a discrete gradient which
approximates the value of $\nabla c$ at cell $K$. Different choices on how to
reconstruct the discrete gradient $\widetilde{\nabla}_{{\mathcal{D}}}c$ then
leads to different schemes. For example, we can take
$\widetilde{\nabla}_{{\mathcal{D}}}c_{K}$ to be $\overline{\nabla}_{K}c$
defined as in (3). For schemes with cell-centered unknowns, the value of
$c_{\sigma}$ in (3) needs to be chosen. Taking the value of $c_{\sigma}$ from
the upwind direction then leads to the second-order upwind scheme. On the
other hand, taking the value of $c_{\sigma}$ from the downwind direction leads
to the second-order centered scheme. To extend the description of the fluxes
(13) onto non-Cartesian meshes, we treat $E$ as a generic neighboring cell
that shares the edge $\sigma$ with $K$. As with Cartesian meshes, using a
second-order cell-centered scheme leads to a much wider stencil. The main
point to consider is how the term $\widetilde{\nabla}_{{\mathcal{D}}}c_{K}$ in
(13) is computed on generic meshes. Here, we use $\overline{\nabla}_{K}c$
defined as in (3) as our discrete gradient. This gives us

$\widetilde{\nabla}_{{\mathcal{D}}}c_{K}=\frac{1}{|K|}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|(c_{\sigma}^{\mathrm{up}}-c_{K})\mathbf{n}_{K,\sigma},$
(14)

where the value of $c_{\sigma}^{\mathrm{up}}$ is then taken from the upwind
direction. That is, we take $c_{\sigma}^{\mathrm{up}}=c_{K}$ if
$\mathbf{V}_{K,\sigma}>0$ and $c_{\sigma}^{\mathrm{up}}$ from the neighboring
cell otherwise.

As an example, consider a cell $K$ of a triangular mesh with edges
$\sigma_{1},\sigma_{2},\sigma_{3}$ being shared with cells $E,W,S$
respectively. Taking a velocity field $\mathbf{V}$ such that
$\mathbf{V}_{K,\sigma_{1}}>0$ and $\mathbf{V}_{K,\sigma_{i}}<0$ for the other
edges $\sigma_{i},i=2,3$ of $K$ (see Figure 5), the discrete derivative is
then computed to be

$\widetilde{\nabla}_{{\mathcal{D}}}c_{K}=\frac{1}{|K|}\big{(}|\sigma_{2}|(c_{W}-c_{K})\mathbf{n}_{K,\sigma_{2}}+|\sigma_{3}|(c_{S}-c_{K})\mathbf{n}_{K,\sigma_{3}}\big{)}.$

Figure 5. Cell values needed for constructing $\overline{\nabla}_{K}c$ for
cell-centered schemes.

Since the construction of the advective fluxes $F_{K,\sigma}^{A}$ on an edge
$\sigma$ of cell $K$ in (13) requires a linear reconstruction both in cell $K$
and from its neighboring cell, this requires information not only from the
neighboring cells of $K$, but also from the neighbors of its neighbors.
Moreover, this generalisation of the second-order upwind fluxes onto arbitrary
polygonal meshes comes with an additional computational cost/storage. In
particular, the sign of $\mathbf{V}\cdot\mathbf{n}_{K,\sigma}$ needs to be
computed/stored for each of the edges of cell $K$ and its neighbors in order
to determine the value of $c_{\sigma}^{\mathrm{up}}$ in computing the discrete
gradient (14).

###  4.2. Fully local second-order fluxes

In this section, we introduce the idea of having a fully local second-order
scheme. In order to do so, we employ the concept of hybridisation, as in [1].
This was done for first order upwind schemes in [3] by the introduction of
additional unknowns, one on each face, leading to the hybridised first-order
upwind fluxes

$F_{K,\sigma}^{A}=c_{K}\mathbf{V}_{K,\sigma}^{+}-c_{\sigma}\mathbf{V}_{K,\sigma}^{-}.$
(15)

By doing so, the direct dependence of the fluxes on values from neighboring
cells is eliminated. The main idea for this section is to extend the
hybridised first-order upwind fluxes into second-order fluxes by writing

$F_{K,\sigma}^{A}=\bigg{(}c_{K}+\widetilde{\nabla}_{{\mathcal{D}}}c_{K}\cdot({\boldsymbol{x}}-{\boldsymbol{x}}_{K})\bigg{)}\mathbf{V}_{K,\sigma}^{+}-c_{\sigma}\mathbf{V}_{K,\sigma}^{-}.$
(16)

One of the main advantages of using a hybridised scheme is, if information is
expected to have arrived from another cell, we keep the value of $c$ at
$\sigma$ to be implicit. Moreover, for second-order schemes, if the discrete
gradient $\widetilde{\nabla}_{\mathcal{D}}c$ is taken such that
$(\widetilde{\nabla}_{\mathcal{D}}c)_{|K}=\overline{\nabla}_{K}c$ as defined
in (3), then the values of $c_{\sigma}$ in the discrete gradient are
determined naturally, without having to choose between the values in the
upwind and downwind direction. This also allows us to save in terms of the
computational cost associated with computing
$(\widetilde{\nabla}_{\mathcal{D}}c)$.

On the other hand, since we introduced additional unknowns along the edges, we
also need to introduce additional equations (one corresponding to each edge).
Denoting by $N_{K}$ and $N_{e}$ denote the number of cells and edges in the
mesh, respectively, we then have a scheme with $N_{K}+N_{e}$ equations and
$N_{K}+N_{e}$ unknowns. Using definitions (2) and (16) for the diffusive and
advective fluxes, and denoting

$F_{K,\sigma}=F_{K,\sigma}^{D}+F_{K,\sigma}^{A},$

the first $N_{k}$ equations are given by the balance of fluxes:

$\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|F_{K,\sigma}=|K|f_{K}.$ (17a)
Following this, we impose the conservation of fluxes for interior edges. That
is, for cells $K$ and $L$ that share a common edge $\sigma$,
$F_{K,\sigma}+F_{L,\sigma}=0.$ (17b) Finally, on each boundary edge, an
equation is needed for imposing the Dirichlet boundary conditions. In
particular, for each $\sigma\in{\mathcal{E}}_{\mathrm{ext}}$, we take
$c_{\sigma}$ to be the average value of $g$ over $\sigma$, given by
$c_{\sigma}=\frac{1}{|\sigma|}\int_{\sigma}g\,\,ds.$ (17c)

Note here that the scheme (17a)-(17c) consists of $N_{e}$ more equations than
cell-centered schemes. The main advantage, however, of the hybridised
formulation is that the fluxes no longer depend directly on values from
neighboring cells, which enables us to use these second-order fluxes even on
cells at the boundary of the domain. Moreover, we can employ static
condensation for solving the system of equations, which essentially leads to
solving only a system of $N_{e}$ equations in $N_{e}$ unknowns.

##  5\. Convergence of the scheme

In this section, we study the convergence of a family of hybridised second-
order finite volume schemes, with diffusive and advective fluxes as in (2) and
(16), respectively. Since the diffusive fluxes used here come from the HMM
method, the convergence analysis for these terms are more or less the same as
those found in [3, 6, 9]. For our proofs, we focus on the advective fluxes,
and how to deal with the correction term. For simplicity of exposition, we
assume that the discrete gradient for the correction term in (16) is taken
from (6). A typical norm used for measuring errors is a discrete $L^{2}$ norm.
For this, we start by defining a function reconstruction
$\Pi_{{\mathcal{D}}_{h}}:X_{{\mathcal{D}}_{h}}\rightarrow L^{2}(\Omega)$, such
that for all $q_{h}\in X_{{\mathcal{D}}_{h}}$

$(\Pi_{{\mathcal{D}}_{h}}q_{h})_{|K}=q_{K},\quad\forall
K\in{\mathcal{M}}_{h}.$ (18)

The discrete $L^{2}$ norm is then defined as

$\left\|q_{h}\right\|_{L^{2}(\Omega)}:=\left\|\Pi_{{\mathcal{D}}_{h}}q_{h}\right\|_{L^{2}(\Omega)}=\bigg{(}\sum_{K\in{\mathcal{M}}_{h}}|K|q_{K}^{2}\bigg{)}^{1/2},\quad\forall
q_{h}\in X_{{\mathcal{D}}_{h}}.$ (19)

Following the ideas in [3, 5], it is also useful to perform analysis for
hybridised schemes by using a discrete $H^{1}$-like norm on
$X_{{\mathcal{D}}_{h}}$:

$\left\|q_{h}\right\|_{1,{\mathcal{D}}_{h}}=\bigg{(}\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{\mathcal{E}}_{K}}\frac{|\sigma|}{d_{K,\sigma}}|q_{K}-q_{\sigma}|^{2}\bigg{)}^{1/2},\quad\forall
q_{h}\in X_{{\mathcal{D}}_{h}}.$ (20)

We start by specifying the following estimate on the discrete gradient
$\nabla_{{\mathcal{D}}}$.

######  Lemma 5.1.

Let $q_{h}\in X_{{\mathcal{D}}_{h}}$, then

$\left\|\nabla_{{\mathcal{D}}_{h}}q_{h}\right\|_{L^{2}(\Omega)}\lesssim\left\|q_{h}\right\|_{1,{\mathcal{D}}_{h}},$
(21)

where $a\lesssim b$ means that there is a constant $C$, independent of $h$,
such that $a\leq Cb$.

Proof.  Using the definition (6) of the discrete gradient and the
orthogonality of the stabilisation term (7), we have that

$\displaystyle\int_{\Omega}\nabla_{{\mathcal{D}}_{h}}q_{h}\cdot\nabla_{{\mathcal{D}}_{h}}q_{h}\,\,d{\boldsymbol{x}}$
$\displaystyle=\sum_{K\in{\mathcal{M}}_{h}}\int_{K}\nabla_{{\mathcal{D}}_{h}}q_{h}\cdot\nabla_{{\mathcal{D}}_{h}}q_{h}\,\,d{\boldsymbol{x}}$
$\displaystyle=\sum_{K\in{\mathcal{M}}_{h}}\int_{K}\overline{\nabla}_{K}q_{h}\cdot\overline{\nabla}_{K}q_{h}\,\,d{\boldsymbol{x}}+\sum_{K\in{\mathcal{M}}_{h}}\int_{K}S_{k}(q_{h})\cdot
S_{K}(q_{h})\,\,d{\boldsymbol{x}}.$

We start by looking at the consistent term and use the definition (3) to
obtain

$\displaystyle\int_{\Omega}\overline{\nabla}_{{\mathcal{D}}_{h}}q_{h}\cdot\overline{\nabla}_{{\mathcal{D}}_{h}}q_{h}\,\,d{\boldsymbol{x}}$
$\displaystyle=\sum_{K\in{\mathcal{M}}_{h}}\int_{K}\overline{\nabla}_{K}q_{h}\cdot\overline{\nabla}_{K}q_{h}\,\,d{\boldsymbol{x}}$
$\displaystyle=\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|(q_{\sigma}-q_{K})\overline{\nabla}_{K}q_{h}\cdot\mathbf{n}_{K,\sigma}.$

Applying Cauchy-Schwarz, we then have

$\displaystyle\left\|\overline{\nabla}_{{\mathcal{D}}_{h}}q_{h}\right\|_{L^{2}(\Omega)}^{2}$
$\displaystyle\leq\sum_{K\in{\mathcal{M}}_{h}}\bigg{(}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}\frac{|\sigma|}{d_{K,\sigma}}(q_{\sigma}-q_{K})^{2}\bigg{)}^{1/2}\bigg{(}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|d_{K,\sigma}(\overline{\nabla}_{K}q_{h}\cdot\mathbf{n}_{K,\sigma})^{2}\bigg{)}^{1/2}$
$\displaystyle\leq\bigg{(}\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}\frac{|\sigma|}{d_{K,\sigma}}(q_{\sigma}-q_{K})^{2}\bigg{)}^{1/2}\bigg{(}\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|d_{K,\sigma}(\overline{\nabla}_{K}q_{h}\cdot\mathbf{n}_{K,\sigma})^{2}\bigg{)}^{1/2}$
$\displaystyle\leq\left\|q_{h}\right\|_{1,{\mathcal{D}}_{h}}\bigg{(}\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|d_{K,\sigma}(\overline{\nabla}_{K}q_{h}\cdot\mathbf{n}_{K,\sigma})^{2}\bigg{)}^{1/2}.$
(22)

Using the fact that

$(\overline{\nabla}_{K}q_{h}\cdot\mathbf{n}_{K,\sigma})^{2}\leq\overline{\nabla}_{K}q_{h}\cdot\overline{\nabla}_{K}q_{h}$

and

$\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|d_{K,\sigma}=d|K|,$

we then have

$\displaystyle\bigg{(}\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|d_{K,\sigma}(\overline{\nabla}_{K}q_{h}\cdot\mathbf{n}_{K,\sigma})^{2}\bigg{)}^{1/2}$
$\displaystyle\leq\sqrt{d}\bigg{(}\sum_{K\in{\mathcal{M}}_{h}}|K|\overline{\nabla}_{K}q_{h}\cdot\overline{\nabla}_{K}q_{h}\bigg{)}^{1/2}$
$\displaystyle=\sqrt{d}\left\|\overline{\nabla}_{{\mathcal{D}}_{h}}q_{h}\right\|_{L^{2}(\Omega)},$

which upon substitution to (22), leads to

$\left\|\overline{\nabla}_{{\mathcal{D}}_{h}}q_{h}\right\|_{L^{2}(\Omega)}^{2}\leq
d\left\|q_{h}\right\|_{1,{\mathcal{D}}_{h}}^{2}.$ (23)

Now, for the stabilisation term, we have from (4), (5), and (8) that

$\displaystyle\sum_{K\in{\mathcal{M}}_{h}}\int_{K}S_{k}(q_{h})\cdot
S_{K}(q_{h})\,\,d{\boldsymbol{x}}$
$\displaystyle=\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{\mathcal{E}}_{K}}|D_{K,\sigma}|S_{K,\sigma}\cdot
S_{K,\sigma}$
$\displaystyle=d\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{\mathcal{E}}_{K}}\frac{|D_{K,\sigma}|}{d_{K,\sigma}^{2}}(q_{\sigma}-q_{K}-\overline{\nabla}_{K}q\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K}))^{2}$
$\displaystyle=\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{\mathcal{E}}_{K}}\frac{|\sigma|}{d_{K,\sigma}}(q_{\sigma}-q_{K}-\overline{\nabla}_{K}q\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K}))^{2}.$

Using the fact that $(a+b)^{2}\leq 2a^{2}+2b^{2}$ for any real numbers $a,b$,
we then have

$\displaystyle\sum_{K\in{\mathcal{M}}_{h}}\int_{K}S_{k}(q_{h})\cdot
S_{K}(q_{h})\,\,d{\boldsymbol{x}}$ $\displaystyle\leq
2\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{\mathcal{E}}_{K}}\frac{|\sigma|}{d_{K,\sigma}}(q_{\sigma}-q_{K})^{2}$
$\displaystyle\,\,+2\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{\mathcal{E}}_{K}}\frac{|\sigma|}{d_{K,\sigma}}(\overline{\nabla}_{K}q\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K}))^{2}.$

Now, the first term is simply
$\left\|q_{h}\right\|_{1,{\mathcal{D}}_{h}}^{2}$, so we only need to deal with
the second term, which we will denote by $T_{2}$. Using Cauchy-Schwarz, the
regularity of the mesh (MR2), and (23) we have that

$\displaystyle T_{2}$ $\displaystyle\leq
2\sum_{K\in{\mathcal{M}}_{h}}(\overline{\nabla}_{K}q\cdot\overline{\nabla}_{K}q)\sum_{\sigma\in{\mathcal{E}}_{K}}\frac{|\sigma|}{d_{K,\sigma}}\big{(}({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K})\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K})\big{)}$
$\displaystyle\lesssim\sum_{K\in{\mathcal{M}}_{h}}|K|(\overline{\nabla}_{K}q\cdot\overline{\nabla}_{K}q)$
$\displaystyle\lesssim\left\|q_{h}\right\|_{1,{\mathcal{D}}_{h}}^{2}.$

This allows us to conclude that

$\sum_{K\in{\mathcal{M}}_{h}}\int_{K}S_{k}(q_{h})\cdot
S_{K}(q_{h})\,\,d{\boldsymbol{x}}\lesssim\left\|q_{h}\right\|_{1,{\mathcal{D}}_{h}}^{2}.$
(24)

Combining the inequalities (23) and (24) then concludes the proof.

###  5.1. A priori estimates

######  Lemma 5.2.

Let us assume that (A1)–(A3) hold. Let ${\mathcal{M}}_{h}$ be an admissible
discretization of $\Omega$ such that
$\theta>\mathrm{regul}({\mathcal{M}}_{h})$, and let $F_{K,\sigma}^{A}$ be the
advective flux of $q_{h}\in X_{{\mathcal{D}}_{h}}$ given by (16) for the
velocity field $\mathbf{V}\in C^{1}(\overline{\Omega})^{d}$. Then there exists
a non-negative constant $C_{1}\geq 0$ that only depends on $\theta$ and
$\mathbf{V}$ such that

$\displaystyle\forall q_{h}\in X_{{\mathcal{D}}_{h}}:$ (25)
$\displaystyle\,\,\,\frac{1}{2}\int_{\Omega}(\Pi_{{\mathcal{D}}_{h}}q_{h})^{2}{\mathop{\rm
div}}(\mathbf{V})\,\,d{\boldsymbol{x}}\leq\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{\mathcal{E}}_{K}}|\sigma|F_{K,\sigma}^{A}(q_{K}-q_{\sigma})+C_{1}h\left\|q_{h}\right\|_{1,{\mathcal{D}}_{h}}^{2}.$

Proof.

We start by taking note that

$\displaystyle\mathbf{V}_{K,\sigma}^{+}$
$\displaystyle=\frac{1}{2}\mathbf{V}_{K,\sigma}+\frac{1}{2}(\mathbf{V}_{K,\sigma}^{+}+\mathbf{V}_{K,\sigma}^{-}),$
$\displaystyle\mathbf{V}_{K,\sigma}^{-}$
$\displaystyle=-\frac{1}{2}\mathbf{V}_{K,\sigma}+\frac{1}{2}(\mathbf{V}_{K,\sigma}^{+}+\mathbf{V}_{K,\sigma}^{-}),$

so that we can write the fluxes (16) as

$F_{K,\sigma}^{A}=\frac{1}{2}(q_{K}+q_{\sigma})\mathbf{V}_{K,\sigma}+\frac{1}{2}(\mathbf{V}_{K,\sigma}^{+}+\mathbf{V}_{K,\sigma}^{-})(q_{K}-q_{\sigma})+\nabla_{\mathcal{D}}q_{K}\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K})\mathbf{V}_{K,\sigma}^{+}.$
(26)

Using (26) and the fact that
$\mathbf{V}_{K,\sigma}+\mathbf{V}_{K^{\prime},\sigma}=0$ for two cells $K$ and
$K^{\prime}$ sharing $\sigma$, we then obtain

$\displaystyle\\!\\!\\!\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{\mathcal{E}}_{K}}|\sigma|F_{K,\sigma}^{A}(q_{K}-q_{\sigma})$
$\displaystyle=\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}\frac{|\sigma|}{2}(q_{K}^{2}-q_{\sigma}^{2})\mathbf{V}_{K,\sigma}+\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}\frac{|\sigma|}{2}(\mathbf{V}_{K,\sigma}^{+}+\mathbf{V}_{K,\sigma}^{-})(q_{K}-q_{\sigma})^{2}$
$\displaystyle\,\,\,+\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|\nabla_{\mathcal{D}}q_{K}\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K})(q_{K}-q_{\sigma})\mathbf{V}_{K,\sigma}^{+}$
$\displaystyle=\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}\frac{|\sigma|}{2}q_{K}^{2}\mathbf{V}_{K,\sigma}+\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}\frac{|\sigma|}{2}(\mathbf{V}_{K,\sigma}^{+}+\mathbf{V}_{K,\sigma}^{-})(q_{K}-q_{\sigma})^{2}$
$\displaystyle\,\,\,+\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|\nabla_{\mathcal{D}}q_{K}\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K})(q_{K}-q_{\sigma})\mathbf{V}_{K,\sigma}^{+}$
$\displaystyle\geq\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}\frac{|\sigma|}{2}q_{K}^{2}\mathbf{V}_{K,\sigma}+\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|\nabla_{\mathcal{D}}q_{K}\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K})(q_{K}-q_{\sigma})\mathbf{V}_{K,\sigma}^{+}.$

Using the fact that

$\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|\mathbf{V}_{K,\sigma}=\int_{K}{\mathop{\rm
div}}(\mathbf{V})\,\,d{\boldsymbol{x}},$

we then have

$\displaystyle\frac{1}{2}$
$\displaystyle\int_{\Omega}(\Pi_{{\mathcal{D}}_{h}}q_{h})^{2}{\mathop{\rm
div}}(\mathbf{V})\,\,d{\boldsymbol{x}}$ (27)
$\displaystyle\\!\\!\\!\leq\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{\mathcal{E}}_{K}}|\sigma|F_{K,\sigma}^{A}(q_{K}-q_{\sigma})+\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|\nabla_{\mathcal{D}}q_{K}\cdot({\boldsymbol{x}}_{K}-{\boldsymbol{x}}_{\sigma})(q_{K}-q_{\sigma})\mathbf{V}_{K,\sigma}^{+}.$

Now, upon applying Cauchy-Schwarz, we obtain

$\displaystyle\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|\nabla_{\mathcal{D}}q_{K}\cdot({\boldsymbol{x}}_{K}-{\boldsymbol{x}}_{\sigma})(q_{K}-q_{\sigma})\mathbf{V}_{K,\sigma}^{+}$
$\displaystyle\,\,\leq\sum_{K\in{\mathcal{M}}_{h}}\bigg{(}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|\bigg{(}\nabla_{\mathcal{D}}q_{K}\cdot({\boldsymbol{x}}_{K}-{\boldsymbol{x}}_{\sigma})\mathbf{V}_{K,\sigma}^{+}\bigg{)}^{2}\bigg{)}^{1/2}\bigg{(}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|(q_{K}-q_{\sigma})^{2}\bigg{)}^{1/2}$
$\displaystyle\,\,\leq\bigg{(}\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}\bigg{(}|\sigma|\nabla_{\mathcal{D}}q_{K}\cdot({\boldsymbol{x}}_{K}-{\boldsymbol{x}}_{\sigma})\mathbf{V}_{K,\sigma}^{+}\bigg{)}^{2}\bigg{)}^{1/2}\bigg{(}\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|(q_{K}-q_{\sigma})^{2}\bigg{)}^{1/2}$
$\displaystyle\,\,\lesssim\bigg{(}\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|(\nabla_{\mathcal{D}}q_{K}\cdot\nabla_{\mathcal{D}}q_{K})\big{(}({\boldsymbol{x}}_{K}-{\boldsymbol{x}}_{\sigma})\cdot({\boldsymbol{x}}_{K}-{\boldsymbol{x}}_{\sigma})\big{)}\bigg{)}^{1/2}h^{1/2}\left\|(q_{h},q_{{\mathcal{E}}_{h}})\right\|_{1,D_{h},\mathcal{E}_{h}}$
$\displaystyle\,\,\lesssim
h\left\|\nabla_{\mathcal{D}}q_{h}\right\|_{L^{2}(\Omega)}\left\|q_{h}\right\|_{1,{\mathcal{D}}_{h}},$

which, together with (21), leads to

$\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|\nabla_{\mathcal{D}}q_{K}\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K})(q_{K}-q_{\sigma})\mathbf{V}_{K,\sigma}^{+}\lesssim
h\left\|q_{h}\right\|_{1,{\mathcal{D}}_{h}}^{2}.$

Substituting the above inequality into (27) then leads us to (25).

######  Lemma 5.3.

Let us assume that (A1)-(A3) hold. Let ${\mathcal{M}}_{h}$ be an admissible
discretization of $\Omega$ such that $\theta>\mathrm{regul}(D_{h})$, and let
$F_{K,\sigma}^{A}$ be the advective flux given by (16) for $\mathbf{V}\in
C^{1}(\Omega)^{d}$. Then, for all solutions $c_{h}$ to the scheme (17a)-(17c)
with diffusive and advective fluxes defined as in (2) and (16), we have

$\left\|c_{h}\right\|_{1,{\mathcal{D}}_{h}}^{2}\lesssim\left\|f\right\|_{L^{2}(\Omega)}\left\|c_{h}\right\|_{L^{2}(\Omega)}+h\left\|c_{h}\right\|_{1,{\mathcal{D}}_{h}}^{2}.$
(28)

Proof.  Take

$\int_{\Omega}f\Pi_{{\mathcal{D}}_{h}}c_{h}\,\,d{\boldsymbol{x}}=\sum_{K\in{\mathcal{M}}_{h}}\int_{K}fc_{K}\,\,d{\boldsymbol{x}}$

and let $c_{h}$ be the solution to the scheme. Then, using the conservation of
fluxes, we have

$\displaystyle\int_{\Omega}f\Pi_{{\mathcal{D}}_{h}}c_{h}\,\,d{\boldsymbol{x}}$
$\displaystyle=\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|(F_{K,\sigma}^{D}+F_{K,\sigma}^{A})c_{K}$
$\displaystyle=\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|F_{K,\sigma}^{D}(c_{K}-c_{\sigma})+\sum_{K\in{\mathcal{M}}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|F_{K,\sigma}^{A}(c_{K}-c_{\sigma}).$

Applying Cauchy-Schwarz then leads us to

$\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}F_{K,\sigma}^{D}(c_{K}-c_{\sigma})+\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}F_{K,\sigma}^{A}(c_{K}-c_{\sigma})\leq\left\|f\right\|_{L^{2}(\Omega)}\left\|c_{h}\right\|_{L^{2}(\Omega)}.$

Using the definition (2) for the diffusive fluxes, we have

$\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|F_{K,\sigma}^{D}(c_{K}-c_{\sigma})=\sum_{K\in{\mathcal{M}}_{h}}\int_{K}\Lambda_{K}\nabla_{{\mathcal{D}}_{h}}c_{h}\cdot\nabla_{{\mathcal{D}}_{h}}c_{h}\,\,d{\boldsymbol{x}}.$

We then use the assumption (A1) on $\Lambda$ and [5, Lemma 13.11] to establish
that

$\left\|c_{h}\right\|_{1,{\mathcal{D}}_{h}}^{2}\lesssim\sum_{K\in{\mathcal{M}}_{h}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|F_{K,\sigma}^{D}(c_{K}-c_{\sigma}).$

Combining this with the inequality (25), and using the fact that ${\mathop{\rm
div}}(\mathbf{V})\geq 0$ then allows us to conclude the proof.

###  5.2. Convergence result

In this section, we show that the numerical solution of the scheme given by
(17a)-(17c), with diffusive and advective fluxes defined as in (2),(16)
converges to the weak solution of (1). For simplicity of exposition, we
consider homogeneous Dirichlet boundary conditions.

######  Definition 5.4 (Weak solution of the advection-diffusion problem (1)).

We say that $c\in H^{1}(\Omega)$ is the weak solution of (1) if for any
$\psi\in H^{1}(\Omega)$

$\int_{\Omega}\Lambda\nabla
c\cdot\nabla\psi\,\,d{\boldsymbol{x}}-\int_{\Omega}c\mathbf{V}\cdot\nabla\psi\,\,d{\boldsymbol{x}}=\int_{\Omega}f\psi\,\,d{\boldsymbol{x}}.$
(29)

######  Theorem 5.5.

Let $c\in H^{1}(\Omega)$ be the weak solution to (1). Under assumptions
(A1)–(A3), let $c_{h}\in X_{{\mathcal{D}}_{h}}$ be the numerical solution to
the scheme (17a)-(17c), with diffusive and advective fluxes constructed as in
(2) and (16). Then, for $h\rightarrow 0$, the following hold:

  1. (1)

$\Pi_{{\mathcal{D}}_{h}}c_{h}\rightarrow c$ in $L^{r}(\Omega)$ for all
$r<\frac{2d}{d-2}$

  2. (2)

$\overline{\nabla}_{{\mathcal{D}}_{h}}c_{h}\rightarrow\nabla c$ in
$L^{2}(\Omega)^{d}$.

Proof.  Using Lemma 5.3, we have
$\left\|c_{h}\right\|_{1,{\mathcal{D}}_{h}}^{2}\lesssim\left\|f\right\|_{L^{2}(\Omega)}\left\|\Pi_{{\mathcal{D}}_{h}}c_{h}\right\|_{L^{2}(\Omega)}$
when $h$ is small. In view of Lemma 9.1, we obtain an upper bound on
$\left\|c_{h}\right\|_{1,{\mathcal{D}}_{h}}$. Then the result of Lemma 9.2
implies the existence of a function $c\in H^{1}_{0}(\Omega)$ such that, up to
a subsequence, $\Pi_{{\mathcal{D}}_{h}}c_{h}\rightarrow c$ in $L^{r}(\Omega)$
for all $r<\frac{2d}{d-2}$ and
$\overline{\nabla}_{{\mathcal{D}}_{h}}c_{h}\rightarrow\nabla c$ weakly in
$L^{2}(\Omega)^{d}$. Now, consider $\varphi\in C_{c}^{\infty}(\Omega)$ and for
$K\in{\mathcal{M}}_{h},\sigma\in{\mathcal{E}}_{h}$ write
$\varphi_{K}=\varphi({\boldsymbol{x}}_{K}),\varphi_{\sigma}=\varphi({\boldsymbol{x}}_{\sigma})$.
Let $c_{h}\in X_{{\mathcal{D}}_{h}}$ be the numerical solution to (17a)-(17c).
We multiply the balance of flux equations (17a) by $\varphi_{K}$ and take the
sum over $K\in{\mathcal{M}}_{h}$ to obtain

$\displaystyle\sum_{K\in{\mathcal{M}}}\int_{K}f\varphi_{K}\,\,d{\boldsymbol{x}}$
$\displaystyle=\sum_{K\in{\mathcal{M}}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|F_{K,\sigma}\varphi_{K}.$

Using the conservation of internal fluxes, we then have

$\displaystyle\int_{\Omega}f\varphi_{h}\,\,d{\boldsymbol{x}}$
$\displaystyle=\sum_{K\in{\mathcal{M}}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|F_{K,\sigma}(\varphi_{K}-\varphi_{\sigma})$
$\displaystyle=\sum_{K\in{\mathcal{M}}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|F_{K,\sigma}^{D}(\varphi_{K}-\varphi_{\sigma})$
$\displaystyle\,\,+\sum_{K\in{\mathcal{M}}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|\bigg{(}\big{(}c_{K}+\nabla_{{\mathcal{D}}_{h}}c_{K}\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K})\big{)}\mathbf{V}_{K,\sigma}^{+}-c_{\sigma}\mathbf{V}_{K,\sigma}^{-}\bigg{)}(\varphi_{K}-\varphi_{\sigma})$
$\displaystyle=\sum_{K\in{\mathcal{M}}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|F_{K,\sigma}^{D}(\varphi_{K}-\varphi_{\sigma})+\sum_{K\in{\mathcal{M}}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|c_{K}\mathbf{V}_{K,\sigma}(\varphi_{K}-\varphi_{\sigma})$
$\displaystyle\,\,+\sum_{K\in{\mathcal{M}}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|(c_{K}-c_{\sigma})\mathbf{V}_{K,\sigma}^{-}(\varphi_{K}-\varphi_{\sigma})$
$\displaystyle\,\,+\sum_{K\in{\mathcal{M}}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|\big{(}\nabla_{{\mathcal{D}}_{h}}c_{K}\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K})\big{)}\mathbf{V}_{K,\sigma}^{+}(\varphi_{K}-\varphi_{\sigma}).$

We then write the right hand side of the above equation as
$T_{1}+T_{2}+T_{3}+T_{4}$. For the term $T_{1},$ we write a Taylor expansion

$\varphi_{\sigma}=\varphi_{K}+\nabla\varphi({\boldsymbol{x}}_{K})\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K})+R_{K,\sigma}(\varphi),$

where $R_{K,\sigma}(\varphi)\lesssim
h^{2}\left\|\nabla^{2}\varphi\right\|_{\infty}$. This leads to

$\displaystyle T_{1}$
$\displaystyle=\sum_{K\in{\mathcal{M}}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|F_{K,\sigma}^{D}\nabla\varphi({\boldsymbol{x}}_{K})\cdot({\boldsymbol{x}}_{K}-{\boldsymbol{x}}_{\sigma})+\sum_{K\in{\mathcal{M}}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|F_{K,\sigma}^{D}R_{K,\sigma}(\varphi)$
$\displaystyle=T_{1,1}+T_{1,2}.$

Using the definition (2) of the diffusive flux and the orthogonality of the
stabilisation term (7), we then get

$\displaystyle T_{1,1}$
$\displaystyle=\sum_{K\in{\mathcal{M}}_{h}}\int_{K}\nabla_{{\mathcal{D}}_{h}}c_{h}\cdot\nabla\varphi({\boldsymbol{x}}_{K})\,\,d{\boldsymbol{x}}$
$\displaystyle=\int_{\Omega}\overline{\nabla}_{{\mathcal{D}}_{h}}c_{h}\cdot(\nabla\varphi)_{h}\,\,d{\boldsymbol{x}},$

where $(\nabla\varphi)_{h}$ is defined such that
$\big{(}(\nabla\varphi)_{h}\big{)}_{|K}=\varphi({\boldsymbol{x}}_{K})$ for
each $K\in{\mathcal{M}}_{h}$. Using the weak convergence of
$\overline{\nabla}_{{\mathcal{D}}_{h}}c_{h}$ and the strong convergence of
$(\nabla\varphi)_{h}$, we have that as $h\rightarrow 0$,

$T_{1,1}\rightarrow\int_{\Omega}\Lambda\nabla
c\cdot\nabla\varphi\,\,d{\boldsymbol{x}}.$ (30)

Also, as $h\rightarrow 0$, $T_{1,2}\rightarrow 0$. Now, considering $T_{2}$,
we have that

$\displaystyle T_{2}$
$\displaystyle=\sum_{K\in{\mathcal{M}}}c_{K}\varphi_{K}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|\mathbf{V}_{K,\sigma}-\sum_{K\in{\mathcal{M}}}c_{K}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|\mathbf{V}_{K,\sigma}\varphi_{\sigma}$
$\displaystyle=\sum_{K\in{\mathcal{M}}}\int_{K}c_{K}\varphi_{K}{\mathop{\rm
div}}(\mathbf{V})\,\,d{\boldsymbol{x}}-\sum_{K\in{\mathcal{M}}}c_{K}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}\int_{\sigma}\varphi\mathbf{V}\cdot\mathbf{n}_{K,\sigma}\,\,ds$
$\displaystyle\,\,+\sum_{K\in{\mathcal{M}}}c_{K}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}\int_{\sigma}(\varphi-\varphi_{\sigma})\mathbf{V}\cdot\mathbf{n}_{K,\sigma}\,\,ds$
$\displaystyle=\int_{\Omega}\Pi_{{\mathcal{D}}_{h}}c_{h}\varphi_{h}{\mathop{\rm
div}}(\mathbf{V})\,\,d{\boldsymbol{x}}-\int_{\Omega}\Pi_{{\mathcal{D}}_{h}}c_{h}{\mathop{\rm
div}}(\varphi\mathbf{V})\,\,d{\boldsymbol{x}}$
$\displaystyle\,\,+\sum_{K\in{\mathcal{M}}}c_{K}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}\int_{\sigma}(\varphi-\varphi_{\sigma})\mathbf{V}\cdot\mathbf{n}_{K,\sigma}\,\,ds.$

Due to the strong convergence of $\Pi_{{\mathcal{D}}_{h}}c_{h}$, we see that
the first two terms on the right hand side converges to
$\int_{\Omega}c\varphi{\mathop{\rm
div}}(\mathbf{V})\,\,d{\boldsymbol{x}}-\int_{\Omega}c{\mathop{\rm
div}}(\varphi V)\,\,d{\boldsymbol{x}}$. We now consider the third term. Due to
the fact that $\mathbf{V}_{K,\sigma}$ is conservative, we may write

$\displaystyle\bigg{|}\sum_{K\in{\mathcal{M}}}c_{K}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}\int_{\sigma}(\varphi-\varphi_{\sigma})\mathbf{V}\cdot\mathbf{n}_{K,\sigma}\,\,ds\bigg{|}$
$\displaystyle=\bigg{|}\sum_{K\in{\mathcal{M}}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}(c_{K}-c_{\sigma})\int_{\sigma}(\varphi-\varphi_{\sigma})\mathbf{V}\cdot\mathbf{n}_{K,\sigma}\,\,ds\bigg{|}$
(31) $\displaystyle\lesssim
h\left\|\nabla\varphi\right\|_{\infty}\sum_{K\in{\mathcal{M}}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma||c_{K}-c_{\sigma}|.$

Using Cauchy-Schwarz and the boundedness of
$\left\|c_{h}\right\|_{1,{\mathcal{D}}_{h}}$, we have

$\sum_{K\in{\mathcal{M}}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma||c_{K}-c_{\sigma}|\leq(d|\Omega|)^{1/2}\left\|c_{h}\right\|_{1,{\mathcal{D}}_{h}}\lesssim
1.$

Thus, as $h\rightarrow 0$, we have

$T_{2}\rightarrow\int_{\Omega}c\varphi{\mathop{\rm
div}}(\mathbf{V})\,\,d{\boldsymbol{x}}-\int_{\Omega}c{\mathop{\rm
div}}(\varphi\mathbf{V})\,\,d{\boldsymbol{x}}.$ (32)

For $T_{3}$, we use the boundedness $|\mathbf{V}_{K,\sigma}|\lesssim 1$ and an
argument similar to (31) in order to establish that as $h\rightarrow 0$,

$T_{3}\rightarrow 0.$ (33)

Finally, we consider the term $T_{4}$. Here, we have

$\displaystyle
T_{4}\lesssim\left\|\nabla\varphi\right\|_{\infty}\sum_{K\in{\mathcal{M}}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|\big{|}\nabla_{{\mathcal{D}}_{h}}c_{K}\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K})\big{|},$

By Cauchy-Schwarz, we then have

$\displaystyle\sum_{K\in{\mathcal{M}}}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|\big{|}\nabla_{{\mathcal{D}}_{h}}c_{K}\cdot({\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K})\big{|}$
$\displaystyle\,\,\,\leq\sum_{K\in{\mathcal{M}}}(\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma||{\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K}|\nabla_{{\mathcal{D}}_{h}}c_{K}\cdot\nabla_{{\mathcal{D}}_{h}}c_{K})^{1/2}(\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma||{\boldsymbol{x}}_{\sigma}-{\boldsymbol{x}}_{K}|)^{1/2}$
$\displaystyle\,\,\,\lesssim
h\sum_{K\in{\mathcal{M}}}\int_{K}\nabla_{{\mathcal{D}}_{h}}c_{K}\cdot\nabla_{{\mathcal{D}}_{h}}c_{K}\,\,d{\boldsymbol{x}},$

which, together with (21), leads to

$\displaystyle T_{4}\lesssim
h\left\|\nabla\varphi\right\|_{\infty}\left\|c_{h}\right\|_{1,{\mathcal{D}}_{h}}.$

Hence, as $h\rightarrow 0$,

$T_{4}\rightarrow 0.$ (34)

Combining the results in (30), (32), (33), (34) then shows us that indeed the
numerical solution converges to the weak solution (29) of the advection-
diffusion problem.

Following Step 3 of the proof in [3, Theorem 3.7], we can establish that the
weak convergence $\overline{\nabla}_{{\mathcal{D}}_{h}}c_{h}\rightarrow\nabla
c$ is, in fact, strong.

We note here that key to the convergence analysis is the property (21) of the
gradient used for the linear term in (16). Hence, flux-limited second-order
upwind fluxes, such as those in [4], are also covered by the analysis
presented above. In particular, we note that these fluxes take the form

$F_{K,\sigma}^{A}=\bigg{(}c_{K}+\phi_{K}(c)\widetilde{\nabla}_{{\mathcal{D}}}c_{K}\cdot({\boldsymbol{x}}-{\boldsymbol{x}}_{K})\bigg{)}\mathbf{V}_{K,\sigma}^{+}-c_{\sigma}\mathbf{V}_{K,\sigma}^{-},$
(35)

where $\phi_{K}:X_{\mathcal{D}}\rightarrow[0,1]$. Since $\phi$ is bounded
between $0$ and $1$, $\phi\widetilde{\nabla}_{\mathcal{D}}q$ satisfies the
estimate (21) for any $q\in X_{\mathcal{D}}$; thus, the convergence results
still hold when the advective fluxes (16) are replaced with (35).

##  6\. Numerical tests

In this section, we present numerical tests of the fully local second-order
upwind scheme defined by (17a)-(17c) over the domain $\Omega$, with diffusive
and advective fluxes defined as in (2) and (16). This will be compared with
the hybridised upwind scheme (with advective fluxes as in (15)), and the cell-
centered second-order upwind scheme (with advective fluxes (13)). We note
however that for cell-centered schemes, the advective fluxes (13) cannot be
straightforwardly computed near the boundary of the domain; hence, we switch
into a first-order upwind scheme near the boundary of the domain. This will
lead to a system with $N_{K}$ \+ $N_{e_{\mathrm{ext}}}$ equations and
unknowns, where $N_{e_{\mathrm{ext}}}$ is the number of boundary edges.

###  6.1. 1D test: $\epsilon$-sensitivity

We start by performing a test over a one dimensional domain $\Omega=(0,1)$.
Here, we check for the $\epsilon$-sensitivity of the schemes. This test checks
that the numerical diffusion introduced by the scheme (if any) is not too
much. In particular, if the actual solution of the problem contains a thin
boundary layer, we also expect the numerical solution to observe the same
property. Consider the ODE

$c^{\prime}(x)-\epsilon c^{\prime\prime}(x)=0\quad\mbox{for}\quad x\in(0,1),$

with Dirichlet boundary conditions

$c(0)=1,c(1)=0.$

The solution for this differential equation can be calculated exactly, and is
given by

$c(x)=\frac{1}{e^{-1/\epsilon}-1}\bigg{(}e^{1/\epsilon(x-1)}-1\bigg{)}.$

Here, the exact solution $c$ has a boundary layer, which is controlled by the
diffusion parameter $\epsilon$. A numerical scheme with a good
$\epsilon$-sensitivity should allow us to capture the boundary layer, even
when it can only be resolved up to 1 grid cell on the mesh. For this test, we
consider a mesh with 100 equidistant cells. Hence, the size of each cell is
given to be $h=0.01$. We then consider the diffusion parameters
$\epsilon=2^{-4},2^{-6},2^{-8},2^{-10}$. In Figures 6-8, left, the numerical
solutions (dashed lines) are plotted against the exact solution (solid lines)
for different values of $\epsilon$. We then plot on the right of Figures 6-8,
the pointwise error values $c(x)-\Pi_{{\mathcal{D}}}c(x)$. Since the numerical
solutions only deviate from the actual solution near the boundary layer, we
zoom in and present the plots for $x\in[0.5,1]$.

|
---|---
Figure 6.  Exact solution against numerical solution, hybridised upwind scheme. (left: solution plots, right: errors). |
---|---
Figure 7.  Exact solution against numerical solution, cell-centered second-order scheme. (left: solution plots, right: errors). |
---|---
Figure 8.  Exact solution against numerical solution, hybridised second-order
scheme. (left: solution plots, right: errors).

We start by looking at the first-order upwind scheme. We note in Figure 6,
right that the quantity $c(x)-\Pi_{{\mathcal{D}}}c(x)$ is always nonnegative.
This is expected from first-order upwind schemes due to the numerical
diffusion it introduces, which leads to the smoothening of the solution and
widening of the boundary layer. This is illustrated in Figure 6, left. We now
move on to the cell-centered second-order scheme. Here, we see a better
agreement between the numerical solution and the exact solution in the
interior of the domain. However, since we switch to a first-order upwind
scheme near the boundary of the domain, we see that the solution and error
plots in Figure 7 are similar to those in Figure 6. Now, upon looking at the
hybridised second-order scheme in Figure 8, we see that for $\epsilon\geq
2^{-8}$, the numerical solutions obtained via the hybridised second-order
scheme behave in a manner that is very similar to the exact solution. We note
however, that for $\epsilon=2^{-10}$, the numerical solution obtained from the
hybridised second-order scheme exhibits an overshoot. This is expected, since
second-order linear schemes do not guarantee monotonicity and stability of the
numerical solution. One way to resolve this is to do, as with the cell-
centered scheme, switching into a first-order upwind scheme near the boundary
of the domain. However, this will only result to numerical solutions similar
to that in Figure 7, and will not be helpful in resolving the boundary layer.
Another way to resolve the overshoot is by introducing an artificial vanishing
diffusion term. In order to do so, we consider, for the hybridised second-
order scheme, a diffusion parameter $\epsilon+h^{1.5}$ instead. This leads to
a significant improvement in the numerical results, as can be seen in Figure
9. In particular, we are now able to capture the boundary layer, with only
small overshoots (less than 5%) in the numerical solution.

|
---|---
Figure 9.  Exact solution against numerical solution, hybridised second-order
scheme with vanishing diffusion. (left: solution plots, right: errors).

To summarise, this test shows us that the numerical diffusion introduced by
the hybridised first-order upwind scheme smoothens the solution, leading to a
widening of the boundary layer. On the other hand, with the introduction of
vanishing diffusion, the numerical solution provided by the hybridised second-
order scheme is able to capture the boundary layer properly, with only minimal
overshoots.

###  6.2. 2D tests

We now proceed to tests in 2D, which will be done on the domain
$\Omega=(0,1)\times(0,1)$. The numerical tests will be done on different mesh
types, starting with regular Cartesian and triangular meshes (see Figure 10),
followed by distorted meshes, which include moved Cartesian, moved triangular,
and Kershaw type meshes [10, 12] (see Figure 11). These mesh types will be
denoted by ${\mathcal{M}}_{1},{\mathcal{M}}_{2},\dots\mathcal{M}_{5}$,
respectively. Here, the moved Cartesian and moved triangular meshes in Figure
11 are constructed following the guidelines provided in [13]. That is,
starting with a uniform Cartesian and triangular mesh as in Figure 10, if the
maximum diameter of the cells are given by $h$, then the internal nodes
$(x,y)$ are perturbed randomly by taking

$\hat{x}:=x+0.4\beta_{x}h,\quad\hat{y}:=y+0.4\beta_{y}h,$

where $\beta_{x},\beta_{y}$ are random values between $-0.5$ and $0.5$.

|
---|---
Figure 10. Mesh types: ${\mathcal{M}}_{1}$ Cartesian (left); ${\mathcal{M}}_{2}$ triangular (right).  |  |
---|---|---
Figure 11. Mesh types: ${\mathcal{M}}_{3}$ moved Cartesian (left);
${\mathcal{M}}_{4}$ moved triangular (middle); ${\mathcal{M}}_{5}$ Kershaw
(right).

For the convergence tests, we measure the relative solution error

$E_{c}:=\dfrac{\left\|\Pi_{{\mathcal{D}}_{h}}c_{h}-c\right\|_{L^{2}(\Omega)}}{\left\|c\right\|_{L^{2}(\Omega)}},$

and the relative error in the discrete gradient

$E_{g}:=\dfrac{\left\|\nabla_{{\mathcal{D}}_{h}}c_{h}-\nabla
c\right\|_{L^{2}(\Omega)}}{\left\|c\right\|_{H^{1}(\Omega)}}.$

In order to have a more detailed comparison, we also give the number of DOFs
required to implement each of the schemes for certain mesh types.

Table 1. Mesh size $h$ and DOFs for hybridised and cell-centered schemes for each refinement level $r$, mesh ${\mathcal{M}}_{1}$. $r$ | $h$ | $N_{K}$ | $N_{e}$ | DOFs(hybridised) | DOFs(cell-centered)
---|---|---|---|---|---
1 | 3.535e-01 | 16 | 40 | 56 | 32
2 | 1.767e-01 | 64 | 144 | 208 | 96
3 | 8.838e-02 | 256 | 544 | 800 | 320
4 | 4.419e-02 | 1024 | 2112 | 3136 | 1152
5 | 2.209e-02 | 4096 | 8320 | 12416 | 4352
6 | 1.104e-02 | 16384 | 33024 | 49408 | 16896

####  6.2.1. Convergence test, smooth solution

We start with performing a convergence test for the advection-diffusion
equation (1) with prescribed solution

$c(x,y)=\sin(\pi x)\sin(\pi y).$

Here, we set the velocity field $\mathbf{V}=[1,2]$ and in order to have an
anisotropic advection-dominated problem, we set the diffusion tensor

$\Lambda=\begin{bmatrix}1.5\times 10^{-4}&10^{-6}\\\
10^{-6}&10^{-8}\end{bmatrix}.$

$10^{-2}$$10^{-1}$$10^{0}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$2$h$$E_{c}$ | $10^{-2}$$10^{-1}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$1$h$$E_{g}$
---|---
Figure 12. Convergence plot, test 6.2.1. $\circ-$hybridised second-order,
$\square-$hybridised first-order, $\diamond-$cell-centered second-order, blue:
${\mathcal{M}}_{1}$, red: ${\mathcal{M}}_{3}$. The slopes of the triangles are
order $h^{2}$ and $h$ on the left, $h$ and $h^{0.5}$ on the right.

As can be seen in Figures 12 and 13, the hybridised second-order scheme
attains second-order convergence in the solution and first-order convergence
in the gradient, except on Cartesian type meshes, where the gradient converges
at a rate of approximately $h^{0.7}$. This is due to the fact that the error
in the gradient is much smaller than the error in the solution profile; once
the error in the solution is smaller than that of the gradient, we expect to
observe first-order convergence for the gradient. Firstly, we note that the
hybridised second-order scheme is an improvement over the hybridised (first-
order) upwind scheme, which has a solution profile that converges with $O(h)$.
Moreover, for all mesh types, the solution of the hybridised second-order
scheme on the second refinement level is already much better than that of the
hybridised upwind scheme on the finest mesh.

We now look at the second-order cell-centered scheme. We see that in general,
this gives an improvement over that of the hybridised first-order scheme, i.e.
the second-order cell-centered scheme provides a more accurate solution with
the same mesh size $h$. However, no significant improvement is observed on the
accuracy in terms of the gradient reconstruction. Upon comparing the cell-
centered scheme with the hybridised second-order scheme, we see that the
hybridised scheme performs better, both in terms of mesh size $h$ and in terms
of the number of DOFs needed to solve the system. In particular, for mesh type
${\mathcal{M}}_{1}$, the solution of the second-order cell-centered scheme on
the finest mesh, which requires solving a system with 16896 unknowns, gives a
relative error of 5.903e-04. However, the hybridised second-order scheme
already achieves a relative error of 7.997e-04 by solving a system of only
3136 unknowns. A similar observation in terms of the advantages (in DOFs and
mesh size) can also be made on the other meshes. We also note here that
although the discrete gradient for the cell-centered second-order scheme
converges with an order $h^{0.5}$ on ${\mathcal{M}}_{1}$, no convergence is
observed in the other meshes. This can be explained by the fact that for
advection-dominated problems, the advective fluxes are the dominant factors,
and hence, the conservation of fluxes (17b) approximately imposes that

$F_{K,\sigma}^{A}+F_{L,\sigma}^{A}=0.$

This leads to taking $c_{\sigma}$ from the upwind direction, thus making the
discrete gradient (3) equivalent to (14), which also explains why the error in
the gradient $E_{g}$ for both the hybridised upwind scheme and the cell-
centered second-order scheme are quite close to each other.

$10^{-1.5}$$10^{-1}$$10^{-0.5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$2$h$$E_{c}$ | $10^{-1.5}$$10^{-1}$$10^{-0.5}$$10^{-2}$$10^{-1}$$10^{0}$1$h$$E_{g}$
---|---
Figure 13. Convergence plot, test 6.2.1. $\circ-$hybridised second-order,
$\square-$hybridised first-order, $\diamond-$cell-centered second-order,
green: ${\mathcal{M}}_{2}$, brown: ${\mathcal{M}}_{4}$, black:
${\mathcal{M}}_{5}$. The slopes of the triangles are order $h^{2}$ and $h$ on
the left, $h$ and $h^{0.5}$ on the right.

Upon having a closer look at Figures 12 and 13, we observe that the second-
order cell-centered scheme was able to attain second-order convergence on the
regular Cartesian mesh ${\mathcal{M}}_{1}$, whilst only first-order
convergence on the other types of meshes. This can be explained more clearly
by looking at the gradient in Figures 12 and 13, right. Here, we see that the
approximate gradient (14) converges for Cartesian type meshes, whereas it does
not converge on the other types of meshes. Hence, the linear term
$\widetilde{\nabla}_{\mathcal{D}}c_{K}\cdot({\boldsymbol{x}}-{\boldsymbol{x}}_{K})$
added onto the advective flux (16) is not accurate and thus does not help
improve the convergence of the solution. This means that an improvement over
the formulation $\eqref{eq:discGrad_upwind}$ of the discrete gradient
$\widetilde{\nabla}_{\mathcal{D}}c_{K}$, such as that proposed in [2], would
be needed in order to apply a second-order cell-centered scheme on generic
meshes. This is not straightforward to implement; however, the observations
made on the Cartesian meshes ${\mathcal{M}}_{1}$ are sufficient to support the
claim that hybridised second-order schemes perform better than cell-centered
schemes.

####  6.2.2. Convergence test, solution with boundary layers

We now perform a test for a solution with a boundary layer. For this test
case, we prescribe an exact solution

$c(x,y)=\bigg{(}x-e^{\frac{2(x-1)}{\nu}}\bigg{)}\bigg{(}y^{2}-e^{\frac{3(y-1)}{\nu}}\bigg{)}.$

Here, we set the velocity field $\mathbf{V}=[2,3]$ and take
$\Lambda=\nu\mathbf{I}$, where $\mathbf{I}$ is the identity matrix. Setting
$\nu=10^{-4}$ leads to an advection-dominated problem for which the solution
is characterized by a boundary layer near the top and right side of the
domain. As in [3, 13, 14], the goal of our numerical tests is to demonstrate
that the scheme has good convergence properties and produces numerical
solutions without oscillations in a subdomain outside the boundary layer.
Hence, we measure only the errors in the subdomain $[0,0.8]\times[0,0.8]$.

$10^{-2}$$10^{-1}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$2$h$$E_{c}$ | $10^{-2}$$10^{-1}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$1$h$$E_{g}$
---|---
Figure 14. Convergence plot for the hybridised second-order scheme on
different mesh types (blue: ${\mathcal{M}}_{1}$, green: ${\mathcal{M}}_{2}$,
red: ${\mathcal{M}}_{3}$, brown: ${\mathcal{M}}_{4}$, black:
${\mathcal{M}}_{5}$), test 6.2.2. The slopes of the triangles are order
$h^{2}$ and $h$ on the left, $h$ and $h^{0.5}$ on the right.

As can be seen in Figure 14, the hybridised second-order scheme provides
numerical solutions that are second-order convergent, and gradients that are
first-order convergent, regardless of the mesh. Comparison with the hybridised
upwind scheme and the cell-centered second-order scheme yielded similar
results as test 6.2.1, and are no longer presented here.

We now study the shock-capturing behavior of the hybridised second-order
scheme by plotting the numerical solution in Figure 15. We only plot on mesh
types ${\mathcal{M}}_{1}$ and ${\mathcal{M}}_{3}$ (one regular, one
irregular), as the numerical solution on other types of meshes exhibit a
similar behavior.

|
---|---
Figure 15. Solution profile, hybridised second-order scheme, test 6.2.2.
(left: ${\mathcal{M}}_{1}$, right: ${\mathcal{M}}_{3}$).

As can be seen, non-physical oscillations develop at the boundary of the
domain; this is expected for second-order schemes which are not flux-limited
or total variation diminishing. Moreover, these oscillations are worse on
distorted meshes than on regular meshes. In comparison, the hybridised first-
order upwind scheme provides solutions which are bounded between 0 and 1,
without any non-physical oscillations (see Figure 16).

|
---|---
Figure 16. Solution profile, hybridised upwind scheme, test 6.2.2. (left:
${\mathcal{M}}_{1}$, right: ${\mathcal{M}}_{3}$).

In order to mitigate the non-physical oscillations for the hybridised second-
order schemes, we want, as in test 6.1, to introduce an artificial vanishing
diffusion term. To this end, employ an idea that is similar to that in [7].
That is, if $\mathbf{V}_{K}$ is an approximation of $\mathbf{V}$ at cell $K$
and if $\Lambda_{K}$ is an approximation of $\Lambda$ at cell $K$ with
diagonalisation $\Lambda_{K}=U_{K}^{\prime}D_{K}U_{K}$, we consider

$\widetilde{\Lambda}_{K}=U_{K}^{\prime}(D_{K}+|\mathbf{V}_{K}|h^{1.5})U_{K},$
(36)

and use $\widetilde{\Lambda}_{K}$ instead in the definition of the diffusive
fluxes (2). As can be seen in Figure 17, introducing a vanishing diffusion
term and using the modified diffusion tensor (36) allows us to mitigate the
non-physical oscillations. Moreover, the solutions here are now bounded
between 0 and 1. We note, however, that since the artificial diffusion
vanishes at a rate of $h^{1.5}$, our scheme reduces to order 1.5; this is not
the optimal order 2 convergence, but still offers an improvement over the
hybridised first-order scheme. Alternatively, the use of nonlinear flux-
limited schemes, as in [4, 14] can also mitigate the non-physical
oscillations, whilst preserving second-order accuracy.

|
---|---
Figure 17. Solution profile, hybridised second-order scheme with vanishing
diffusion, test 6.2.2. (left: ${\mathcal{M}}_{1}$, right:
${\mathcal{M}}_{3}$).

####  6.2.3. Strongly anisotropic heterogeneous and convection-dominated case

Finally, we present a numerical test which involves a strongly heterogeneous
and anisotropic diffusion tensor, as described in [3, 8]. For this test, an
exact analytic solution is not available, so we comment on the qualitative
properties of the numerical solution. Here, homogeneous Dirichlet boundary
conditions are imposed, and we use a source term
$f(x,y)=10^{-2}\exp(-(r-0.35)^{2}/0.005)$, where
$r^{2}=(x-0.5)^{2}+(y-0.5)^{2}$. The diffusion tensor is piecewise constant,
defined in the following subdomains:
$\Omega_{1}=(0,2/3)\times(0,2/3),\Omega_{2}=(2/3,1)\times(0,2/3),\Omega_{3}=(2/3,1)\times(2/3,1),\Omega_{4}=(0,2/3)\times(2/3,1)$,
with

$\Lambda=\begin{bmatrix}10^{-6}&0\\\ 0&1\end{bmatrix}\qquad\mbox{ in
}\Omega_{1}\mbox{ and }\Omega_{3},$

and

$\Lambda=\begin{bmatrix}1&0\\\ 0&10^{-6}\end{bmatrix}\qquad\mbox{ in
}\Omega_{2}\mbox{ and }\Omega_{4}.$

The velocity field considered is
$\mathbf{V}=(40x(2y-1)(x-1),-40y(2x-1)(y-1))^{T}$, which simulates a
counterclockwise rotation. Figure 19 shows the numerical solution obtained
from the hybridised second-order scheme on a regular and distorted mesh
(${\mathcal{M}}_{1}$ and ${\mathcal{M}}_{3}$), respectively. Here, the
distorted mesh ${\mathcal{M}}_{3}$ is modified so that it matches the
discontinuities (see Figure 18). The numerical solutions obtained on other
types of meshes exhibit similar behaviors.

|
---|---
Figure 18. Diffusion tensor and mesh, test 6.2.3.  |
---|---
Figure 19. Solution profile, test 6.2.3. (left: ${\mathcal{M}}_{1}$, right:
${\mathcal{M}}_{3}$).

Here, we observe maximum values of $7.3\times 10^{-4}$ and $7.6\times 10^{-4}$
for the regular and distorted meshes, respectively. The scheme works well on
regular meshes in the sense that the solution profile is similar to those
observed in [3, 8]. Also, the maximum value of $7.3\times 10^{-4}$ is very
close to $6.9\times 10^{-4}$ in the literature. However, it can be seen in
Figure 19, right, that some spurious oscillations are present on the distorted
mesh. As with test 6.2.2, such a problem is not encountered when using a
hybridised upwind scheme. These oscillations can be mitigated by either
refining the mesh, or using the modified diffusion tensor (36), as seen in
Figure 20. The main advantage, however, of using the modified diffusion tensor
over mesh refinement is that a better quality of the solution profile is
obtained without having to introduce additional DOFs for solving the system.

|
---|---
Figure 20. Solution profile, test 6.2.3, ${\mathcal{M}}_{3}$. (left: vanishing
diffusion, right: mesh refinement).

##  7\. Conclusion

In this work, we proposed a fully local hybridised second-order finite volume
scheme for advection-diffusion equations. We then presented a convergence
proof for these hybridised second-order schemes, which also cover flux-limited
variants. Numerical results were then presented to compare the hybridised
second-order scheme with the hybridised upwind scheme, and the classical cell-
centered second-order scheme. Firstly, we note that the hybridised second-
order scheme provided solutions which are second-order convergent, and
gradients which are first-order convergent. It was also shown that the
hybridised second-order scheme had a good $\epsilon$-sensitivity. Upon
comparison with the cell-centered scheme, the hybridised scheme can achieve
the same level or even better accuracy on much coarser meshes. Moreover, the
hybridised scheme can straightforwardly be extended onto generic meshes,
whereas further improvement on the approximate gradient (14) is needed in
order to apply the cell-centered schemes onto generic meshes, otherwise the
convergence reduces to first order. Another advantage of the hybridised scheme
over the cell-centered second-order schemes is that the stencil only depends
on local values, and does not need information from neighboring cells, which
allows it to be straightforwardly implemented near the boundaries of the
domain. Moreover, static condensation can be employed in order to implement
the hybridised scheme efficiently. Comparison with hybridised upwind schemes
also shows the advantage of the hybridised second-order schemes in terms of
the convergence of the solution and the gradient. However, in some instances,
for solutions with boundary layers, the solution from the hybridised second-
order scheme exhibited non-physical oscillations. This is expected from linear
second-order schemes. In comparison, the hybridised upwind scheme does not
encounter such problems. In order to mitigate these non-physical oscillations,
we propose an idea which involves introducing artificial vanishing diffusion
to the hybridised second-order scheme. This reduces the order of convergence
to 1.5, but still offers an improvement over the hybridised upwind scheme. One
prospect for future work would involve trying to determine when artificial
diffusion is needed, and locate on which cells it needs to be introduced, so
that the second-order accuracy is preserved over regions for which the
solution is smooth. Another avenue for future work would involve extending
these ideas to hybrid high order schemes [15], with the aim of obtaining order
$h^{(k+2)}$ estimates for polynomials of degree $k$ in advection-dominated
regimes. We also aim to extend these ideas to time-dependent advection
diffusion equations.

##  8\. Acknowledgements

The author would want to thank Prof. Jérôme Droniou and Prof. Barry Koren for
the discussions and advice, which helped improve the presentation of the
paper.

##  9\. Appendix

In this section, we present without proof two lemmas from [3, 9] which are
used in the convergence proof of Theorem 5.5.

######  Lemma 9.1 (Discrete Sobolev inequality).

Let ${\mathcal{M}}_{h}$ be an admissible discretisation of $\Omega$ satisfying
assumptions (MR1) and (MR2). Let $\theta>0$ and such that
$\theta<\frac{d_{K,\sigma}}{d_{K^{\prime},\sigma}}<\theta^{-1}$ for all
$\sigma\in{\mathcal{E}}_{h,int}$. Let $r=\frac{2d}{d-2}$ if $d>2$ and
$r<\infty$ if $d=2$. Then there exists a real positive constant $C$ that only
depends on $\Omega,\theta$ and $r$ such that, for all $q_{h}\in
X_{{\mathcal{D}}_{h}}$, we have $\left\|q_{h}\right\|_{L^{r}(\Omega)}\leq
C\left\|q_{h}\right\|_{1,{\mathcal{D}}_{h}}$ .

######  Lemma 9.2 (Discrete Rellisch theorem).

Let $\Lambda$ be a diffusion tensor satisfying hypothesis (A1). Let
$({\mathcal{M}}_{h})_{h\rightarrow 0}$ be a family of admissible
discretizations of $\Omega$ with mesh size $h$ tending to 0 and satisfying the
regularity assumptions (MR1) and (MR2). Let $c_{h}\in X_{{\mathcal{D}}_{h}}$
be a numerical scalar field such that
$\left\|c_{h}\right\|_{1,{\mathcal{D}}_{h}}$ remains bounded as $h\rightarrow
0$. Then there exists a scalar field $c\in H^{1}(\Omega)$ such that, up to a
subsequence as $h\rightarrow 0$, the following hold:

  1. (1)

$c_{h}\rightarrow c$ in $L^{r}(\Omega)$ for all $r<\frac{2d}{d-2}$;

  2. (2)

$\overline{\nabla}_{{\mathcal{D}}_{h}}c_{h}\rightarrow\nabla c$ weakly in
$L^{2}(\Omega)^{d}$.

## References

  * [1] Arnold, D. N. and Brezzi, F.  Mixed and nonconforming finite element methods : implementation, postprocessing and error estimates.  ESAIM: M2AN, 19(1):7–32, 1985.
  * [2] T. Barth and D. Jespersen.  The design and application of upwind schemes on unstructured meshes.  In 27th Aerospace Sciences Meeting. 1989.
  * [3] L. Beirao da Veiga, J. Droniou, and G. Manzini.  A unified approach for handling convection terms in finite volumes and mimetic discretization methods for elliptic problems.  IMA Journal of Numerical Analysis, 31(4):1357–1401, 2011.
  * [4] E. Bertolazzi and G. Manzini.  A second-order maximum principle preserving finite volume method for steady convection-diffusion problems.  SIAM Journal on Numerical Analysis, 43(5):2172–2199, 2005.
  * [5] J. Droniou, R. Eymard, T. Gallouët, C. Guichard, and R. Herbin.  The gradient discretisation method, volume 82 of Mathematics & Applications.  Springer, 2018.
  * [6] J. Droniou, R. Eymard, T. Gallouët, and R. Herbin.  A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods.  Math. Models Methods Appl. Sci., 20(2):265–295, 2010.
  * [7] J. Droniou, R. Eymard, A. Prignet, and K. S. Talbot.  Unified convergence analysis of numerical schemes for a miscible displacement problem.  Found. Comput. Math., 19(2):333–374, Apr. 2019.
  * [8] A. Ern, A. F. Stephansen, and P. Zunino.  A discontinuous Galerkin method with weighted averages for advection―diffusion equations with locally small and anisotropic diffusivity.  IMA J. Numer. Anal., 29:235–256, 2009.
  * [9] R. Eymard, T. Gallouët, and R. Herbin.  Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes sushi: a scheme using stabilization and hybrid interfaces.  IMA Journal of Numerical Analysis, 30(4):1009–1043, 2010.
  * [10] R. Herbin and F. Hubert.  Benchmark on discretization schemes for anisotropic diffusion problems on general grids.  In Finite volumes for complex applications V, pages 659–692. ISTE, London, 2008.
  * [11] M. Holmes.  Introduction to numerical methods in differential equations.  2006\.
  * [12] D. S. Kershaw.  Differencing of the diffusion equation in lagrangian hydrodynamic codes.  Journal of Computational Physics, 39(2):375 – 395, 1981.
  * [13] K. Lipnikov, D. Svyatskiy, and Y. Vassilevski.  A monotone finite volume method for advection–diffusion equations on unstructured polygonal meshes.  Journal of Computational Physics, 229(11):4017 – 4032, 2010.
  * [14] G. Manzini and A. Russo.  A finite volume method for advection–diffusion problems in convection-dominated regimes.  Computer Methods in Applied Mechanics and Engineering, 197(13):1242 – 1261, 2008.
  * [15] D. A. D. Pietro, J. Droniou, and A. Ern.  A discontinuous-skeletal method for advection-diffusion-reaction on general meshes.  SIAM J. Numerical Analysis, 53:2135–2157, 2015.
  * [16] H. S. Price, R. S. Varga, and J. E. Warren.  Application of oscillation matrices to diffusion-convection equations.  Journal of Mathematics and Physics, 45(1-4):301–311, 1966.