# A skeletal high-order structure preserving scheme for advection-diffusion equations

Julien Moatti

► To cite this version:

Julien Moatti. A skeletal high-order structure preserving scheme for advection-diffusion equations. Finite Volumes for Complex Applications X, Oct 2023, Strasbourg, France. pp.345-354, 10.1007/978-3-031-40864-9_29. hal-04036599v2

HAL Id: hal-04036599

https://hal.science/hal-04036599v2

Submitted on 20 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A skeletal high-order structure preserving scheme for advection-diffusion equations

Julien Moatti

Inria, Univ. Lille, CNRS, UMR 8524 - Laboratoire Paul Painlevé, F-59000 Lille,
France
julien.moatti@inria.fr

**Abstract.** We introduce a nonlinear structure preserving high-order scheme for anisotropic advection-diffusion equations. This scheme, based on Hybrid High-Order methods, can handle general meshes. It also has an entropy structure, and preserves the positivity of the solution. We present some numerical simulations showing that the scheme converges at the expected order, while preserving positivity and long-time behaviour.

**Keywords:** Anisotropic advection-diffusion equations, general meshes, high-order schemes, structure preserving methods.

## 1 Motivations and context

We are interested in the discretisation of a linear advection-diffusion equation on general meshes with a high-order scheme. Let $\Omega$ be an open, bounded, connected polytopal subset of $\mathbb{R}^d$, $d \in \{2, 3\}$. We consider the following problem with homogeneous Neumann boundary conditions: find $u: \mathbb{R}_+ \times \Omega \to \mathbb{R}$ solution to

$$ egin{cases} 
\partial_t u - \operatorname{div}(\Lambda(\nabla u + u\nabla\phi)) = 0 & \text{in } \mathbb{R}_+ \times \Omega, \\
\Lambda(\nabla u + u\nabla\phi) \cdot n = 0 & \text{on } \mathbb{R}_+ \times \partial\Omega, \\
u(0, \cdot) = u^{in} & \text{in } \Omega,
\end{cases} 
\quad (1) $$

where $n$ is the unit normal vector to $\partial\Omega$ pointing outwards from $\Omega$. We assume that the data satisfy: (i) $\Lambda \in L^\infty(\Omega; \mathbb{R}^{d \times d})$ is a uniformly elliptic diffusion tensor: there exists $\lambda_b > 0$ such that, for a.e. $x$ in $\Omega$, $\Lambda(x)\xi \cdot \xi \ge \lambda_b|\xi|^2$ for all $\xi \in \mathbb{R}^d$; (ii) $\phi \in C^1(\bar{\Omega})$ is a regular potential; (iii) $u^{in} \in L^1(\Omega)$ is a non-negative initial datum, such that $\int_{\Omega} u^{in} \log(u^{in}) < \infty$. The solutions to (1) enjoy some specific and well-known properties. First the mass is preserved along time, i.e. for almost every $t > 0$, $\int_{\Omega} u(t) = \int_{\Omega} u^{in} = M$ where $M > 0$ is the initial mass. Second, the solution is positive for $t > 0$. Last, the solution has a specific long-time behaviour: it converges exponentially fast when $t \to \infty$ towards the thermal equilibrium $u^\infty$, solution to the stationary problem associated to (1), defined as $u^\infty = \frac{M}{\int_{\Omega} e^{-\phi}} e^{-\phi}$.

In order to get a reliable numerical approximation of such problems, one has to preserve these structural properties at the discrete level. It is well-known
that two-point finite volume methods are structure preserving (see [2] for the long-time behaviour), but these methods can only be used for isotropic tensors on meshes satisfying some orthogonality conditions. On the other hand, finite volume methods (using auxiliary unknowns) for anisotropic problems on general meshes were introduced in the past twenty years, but none of these linear methods preserve the positivity of the solutions (see [7]). A possible alternative was proposed in [1], with the introduction and analysis of a nonlinear positivity preserving Vertex Approximate Gradient VAG scheme. Following these ideas, a nonlinear Hybrid Finite Volume (HFV) scheme was designed in [3].

All the schemes discussed above are at most of order two in space (in $L^2$ norm). The aim of this paper is to introduce a high-order scheme preserving the three structural properties discussed above. Since the HFV method coincides with the low-order version of the Hybrid High-Order (HHO) scheme introduced in [5], we propose an HHO generalisation of the scheme introduced in [3]. Numerical results indicate that this scheme offers a better efficiency in terms of computational cost than low order schemes.

# 2 Discrete setting and scheme

## 2.1 Mesh

We define a discretisation of $\Omega$ as a pair $\mathcal{D} = (\mathcal{M}, \mathcal{E})$, where:

- the mesh $\mathcal{M}$ is a partition of $\Omega$ into cells, i.e., a finite family of nonempty disjoint open polytopal subsets $K$ of $\Omega$ such that $\bar{\Omega} = \bigcup_{K \in \mathcal{M}} \bar{K}$,

- the set of faces $\mathcal{E}$ is a partition of the mesh skeleton $\bigcup_{K \in \mathcal{M}} \partial K$ into faces $\sigma$ which are subsets contained in hyperplanes of $\bar{\Omega}$. We denote by $\mathcal{E}_K$ the set of faces of the cell $K$, and we define $n_{K,\sigma} \in \mathbb{R}^d$ as the unit normal vector to $\sigma$ pointing outwards from $K$.

The diameter of a subset $X \subset \bar{\Omega}$ is denoted by $h_X = \sup\{|x - y| \mid (x, y) \in X^2\}$. We define the mesh size of $\mathcal{D}$ as $h_{\mathcal{D}} = \sup\{h_K \mid K \in \mathcal{M}\}$. We refer to [6, Section 1.1] for more detailed statements about the mesh and its regularity.

## 2.2 Polynomials, discrete unknowns and discrete operators

In the following, $k$ is a fixed non-negative integer. First, we introduce polynomial spaces on a subset $X \subset \bar{\Omega}: \mathbb{P}^k(X)$ and $\mathbb{P}^k(X)^d$ denote respectively the spaces of polynomial functions $X \to \mathbb{R}$ and polynomial vector fields $X \to \mathbb{R}^d$ of degree at most $k$. Given $Y \subset \bar{X}$, we also define the $L^2$-projector $\Pi_Y^k: C^0(\bar{X}) \to \mathbb{P}^k(Y)$ by the relation $\forall w \in \mathbb{P}^k(Y), \int_Y \Pi_Y^k(v)w = \int_Y vw$.

We now introduce the set of discrete unknowns corresponding to the mixed-order HHO method [4,6], with face unknowns of degree $k$ and (enriched) cells unknowns of degree $k + 1$:

$$ \underline{V}_{\mathcal{D}}^{k,k+1} = \left\{ \underline{v}_{\mathcal{D}} = ((v_K)_{K \in \mathcal{M}}, (v_{\sigma})_{\sigma \in \mathcal{E}}) \mid \begin{array}{l} \forall K \in \mathcal{M}, v_K \in \mathbb{P}^{k+1}(K) \\ \forall \sigma \in \mathcal{E}, v_{\sigma} \in \mathbb{P}^k(\sigma) \end{array} \right\}. $$
Given a cell $K \in \mathcal{M}$, we let $\underline{V}_K^{k,k+1} = \mathbb{P}^{k+1}(K) \times \prod_{\sigma \in \mathcal{E}_K} \mathbb{P}^k(\sigma)$ be the restriction of $\underline{V}_D^{k,k+1}$ to $K$, and for any generic discrete unknown $\underline{v}_D \in \underline{V}_D^{k,k+1}$ we denote by $\underline{v}_K = (v_K, (v_\sigma)_{\sigma \in \mathcal{E}_K}) \in \underline{V}_K^{k,k+1}$ its local restriction to the cell $K$. Given any $\underline{v}_D \in \underline{V}_D^{k,k+1}$, we associate two piecewise polynomial functions $v_\mathcal{M} : \Omega \to \mathbb{R}$ and $v_\mathcal{E} : \bigcup_{K \in \mathcal{M}} \partial K \to \mathbb{R}$ such that

$$v_{\mathcal{M}|K} = v_K \text{ for all } K \in \mathcal{M} \text{ and } v_{\mathcal{E}|\sigma} = v_\sigma \text{ for all } \sigma \in \mathcal{E}.$$

We also introduce $\underline{1}_D \in \underline{V}_D^{k,k+1}$ the discrete element such that $1_K = 1$ for any cell $K \in \mathcal{M}$ and $1_\sigma = 1$ for any face $\sigma \in \mathcal{E}$.

Now, given a cell $K \in \mathcal{M}$, we define a local discrete gradient operator $G_K^k : \underline{V}_K^{k,k+1} \to \mathbb{P}^k(K)^d$ such that, for any $\underline{v}_K \in \underline{V}_K^{k,k+1}$, $G_K^k(\underline{v}_K)$ satisfies

$$\int_K G_K^k(\underline{v}_K) \cdot \tau = \int_K \nabla v_K \cdot \tau + \sum_{\sigma \in \mathcal{E}_K} \int_\sigma (v_\sigma - v_K) \tau \cdot n_{K,\sigma} \quad \forall \tau \in \mathbb{P}^k(K)^d. \quad (2)$$

For any face $\sigma \in \mathcal{E}_K$, we also define the jump operator $J_{K,\sigma} : \underline{V}_K^{k,k+1} \to \mathbb{P}^k(\sigma)$ by

$$J_{K,\sigma}(\underline{v}_K) = \Pi_\sigma^k(v_K) - v_\sigma. \quad (3)$$

## 2.3 Scheme

Following the ideas from [1,3] our scheme relies on a nonlinear reformulation of Problem (1). To do so, we introduce the logarithm potential $\ell = \log(u)$ and the quasi-Fermi potential $w = \ell + \phi$. At least formally, one has the following relation:

$$\nabla u + u \nabla \phi = u \nabla (\log(u) + \phi) = e^\ell \nabla w. \quad (4)$$

The scheme relies on this formulation. We will discretise the potentials as polynomials, i.e. approximate $\ell$ and $w$ as discrete unknowns in $\underline{V}_D^{k,k+1}$. Then, mimicking the relation $u = e^\ell$, we will reconstruct the density thus ensuring its positivity. Therefore, a solution $(\underline{\ell}_D^n)_{n \ge 1}$ to the scheme (9) corresponds to an approximation of the logarithms of the solution $u$ (density).

More specifically, for a given discretisation $\underline{\ell}_D \in \underline{V}_D^{k,k+1}$ of the potential $\ell$, one associates a discrete density $\underline{u}_D = (u_\mathcal{M}, u_\mathcal{E})$ defined as a pair of piecewise smooth functions where $u_\mathcal{M} : \Omega \to \mathbb{R}$ corresponds to the cells unknowns and $u_\mathcal{E} : \bigcup_{K \in \mathcal{M}} \partial K \to \mathbb{R}$ corresponds to the face unknowns, defined as

$$u_\mathcal{M} = \exp(\underline{\ell}_\mathcal{M}) \text{ and } u_\mathcal{E} = \exp(\underline{\ell}_\mathcal{E}). \quad (5)$$

Note that a discrete density $\underline{u}_D$ is not a collection of polynomials (which is highlighted by the use of the wave under $u$), but it enjoys positivity, both on cells and faces, since it is defined as the exponential of real functions.

Our scheme is based on local contributions on cells, split into a consistent term and a stabilisation term. Given $K \in \mathcal{M}$ and $\eta_l > 0$, the classical discrete counterpart of $(w, v) \mapsto \int_K \Lambda \nabla w \cdot \nabla v$ is the bilinear form (see [4, Section 3.2.1])

$$a_K : (\underline{w}_K, \underline{v}_K) \mapsto \int_K \Lambda G_K^k(\underline{w}_K) \cdot G_K^k(\underline{v}_K) + \eta_l \sum_{\sigma \in \mathcal{E}_K} \frac{\Lambda_{K,\sigma}}{h_\sigma} \int_\sigma J_{K,\sigma}(\underline{w}_K) J_{K,\sigma}(\underline{v}_K),$$
where $\Lambda_{K\sigma} = \|\Lambda|_K n_{K\sigma} \cdot n_{K\sigma}\|_{L^\infty(\sigma)}$. Similarly, given $\eta_{nl} > 0$, we define a local discretisation of $(\ell, w, v) \mapsto \int_K e^\ell \Lambda \nabla w \cdot \nabla v$ as a sum of nonlinear consistent (6a) and stabilisation (6b) contributions:

$$C_K(\underline{\ell}_K, \underline{w}_K, \underline{v}_K) = \int_K e^{\ell_K} \Lambda G_K^k(\underline{w}_K) \cdot G_K^k(\underline{v}_K), \quad (6a)$$

$$S_K(\underline{\ell}_K, \underline{w}_K, \underline{v}_K) = \eta_{nl} \sum_{\sigma \in \mathcal{E}_K} \frac{\Lambda_{K\sigma}}{h_\sigma} \int_\sigma \frac{e^{\Pi_\sigma^k(\ell_K)} + e^{\ell_\sigma}}{2} J_{K,\sigma}(\underline{w}_K) J_{K,\sigma}(\underline{v}_K). \quad (6b)$$

We can now define a local application $\mathcal{T}_K : \underline{V}_K^{k,k+1} \times \underline{V}_K^{k,k+1} \times \underline{V}_K^{k,k+1} \to \mathbb{R}$ by
$$ \mathcal{T}_K(\underline{\ell}_K, \underline{w}_K, \underline{v}_K) = C_K(\underline{\ell}_K, \underline{w}_K, \underline{v}_K) + S_K(\underline{\ell}_K, \underline{w}_K, \underline{v}_K) + \varepsilon h_K^{k+2} a_K(\underline{w}_K, \underline{v}_K), \quad (7) $$

where $\varepsilon$ is a non-negative parameter. At the global level, we define $\mathcal{T}_D : \underline{V}_D^{k,k+1} \times \underline{V}_D^{k,k+1} \times \underline{V}_D^{k,k+1} \to \mathbb{R}$ by summing the local contributions:

$$\mathcal{T}_D(\underline{\ell}_D, \underline{w}_D, \underline{v}_D) = \sum_{K \in \mathcal{M}} \mathcal{T}_K(\underline{\ell}_K, \underline{w}_K, \underline{v}_K). \quad (8)$$

We let $\underline{\phi}_D \in \underline{V}_D^{k,k+1}$ be the interpolate of $\phi$: for any $K \in \mathcal{M}$, $\phi_K = \Pi_K^{k+1}(\phi)$ and for all $\sigma \in \mathcal{E}$, $\phi_\sigma = \Pi_\sigma^k(\phi)$. Now, using a backward Euler discretisation in time with time step $\Delta t > 0$, we introduce the following scheme for (1):

find $(\underline{\ell}_D^n)_{n \ge 1} \in (\underline{V}_D^{k,k+1})^{\mathbb{N}^*}$ such that

$$\left\{ \int_{\Omega} \frac{u_{\mathcal{M}}^{n+1} - u_{\mathcal{M}}^n}{\Delta t} v_{\mathcal{M}} = - \mathcal{T}_D(\underline{\ell}_D^{n+1}, \underline{\ell}_D^{n+1} + \underline{\phi}_D, \underline{v}_D) \quad \forall \underline{v}_D \in \underline{V}_D^{k,k+1}, \quad (9a)$$

$$u_K^0 = u^{in}|_K \quad \forall K \in \mathcal{M}. \quad (9b)$$

Given a solution $(\underline{\ell}_D^n)_{n \ge 1}$ to the scheme (9), as discussed above, we associate a sequence of positive discrete densities $(u_D^n)_{n \ge 1}$.

**Remark 1 (Parameter $\varepsilon$).** Note that $\mathcal{T}_D$ is to be understood as a discretisation of $(\ell, w, v) \mapsto \int_\Omega (e^\ell + \varepsilon) \Lambda \nabla w \cdot \nabla v$, with $\varepsilon \sim \varepsilon h_D^{k+2}$ a small parameter. The $\varepsilon$ perturbation is used in order to show the existence result of Proposition 2 and can be seen as a kind of stabilisation. The scaling factor $h_K^{k+2}$ in (7) is used to get the expected order of convergence. In practice, numerical results for $\varepsilon = 1$ and $\varepsilon = 0$ are almost the same. The influence of this term will be investigated in future works.

We define the discrete thermal equilibrium as $u_D^\infty = (\rho e^{-\phi_\mathcal{M}}, \rho e^{-\phi_\varepsilon})$, with $\rho = M/\int_\Omega e^{-\phi_\mathcal{M}}$. One can show that $u_D^\infty$ (and the associated logarithm potential $\underline{\ell}_D^\infty \in \underline{V}_D^{k,k+1}$) is the only stationary solution to (9) with mass $M$.

## 3 Main features of the scheme

In this section, we present some results regarding the analysis of the scheme (9). Given $\underline{\ell}_D \in \underline{V}_D^{k,k+1}$ a discrete logarithm, we associate a discrete quasi-Fermi
potential defined as $\underline{w}_D = \underline{\ell}_D + \underline{\phi}_D - \log(\rho)\underline{1}_D$. By definition of $\rho$, one has $w_M = \log(\frac{u_M}{u_M^\infty})$. Note that, for any $(\underline{\ell}_D, \underline{v}_D) \in \underline{V}_D^{k,k+1} \times \underline{V}_D^{k,k+1}$, we have $\mathcal{T}_D(\underline{\ell}_D, \underline{\ell}_D + \underline{\phi}_D, \underline{v}_D) = \mathcal{T}_D(\underline{\ell}_D, \underline{w}_D, \underline{v}_D)$. We now state our fundamental a priori results.

**Proposition 1 (Fundamental a priori relations)** Let $(\underline{\ell}_D^n)_{n \ge 1}$ be a solution to the scheme (9), and $(\underline{w}_D^n)_{n \ge 1}$ be the associated reconstructed discrete density. Then, the following a priori results hold:

(i) the mass is preserved along time: $\forall n \in \mathbb{N}^*$, $\int_{\Omega} u_M^n = \int_{\Omega} u^{in} = M$,

(ii) a discrete entropy/dissipation relation holds: $\forall n \in \mathbb{N}$, $\frac{\mathbb{E}^{n+1} - \mathbb{E}^n}{\Delta t} \le -\mathbb{D}^{n+1}$, where the discrete entropy and dissipation are defined by $\mathbb{E}^n = \int_{\Omega} u_M^\infty \Phi_1(\frac{u_M^n}{u_M^\infty})$ and $\mathbb{D}^n = \mathcal{T}_D(\underline{\ell}_D^n, \underline{w}_D^n, \underline{w}_D^n) \ge 0$ with $\Phi_1: s \mapsto s \log(s)-s+1$ (and $\Phi_1(0) = 1$).

*Proof*. Using $\underline{1}_D$ as a test function in (9a), alongside with (9b), we get the mass conservation identity (i). To get (ii), we test (9a) with $\underline{w}_D^{n+1}$, and we use the convexity of $\Phi_1$ alongside with the expression of $w_M^{n+1}$.

Note that the previous results hold for any $\varepsilon \ge 0$. Following the ideas of [1,3], the entropy/dissipation relation should allow one to analyse the long-time behaviour of the discrete solutions and to get convergence results. These aspects will be the topics of future works. We now state an existence result, which holds only for positive $\varepsilon$. The proof follows the strategy used in [3].

**Proposition 2 (Existence of solutions)** Assume that the stabilisation parameter $\varepsilon$ in (7) is positive. Then, there exists at least one solution $(\underline{\ell}_D^n)_{n \ge 1}$ to the scheme (9). The associated densities $(\underline{w}_D^n)_{n \ge 1}$ are positive functions.

# 4 Numerical results

The numerical scheme (9) requires to solve a nonlinear system of equations at each time step. To do so, we use a Newton method, with an adaptive time stepping strategy: if the Newton method does not converge, we try to compute the solution for a smaller time step $0.5 \times \Delta t$. If the method converges, we use for the subsequent time step the value $2 \times \Delta t$. The maximal time step allowed is the initial time step. Each time a linear system has to be solved we perform a static condensation (see [6, Appendix B.3.2]) in order to eliminate (locally) the cell unknowns. Note that the local computations are not implemented in parallel, but only sequentially. In the sequel, we use the following stabilisation parameters: $\varepsilon = \eta_{nl} = \eta_l = 1$.

The tests considered below (on $\Omega = ]0, 1^2]$) are the same as in [3], to which we refer for more detailed explanations and descriptions. Given a (face) degree $k$, the scheme (9) will be denoted by nlhho_k, whereas the HFV scheme of [3] will be denoted by nlhfv. Note that nlhho_0 hinges on affine cell unknowns, whereas the cell unknowns of nlhfv are constant: these two schemes hence do not coincide, and nlhho_0 is expected to be more costly.
## 4.1 Proof of concept: convergence order and efficiency

Here, we are interested in the convergence of the scheme when $(h_D, \Delta t) \to (0, 0)$. To do so, we set the advective potential and diffusion tensor as $\phi(x, y) = -x$ and $\Lambda = \begin{pmatrix} l_x & 0 \\ 0 & 1 \end{pmatrix}$ for $l_x > 0$. The exact solution is therefore given by

$$u(t, x, y) = C_1 e^{-\alpha t + \frac{\ell}{2}} (2\pi \cos(\pi x) + \sin(\pi x)) + 2C_1 \pi e^{x - \frac{1}{2}},$$

where $C_1 > 0$ and $\alpha = l_x (\frac{1}{4} + \pi^2)$. Note that $u^{in}$ vanishes on $\{x = 1\}$, but for any $t > 0$, $u(t, \cdot) > 0$. Here, our experiments are performed using $l_x = 1$ and $C_1 = 10^{-1}$. We compute the solution on the time interval $[0, 0.1]$, and we

<figure><img src="image_2.png" /><figcaption>Fig. 1: Accuracy of transient solutions. Relative error on triangular meshes.</figcaption></figure>

denote by $(u_D^n)_{1 \le n \le N_f}$ the corresponding discrete density. Then, we compute the relative $L_t^2(L_x^2)$ error on the solution and on the gradient of the solution, defined as

$$\frac{\sqrt{\sum_{n=1}^{N_f} \delta t^n \|u_M^n - u(t^n, \cdot)\|_{L^2(\Omega)}^2}}{\|u\|_{L_t^2(L_x^2)}} \quad \text{and} \quad \frac{\sqrt{\sum_{n=1}^{N_f} \delta t^n \|\mathcal{G}_M(u_D^n) - \nabla u(t^n, \cdot)\|_{L^2(\Omega)}^2}}{\|\nabla u\|_{L_t^2(L_x^2)}}$$

where $\delta t^n = t^n - t^{n-1}$ and the discrete gradient $\mathcal{G}_M(u_D^n)$ is defined by mimicking the continuous relation $\nabla u = e^\ell \nabla \ell$ as a piecewise continuous function satisfying $\mathcal{G}_M(u_D)_K = \exp(\ell_K) G_K^k(\ell_K)$ on $K \in \mathcal{M}$. The $L^2$ norms are computed using quadrature formulas of order $2k+5$. Note that, with the chosen definitions, we do not take into account the time $t = 0$. To plot the error graphs, we do simulations on a triangular mesh family $(D_i)_{1 \le i \le 5}$, such that $h_{D_i}/h_{D_{i+1}} = 2$. Since the time discretisation is of order one, on the i-th mesh of the family, we use a time step of $\Delta t_i = \Delta t_k/2^{(i-1)}$, where $\Delta t_k = 0.05/2^{k+2}$ is the initial time step used on $D_1$.

In Figure 1, we see that the scheme, for face unknowns of degree $k$, converges at order $k+1$ in energy norm and $k+2$ in $L^2$ norm of the density. In Figure 2, we plot the errors as functions of the computing time to get the solution. It is remarkable to see that, even with a low order discretisation in time, significant efficiency gains can be reached by using a high value of $k$. The gain should be even bigger by parallelising the local computations. Of course, the use of higher order
<figure><img src="image_3.png" /><figcaption>Fig. 2: **Accuracy vs. computational cost.** Relative errors on triangular meshes.</figcaption></figure>

time-stepping methods should also lead to significant gains, and this should be investigated in future works. However, the way of getting the entropy dissipation relation is currently unclear for such time discretisations.

## 4.2 Discrete long-time behaviour

We are now interested in the long-time behaviour of discrete solutions. We use the same test-case as before, but with an anisotropic tensor: we set $l_x = 10^{-2}$. We compute the solution on the time interval $[0, 350]$, with $\Delta t = 10^{-1}$, on two Kershaw meshes of sizes 0.02 and 0.006. In Figure 3, we show the evolution

<figure><img src="image_4.png" /><figcaption>Fig. 3: **Long-time behaviour of discrete solutions.** Comparison of the long-time behaviour on Kershaw meshes for $T_f = 350$ and $\Delta t = 0.1$.</figcaption></figure>

along time of the $L^1$ distance between $u_D^n$ and $u^\infty = 2C_1\pi e^{x-\frac{1}{2}}$ computed as $\int_\Omega |u_M^n - u^\infty|$. We observe the exponential convergence towards the steady-state, until some precision is reached. The rates of convergence are similar to the exact one $(\alpha)$, and do not depend on the size of the mesh.

## 4.3 Positivity

This last section is dedicated to assessing the discrete positivity preservation. We set the advection field as $\phi(x, y) = -((x - 0.4)^2 + (y - 0.6)^2)$ and the diffu-
sion tensor as $\Lambda = (^{0.8}_{0}^{0})$. For the initial data, we take $u^{in} = 10^{-3} \mathbb{1}_{B} + \mathbb{1}_{\Omega \setminus B}$, where $B$ is the Euclidean ball $\{(x, y) \in \mathbb{R}^2 | (x - 0.5)^2 + (y - 0.5)^2 \le 0.2^2\}$. We perform simulations on the time interval $[0, 5.10^{-4}]$ with $\Delta t = 10^{-5}$ on a refined tilted hexagonal-dominant mesh (4192 cells). In Table 1, we show the

<table><caption>Table 1: Positivity of discrete solutions.</caption><thead><tr><th></th><th>computing time</th><th>#resol</th><th>mincells</th><th>minfaces</th><th>mincellQN</th><th>minfaceQN</th></tr></thead><tbody><tr><td>nlhfv</td><td>1.77e+01</td><td>175</td><td>9.93e-04</td><td>7.36e-04</td><td>9.93e-04</td><td>7.36e-04</td></tr><tr><td>HMM</td><td>2.20e-01</td><td>50</td><td>-5e-03</td><td>-7.74e-02</td><td>-5e-03</td><td>-7.74e-02</td></tr><tr><td>nlhh0_0</td><td>7.17e+01</td><td>224</td><td>1.00e-03</td><td>1.01e-03</td><td>2.41e-06</td><td>1.01e-03</td></tr><tr><td>nlhh0_1</td><td>4.13e+02</td><td>248</td><td>6.65e-04</td><td>2.05e-05</td><td>1.78e-04</td><td>3.57e-08</td></tr><tr><td>nlhh0_2</td><td>1.45e+03</td><td>251</td><td>9.50e-04</td><td>5.99e-04</td><td>2.67e-07</td><td>1.06e-05</td></tr><tr><td>nlhh0_3</td><td>3.87e+03</td><td>254</td><td>9.85e-04</td><td>8.58e-04</td><td>1.10e-05</td><td>1.79e-05</td></tr></tbody></table>

minimal values reached by the schemes. The values of “mincells” are defined as $\min\{\frac{1}{|K|} \int_K u_M^n | K \in \mathcal{M}, 1 \le n \le N_f\}$, whereas “mincellQN” are the minimal values taken by the densities on the cell quadrature nodes. Analogous definitions hold for the faces. The values of “#resol” correspond to the number of linear systems solved during the computation. Note that the size of these systems depends on the value of $k$. The HMM scheme is a linear one (see [3]), therefore only one LU factorisation was performed to compute the solution, which has 90 (resp. 503) negative cell (resp. face) unknowns.

**Acknowledgements** The author thanks the anonymous reviewers for their remarks and suggestions, as well as Claire Chainais-Hillairet, Maxime Herda and Simon Lemaire for fruitful discussions about this work. The author acknowledges support by the Labex CEMPI (ANR-11-LABX-0007-01).

## References

1. Cancès C., Guichard C.: Numerical analysis of a robust free energy diminishing finite volume scheme for parabolic equations with gradient structure. Found. Comput. Math. 17, 6, 1525-1584 (2017).

2. Chainais-Hillairet C., Herda M.: Large-time behaviour of a family of finite volume schemes for boundary-driven convection-diffusion equations. IMA J. Numer. Anal. 40, 4, 2473-2504 (2020).

3. Chainais-Hillairet C., Herda M., Lemaire S., Moatti J.: Long-time behaviour of hybrid finite volume schemes for advection-diffusion equations: linear and nonlinear approaches. Numer. Math. 151, 4, 963-1016 (2022).

4. Cicuttin, M., Ern, A., Pignet, N.: Hybrid high-order methods. A primer with applications to solid mechanics. Cham: Springer (2021).

5. Di Pietro, D. A., Ern, A., Lemaire, S.: An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators. Comput. Methods Appl. Math. 14, 4, 461-472 (2014).

6. Di Pietro D. A., Droniou, J.: The hybrid high-order method for polytopal meshes. Design, analysis, and applications. Cham: Springer (2020).

7. Droniou, J.: Finite volume schemes for diffusion equations: introduction to and review of modern methods. Math. Models Methods Appl. Sci., 24, 1575-1619 (2014).