Structure preservation in high-order hybrid discretisations of potential-driven advection-diffusion: linear and nonlinear approaches

1.   Motivations and context

We are interested in the polytopal discretisation of a linear potential-driven advection-diffusion equation using high-order schemes. Our goal is to compare an exponentially fitted linear method with a nonlinear approach. Let $\Omega$ be an open, bounded, connected polytopal subset of $\mathbb{R}^d$, $d \in \{2, 3\}$, with Lipschitz boundary. We consider the following anisotropic advection-diffusion problem with homogeneous Neumann boundary conditions: find the density $u : \mathbb{R}_+ \times \Omega \to \mathbb{R}$ solution to

where $n$ is the unit normal vector to $\partial\Omega$ pointing outward from $\Omega$. We assume that the data satisfy:

$({\rm{i}})$ $\Lambda \in L^{\infty}(\Omega)^{d\times d}$ is a uniformly elliptic anisotropy tensor: there exists $\lambda_\flat > 0$ such that, for a.e. $x$ in $\Omega$, $\Lambda (x) \xi \cdot \xi \geq \lambda_\flat|\xi|^2$ for all $\xi \in \mathbb{R}^d$;

$({\rm{ii}})$ $\phi \in W^{1, \infty}(\Omega)$ is a regular potential;

$({\rm{iii}})$ $u^{0} \in L^1(\Omega)$ is a non-negative initial datum, such that $\int_\Omega u^{0} \log\left(u^{0}\right) < \infty$.

The solution to (1.1) enjoys some specific and well-known properties. First, the mass is conserved along time, i.e., for almost every $t > 0$,

where $M > 0$ is the initial mass. Second, the solution is positive:

Last, the solution converges exponentially fast when $t \to \infty$ towards the thermal equilibrium, unique steady solution to (1.1), given by

In order to compute a reliable numerical approximation of Problem (1.1), one should preserve at the discrete level the three above-listed structural properties. In practice, our final target application are drift-diffusion semiconductor models ^[37,50] (and, in particular, anisotropic ones ^[36]). In these models, the electric potential $\phi$ driving the drift is one of the unknowns of the problem, alongside with the densities of charge carriers. It is solution to a Poisson equation. At the PDE level, the thermal equilibrium is defined as the density $u^\infty > 0$ for which the flux $\Lambda(\nabla u^\infty + u^\infty \nabla \phi)$ identically vanishes in $\Omega$. This characterisation implies that the equilibrium quasi-Fermi potential $\log(u^\infty) + \phi$ shall be constant in $\Omega$ ^[43]. At the discrete level, this motivates the following definition.

Definition 1 (Preservation of the thermal equilibrium). A numerical scheme for (1.1) preserves the thermal equilibrium if the corresponding discrete equilibrium quasi-Fermi potential is constant.

Note that this definition implies that the discrete equilibrium density has to be positive. In semiconductor models discretisations, the potential $\phi$ is sought as an element of the discrete space. By Definition 1, preserving the thermal equilibrium then essentially requires to also seek $\log(u^\infty)$ as an element of the latter discrete space. For the schemes we study in this work, the precise meaning of Definition 1 will be made clear in Proposition 1 below.

In the realm of Two-Point Flux Approximation (TPFA) finite volume schemes, the so-called Scharfetter–Gummel fluxes ^[47] are precisely devised so as to preserve the thermal equilibrium. They naturally lead to linear structure-preserving discretisations of the problem (see ^[14,15]). However, TPFA methods can only be used on meshes satisfying orthogonality conditions (with respect to the inner product induced by $\Lambda$, in case $\Lambda$ is symmetric), which essentially restricts their use to isotropic problems. On the other hand, a number of finite volume methods using auxiliary unknowns has been introduced within the past twenty years or so for the discretisation of anisotropic problems on general meshes. One can cite the Discrete Duality Finite Volume (DDFV) method ^[27,38], with additional unknowns attached to a dual mesh, the Vertex Approximate Gradient (VAG) scheme ^[34], with auxiliary unknowns attached to the mesh vertices, or the Mimetic Finite Difference (MFD) and Hybrid Finite Volume (HFV) methods ^[9,33], with auxiliary unknowns attached to the mesh faces. Such methods have proved to be relevant solutions to the anisotropy issue, but none of these linear schemes preserves the positivity of the solutions (see ^[28]). A possible alternative was proposed in ^[13], with the introduction and analysis of a nonlinear positivity-preserving VAG scheme. The design and analysis of this scheme, as well as of its DDFV and HFV counterparts of ^[11,12,16], leverage the entropy structure of Problem (1.1): there exists some physically motivated quantity, called entropy, which decays along time. Reproducing this structure at the discrete level is key to get stability, convergence, and accurate time asymptotics. Other approaches to positivity preservation on general meshes have been explored in the literature. Still in the realm of finite volume methods, one can cite the works ^{[5,29,46,48,49]}. As opposed to ^[13], in which the nonlinearities are introduced at the PDE level then discretised, the latter contributions introduce nonlinearities directly at the discrete level. These nonlinearities, unfortunately, often do not lend themselves to a PDE re-interpretation, making difficult to unravel the potential discrete entropy structures hidden behind. Arbitrary-order positivity-preserving (or, more generally, discrete maximum principle preserving) methods have also been studied in the literature. In the finite element context, one can cite the seminal works ^[17,18] by Ciarlet, as well as the more recent contributions ^[2,35] (see also ^[3] for a comprehensive survey). These approaches are, however, restricted to standard meshes. In addition, only algebraic positivity can usually be enforced, that is positivity of the degrees of freedom, but not of the (piecewise polynomial) functions themselves over the domain. Weak positivity enforcement has also been explored in the Discontinuous Galerkin (DG) framework in ^[41,42]: therein, positivity is enforced on the cell averages of the piecewise polynomial solutions. Turning to pointwise positivity enforcement, let us mention in the DG context the interesting contribution ^[7]. Therein, a nonlinear scheme is introduced for the (reaction-diffusion) Fisher–KPP equation $\partial_t u - \triangle u = u(1-u)$, in which the (positive) densities are defined as $u = {\rm{e}}^\lambda$, with $\lambda$ piecewise polynomial. This scheme is developed so as to preserve the entropy structure of the PDE model. Compared to the high-order DG schemes of ^[41,42], the main improvement lies in the fact that the discrete solutions are positive everywhere. Such a feature allows for a complete analysis of the scheme, including existence, long-time behaviour, and convergence towards a semi-discretised solution. The analysis is based on the properties of a well-chosen stabilisation term, whose expression implies $L^\infty$-norms of the polynomial unknowns over the mesh faces. The results of ^[7], valid on simplicial meshes, have recently been extended (excluding the long-time behaviour) to polytopal meshes in ^[22], still in the DG context. Along the same lines, yet restricted to standard meshes, let us also cite the conforming space-time Galerkin discretisation of ^[8] for cross-diffusion systems.

From the above literature review, it is quite clear that the landscape in terms of positivity-preserving polytopal methods of arbitrarily high approximation order for advection-diffusion problems is relatively scarce. Speaking of pointwise positivity enforcement, the only existing contribution we are aware of is ^[22], in the DG context, and for reaction-diffusion equations. In the present work, our aim is to study an arbitrary-order hybrid polytopal scheme for Problem (1.1), preserving the three structural properties (1.2)–(1.4) listed above. One expected advantage of hybrid methods over DG schemes is a reduction of the number of globally coupled unknowns in the linear systems to be solved at each iteration of the Newton algorithm, which should be all the more substantial that the order of approximation increases. Our (nonlinear) scheme has been briefly introduced, and a first numerical assessment performed, in ^[44]. Our goal in the present article is twofold. First, we want to provide our approach with theoretical foundations. Second, we aim to conduct an extensive numerical validation of our method (convergence orders, efficiency, positivity, large time), including a comparison in terms of structure preservation with a similar (in the spirit) high-order linear scheme. One could indeed expect, at least in practice, that the use of a method (even linear) with sufficiently high order (and thus accuracy), could already constitute in itself a solution to positivity violation issues. The two (linear and nonlinear) methods we consider are built upon the Hybrid High-Order (HHO) technology ^[25,26], as natural extensions of the HFV schemes introduced in ^[16,45]. The linear scheme hinges on the exponential fitting strategy ^[10]. The key idea is the linear change of unknown $u = {\rm{e}}^{-\phi}\rho$, which allows one to reformulate (1.1) as an unconditionally coercive problem in the variable $\rho$. As a by-product of this reformulation, the scheme naturally preserves the thermal equilibrium. The nonlinear scheme relies on the nonlinear change of unknown $u = {\rm{e}}^\ell$ (cf. ^[13]), which is designed so as to preserve the Boltzmann entropy structure of the PDE model and, as a by-product, the positivity of solutions, the thermal equilibrium, and the long-time asymptotics. For the sake of simplicity, both schemes rely on a mixed-order HHO space: given an integer $k\geq 0$, the methods hinge on face unknowns of polynomial degree $k$, and enriched cell unknowns of polynomial degree $k+1$. The main interest of such a discrepancy in the degree between face and cell unknowns is a simplification in the design of the higher-order bulk reconstruction and of the stabilisation, resulting in turn in a simplification of the analysis. In the meantime, such a choice preserves optimal accuracy (order $k+2$ in $L^2$-norm) and frugality (the face unknowns, of degree $k$, are the only globally coupled unknowns). Since we are manipulating mixed-order spaces, following ^[20], we could also refer to our methods as HDG methods. However, we prefer naming them HHO methods for the two following reasons (cf. [19, Section 1.5.2]). First, HDG schemes are developed adopting a mixed-hybrid viewpoint, whereas we adopt here the primal HHO viewpoint. Second, the analysis of HDG methods usually hinges on specific (often simplex-based) projections, whereas our HHO analysis makes here a systematic use of $L^2$-orthogonal projectors (well-defined on polytopal cells). In any case, the two schemes we introduce in this work are new in the HDG/HHO context. Our first theoretical results, stated in Propositions 2 and 3, concern the well-posedness and long-time behaviour of the (linear) exponential fitting scheme. Regarding the nonlinear scheme, we prove the existence of (positive) solutions in Theorem 3.1. These discrete solutions are further proved to converge (in large time) in Proposition 5 towards the discrete equilibrium of the scheme. Note that we could also have compared our nonlinear scheme with the linear HHO method for advection-diffusion of ^[24], which generalises to arbitrary approximation orders the HMM scheme of ^[4] (both introduced for general advection fields). We have not pursued further in this direction, this for two reasons. First, the stability of this scheme hinges on some coercivity assumptions which constrain the variety of potentials that can be considered, and second it does not preserve the thermal equilibrium (see ^[16] in the lowest-order case).

The rest of the article is organised as follows. In Section 2, we first introduce the discrete framework, and describe the two schemes under consideration. Then, in Section 3, we discuss the main properties of the two schemes, and we provide some elements of analysis regarding the well-posedness and discrete long-time behaviours. Last, in Section 4, we discuss the implementation of the nonlinear scheme, and we assess the behaviour of the methods on various test-cases.

2.   Discrete setting and schemes

The two numerical schemes we consider in this article are based on a backward Euler discretisation in time, with uniform time step $\Delta t > 0$. The time discretisation is thus defined as $(t^n)_{n \in \mathbb{N}}$, where $t^n = n \Delta t$. Note that it is straightforward to generalise the discussion below to a variable time step. We focus in this section on space discretisation.

2.1.   Polytopal meshes and anisotropy tensor

In the vein of [23, Definition 1.4], we define a discretisation of $\Omega$ as a couple $\mathcal{D} = (\mathcal{M}, \mathcal{E})$, where:

● The mesh $\mathcal{M}$ is a partition of $\Omega$, i.e., $\mathcal{M}$ is a finite collection of disjoint, open, Lipschitz polytopes $K \subset \Omega$ with $|K|_d > 0$ (the cells) such that $\overline{\Omega} = \bigcup _ {K \in \mathcal{M} } \overline{K}$;

● The set $\mathcal{E}$ is a partition of the mesh skeleton $\partial \mathcal{M} = \bigcup_{K\in \mathcal{M}}\partial K$, i.e., $\mathcal{E}$ is a finite collection of disjoint, connected, relatively open subsets $\sigma$ of $\overline{\Omega}$ with $| \sigma|_{d{-}1} > 0$ (the faces) such that $\partial \mathcal{M} = \bigcup_{ \sigma \in \mathcal{E}} \overline{ \sigma}$. It is assumed that, for all $\sigma \in \mathcal{E}$, $\sigma$ is a Lipschitz polytopal subset of an affine hyperplane of $\mathbb{R}^d$. We assume that, for all $K\in \mathcal{M}$, there exists a subset $\mathcal{E}_K$ of $\mathcal{E}$ (the set of faces of the cell $K$) such that $\partial K = \bigcup_{ \sigma\in \mathcal{E}_K}\overline{ \sigma}$. Finally, we let $n_{K, \sigma} \in \mathbb{R}^d$ be the (constant) unit normal vector to $\sigma \in \mathcal{E}_K$ pointing outward from $K$.

The diameter of a subset $X \subset \overline{\Omega}$ is denoted by $h_X = \sup \lbrace|x-y| \mid (x, y) \in X^2 \rbrace$, and we define the size of $\mathcal{D}$ (the mesh size) as $h_ \mathcal{D} = \max_{K\in \mathcal{M}} h_K$. For further use, we also introduce the smallest cell diameter $h_\flat = \min_{K\in \mathcal{M}}h_K$ of $\mathcal{D}$.

When studying asymptotic behaviours with respect to the mesh size, one has to adopt a measure of regularity for refined families of discretisations. We classically follow [23, Definition 1.9], in which regularity for a refined mesh family is quantified by a uniform (with respect to the mesh size) parameter $\theta\in(0, 1)$, called mesh regularity parameter. This parameter measures the chunkiness of the cells, but also the diameter ratio between the cells and their faces. In what follows, to avoid the proliferation of multiplicative constants, we write $a\lesssim b$ in place of $Ca\leq b$ if $C > 0$ only depends on $\Omega$, on the mesh regularity parameter $\theta$, and (if need be) on $\Lambda$, $\phi$, and the underlying polynomial degree, but is independent of both $h_{ \mathcal{D}}$ (and $h_\flat$) and $\Delta t$.

Remark 1 (Relaxation of the mesh regularity assumptions). Upon replacing the scalings $h_ \sigma$ for $\sigma\in \mathcal{E}_K$ by $h_K$ in the stabilisations/discrete norms below, the analysis performed in this work remains valid under the (much) less stringent mesh regularity assumptions of [30, Assumption 1] (cf. also [23, Definition 1.41]). Contrary to [23, Definition 1.9], these relaxed mesh regularity assumptions allow for small faces and cells with numerous faces, as they may appear in agglomeration-based meshing.

Last, we make an additional regularity assumption on the anisotropy tensor. We assume that

2.2.   Discrete space and operators

For $q\in \mathbb{N}$, and $X$ subset of $\overline{\Omega}$ of Hausdorff dimension $1\leq l\leq d$, we let $\mathbb{P}^q(X)$ denote the vector space of $l$-variate polynomial functions $X \to \mathbb{R}$ of total degree at most $q$. We also define the $L^2(X)$-orthogonal projector $\Pi_X^q : L^1(X) \to \mathbb{P}^q(X)$ such that, given any $v \in L^1(X)$, $\Pi_X^q(v)$ is the only element in $\mathbb{P}^q(X)$ satisfying

Given any $K \in \mathcal{M}$, we also introduce the vector space $\mathbb{P}^q(K)^d$ of $d$-variate polynomial vector fields $K\to \mathbb{R}^d$ of total degree at most $q$, as well as the corresponding $L^2(K)^d$-orthogonal projector (denoted as its scalar version) $\Pi_K^q:L^1(K)^d\to \mathbb{P}^q(K)^d$. For any $\sigma \in \mathcal{E}_K$ and $v \in W^{1, 1}(K)$, we also introduce the shortcut notation

Let $k$ be a given non-negative integer. We introduce the mixed-order HHO space (see ^[19,23]), with face unknowns of degree $k$ and (enriched) cell unknowns of degree $k+1$:

Given a cell $K \in \mathcal{M}$, we let

be the restriction of $\underline{V}_ \mathcal{D}^{k}$ to $K$, and for a generic discrete element $\underline{v}_{ \mathcal{D}} \in \underline{V}_{ \mathcal{D}}^{k}$, we denote by $\underline{v}_K = \big(v_K, (v_{ \sigma})_{ \sigma\in \mathcal{E}_K}\big)\in \underline{V}_K^{k}$ its local restriction to the cell $K$. To any $\underline{v}_{ \mathcal{D}} \in \underline{V}_{ \mathcal{D}}^{k}$, we associate the two piecewise polynomial functions $v_ \mathcal{M} : \Omega \to \mathbb{R}$ and $v_ \mathcal{E} : \partial \mathcal{M} \to \mathbb{R}$ such that

We also let $\underline{1}_ \mathcal{D} \in \underline{V}_ \mathcal{D} ^{k}$ be the discrete element such that $1_K = 1$ for all $K \in \mathcal{M}$ and $1_ \sigma = 1$ for all $\sigma \in \mathcal{E}$. Last, given a cell $K \in \mathcal{M}$, we define the local interpolator $\underline{I}^k_K : W^{1, 1}(K) \to \underline{V}_K^{k}$ such that, for any $v \in W^{1, 1}(K)$,

Similarly, the global interpolator $\underline{I}^k_ \mathcal{D} : W^{1, 1}(\Omega) \to \underline{V}_ \mathcal{D}^{k}$ is defined, for any $v \in W^{1, 1}(\Omega)$, by

As standard in the HHO context, locally to any cell $K\in \mathcal{M}$, we introduce a discrete gradient operator $G_K : \underline{V}_K^{k} \to \mathbb{P}^k(K)^d$ such that, for any $\underline{v}_K \in \underline{V}_K^{k}$, $G_K(\underline{v}_K)\in \mathbb{P}^k(K)^d$ satisfies

This operator is a consistent discrete counterpart of the gradient operator. It satisfies the following commutation property:

Given a face $\sigma \in \mathcal{E}_K$, we also define the jump operator $J_{K, \sigma} : \underline{V}_K^{k} \to \mathbb{P}^k(\sigma)$ such that, for $\underline{v}_K\in \underline{V}_K^k$,

Based on the above ingredients, one can define an HHO counterpart of the local diffusion bilinear form $(z, v) \mapsto \int_K \Lambda \nabla z \cdot \nabla v$. We let $a_K : \underline{V}_K^{k} \times \underline{V}_K^{k} \to \mathbb{R}$ be the bilinear form such that

where $\Lambda_{K, \sigma} = \|\Lambda_{\mid K} n_{K, \sigma} \cdot n_{K, \sigma} \|_{L^\infty(\sigma)}$ (recall the regularity assumption (2.1)). In the context of HDG methods, the linear stabilisation used in (2.4) is often called Lehrenfeld–Schöberl stabilisation, as it was first introduced in ^[39,40]. Classically, one can then define a global bilinear form $a_ \mathcal{D} : \underline{V}_ \mathcal{D}^{k} \times \underline{V}_ \mathcal{D}^{k} \to \mathbb{R}$, discretisation of $(z, v) \mapsto \int_\Omega \Lambda \nabla z \cdot \nabla v$, by summing the local contributions:

For analysis purposes, we introduce a discrete $H^1$-like semi-norm on $\underline{V}^k_ \mathcal{D}$. Given a cell $K \in \mathcal{M}$, we first let, for any $\underline{v}_K\in \underline{V}^k_K$,

Then, at the global level, for any $\underline{v}_ \mathcal{D}\in \underline{V}^k_ \mathcal{D}$, we define

Notice that $|{\cdot}|_{1, \mathcal{D}}$ is not a norm on $\underline{V}_ \mathcal{D}^{k}$, but any $\underline{v}_ \mathcal{D} \in \underline{V}_ \mathcal{D}^{k}$ satisfying $| \underline{v}_ \mathcal{D}|_{1, \mathcal{D}} = 0$ is proportional to $\underline{1}_ \mathcal{D}$. In particular, this implies that $|{\cdot}|_{1, \mathcal{D}}$ is a norm on the zero-mass subspace of $\underline{V}_ \mathcal{D}^{k}$ defined by

Standard HHO analysis implies the following stability estimate:

where the multiplicative constant is proportional to $\lambda_\flat$. In particular, since $|{\cdot}|_{1, \mathcal{D}}$ is a norm on $\underline{V}_{ \mathcal{D}, 0}^{k}$, this estimate implies that $a_ \mathcal{D}$ is coercive on $\underline{V}_{ \mathcal{D}, 0}^{k}$. Finally, we recall the following discrete Poincaré–Wirtinger inequality (cf. [23, Theorem 6.5, $p = q = 2$] in the equal-order case):

2.3.   Exponential fitting scheme

The construction extends the ideas from ^[16]. Our (linear) scheme hinges on the exponential fitting strategy ^[10]. In a nutshell, the exponential fitting approach is based on the following rewriting of the PDE flux: letting $\omega = {\rm{e}}^{-\phi}$, and introducing the Slotboom variable $\rho = u/\omega$, one has

which allows to transform an advection-diffusion equation in $u$ into a purely diffusive, unconditionally coercive (by regularity of $\phi$, $\omega$ is a.e. uniformly bounded away from zero) problem in $\rho$. At the discrete level, the problem is solved in the Slotboom variable, which is sought in $\underline{V}_ \mathcal{D}^{k}$. The discrete density is then defined mimicking the relation $u = \omega \rho$.

In view of (2.9), in order to define our exponential fitting (mixed-order) HHO scheme, we need to introduce a discrete counterpart of the bilinear form $(\rho, v) \mapsto \int_\Omega \omega \Lambda \nabla \rho \cdot \nabla v$. To do so, given $K \in \mathcal{M}$, and leveraging the definition (2.4) of $a_K$, we let $a_K^\omega : \underline{V}_K^{k} \times \underline{V}_K^{k} \to \mathbb{R}$ be such that

where $\Lambda_{ K, \sigma}^\omega = \|\omega \Lambda_{\mid K} n_{K, \sigma} \cdot n_{K, \sigma} \|_{L^\infty(\sigma)}$. At the global level, as previously, we construct the bilinear form $a_ \mathcal{D}^\omega : \underline{V}_ \mathcal{D}^{k} \times \underline{V}_ \mathcal{D}^{k} \to \mathbb{R}$ by summing the local contributions:

We can now introduce the exponential fitting HHO scheme for Problem (1.1): find

such that, for all $n \in \mathbb{N}$,

For any solution $(\underline{\rho}_ \mathcal{D}^n)_{n \geq 1}$ to (2.12), we define a sequence of corresponding densities $(\underline{\mathfrak{u}}_ \mathcal{D}^{\omega, n})_{n \geq 1}$ as follows. To the discrete Slotboom variable $\underline{\rho}_ \mathcal{D} \in \underline{V}_ \mathcal{D}^{k}$ we associate the discrete density

defined, consistently with (2.12b), as the collection of (a priori non-polynomial) functions

The non-polynomial nature of the components of $\underline{\mathfrak{u}}_ \mathcal{D}^\omega$ is, here and in what follows, emphasised by the use of Gothic fonts. Finally, to any discrete density $\underline{\mathfrak{u}}^\omega_ \mathcal{D}$, we associate the two piecewise smooth functions $\mathfrak{u}^\omega_ \mathcal{M} : \Omega \to \mathbb{R}$ and $\mathfrak{u}^\omega_ \mathcal{E} : \partial \mathcal{M} \to \mathbb{R}$ such that

Remark that $\mathfrak{u}^\omega_ \mathcal{M} = \omega\rho_ \mathcal{M}$ in $\Omega$ and $\mathfrak{u}^\omega_ \mathcal{E} = \omega_{\mid\partial \mathcal{M}}\rho_ \mathcal{E}$ on $\partial \mathcal{M}$.

Remark 2 (Non-polynomial definition). Here, we choose to define a discrete density with (a priori) non-polynomial components. One could also think of defining a density with polynomial components, by multiplying (component by component) $\underline{\rho}_ \mathcal{D}$ by $\underline{I}^k_ \mathcal{D}(\omega)$. This is how the solution to the low-order HFV exponential fitting scheme of ^[16] was defined (in that case, both cell/face unknowns were constants).

Remark 3 (Initial condition). Notice that we do not define $\mathfrak{u}^{\omega, 0}_ \mathcal{M}$ in the same way as $\mathfrak{u}^{\omega, n+1}_{ \mathcal{M} }$. Indeed, we directly use in the definition (2.12c) the initial datum $u^{0}$. This choice is motivated by the following observation. Define $\rho^0 = u^0/\omega$, and let $\mathfrak{u}^{\omega, 0}_K = \omega_{\mid K}\Pi^{k+1}_K(\rho^0_{\mid K})$ for all $K\in \mathcal{M}$ in place of (2.12c). Then, the resulting discrete solution's mass is $M^\omega = \sum_{K \in \mathcal{M}} \int_K \mathfrak{u}^{\omega, 0}_K \neq M$ in general.

Recall that $\int_{\Omega}u^0 = M > 0$. Let $\underline{\rho}_ \mathcal{D}^\infty \in \underline{V}_ \mathcal{D}^{k}$ be defined as

One can easily check that $\underline{\rho}_ \mathcal{D}^\infty$ is the only steady solution to the exponential fitting scheme (2.12). Based on (2.13), $\underline{\rho}_ \mathcal{D}^\infty$ is associated to the discrete equilibrium density $\underline{\mathfrak{u}}_ \mathcal{D}^{\omega, \infty}$ such that

It follows that the reconstructed discrete equilibrium density $\mathfrak{u}^{\omega, \infty}_ \mathcal{M}$ (always) coincides with the thermal equilibrium (1.4) in $\Omega$. Such a striking property is, however, to be tempered by Remark 4 below. Following Remark 2, notice that if we had adopted instead a polynomial definition for the discrete densities, we would have obtained that $\underline{\mathfrak{u}}_ \mathcal{D}^{\omega, \infty} = \underline{I}^k_ \mathcal{D}(u^\infty)\in \underline{V}^k_ \mathcal{D}$, as was the case for the low-order exponential fitting HFV method of ^[16] (in that case, both cell/face unknowns were constants). A drawback of such a definition, compared to (2.13), is that the components of $\underline{I}^k_ \mathcal{D}(u^\infty)$ are not necessarily positive functions (note that this issue does not exist in the low- and equal-order HFV case).

Remark 4 (Alternative scheme definition). Let $\underline{\phi}_ \mathcal{D} = \underline{I}_{ \mathcal{D}}^k(\phi)\in \underline{V}_ \mathcal{D}^k$. Another definition of the exponential fitting HHO scheme consists in replacing $\omega$ by ${\rm{e}}^{-\phi_ \mathcal{M}}$ in the expressions of both $a_K^{\omega}$ (see (2.10)) and $\mathfrak{u}_K^{\omega, n+1}$ (see (2.12b)). Then, in place of (2.13), the following definition of discrete densities is adopted:

Such a scheme is somewhat closer to what one would encounter in the context of semiconductor models, since $\phi$ would be unknown, and sought, at the discrete level, in $\underline{V}_ \mathcal{D}^k$. In this case, the discrete equilibrium density $\underline{\mathfrak{u}}_ \mathcal{D}^{\omega, \infty}$ would satisfy, in place of (2.14), the same kind of relations as (2.22) and (2.23) below. Here, we rather choose to exploit the full knowledge we have of the potential $\phi$ to define the scheme.

2.4.   Nonlinear scheme

The construction extends the ideas from ^[16,45]. Our nonlinear scheme relies on a nonlinear reformulation of Problem (1.1) ^[13]. To do so, we introduce the logarithm potential $\ell = \log(u)$ and the quasi-Fermi potential $w = \ell + \phi$. At least formally, if $u$ is positive, one has the following relation on the PDE flux:

We choose to discretise the potentials as piecewise polynomials, i.e., we approximate $\ell$ and $w$ as discrete unknowns in $\underline{V}_ \mathcal{D}^{k}$. Then, mimicking the relation $u = {\rm{e}}^\ell$, each discrete density component is defined as the exponential of a polynomial, thus ensuring its positivity.

In view of (2.16), in order to define our nonlinear HHO scheme, we shall introduce a discrete counterpart of the map $(\ell; w, v)\mapsto \int_\Omega {\rm{e}}^\ell \Lambda \nabla w \cdot \nabla v$. Locally to any cell $K\in \mathcal{M}$, this discrete counterpart is built as the sum of a consistent (2.17a) and a stabilising (2.17b) contributions: for all $\underline{\ell}_K, \underline{w}_K, \underline{v}_K\in \underline{V}_K^k$, we let

We then introduce the local map $\mathcal{T}_K : \underline{V}_K^{k} \times \underline{V}_K^{k} \times \underline{V}_K^{k} \to \mathbb{R}$ such that

where $\varepsilon$ is a non-negative parameter and $a_K$ is the bilinear form defined by (2.4). At the global level, we finally define the map $\mathcal{T}_ \mathcal{D} : \underline{V}_ \mathcal{D}^{k} \times \underline{V}_ \mathcal{D}^{k} \times \underline{V}_ \mathcal{D}^{k} \to \mathbb{R}$ by summing the local contributions:

Remark 5 (Parameter $\varepsilon$). The map $\mathcal{T}_ \mathcal{D}$ is to be understood as a discretisation of $(\ell; w, v) \mapsto \int_\Omega ({\rm{e}}^\ell + \epsilon)\, \Lambda\nabla w \cdot \nabla v$, with $\epsilon$ of magnitude $\varepsilon h_ \mathcal{D}^{k+2}$. At the theoretical level, this $\epsilon$-perturbation of the model is necessary, at the moment, to show the existence result of Theorem 3.1. From a more practical viewpoint, the sensitivity of the method with respect to $\epsilon$ is not completely understood yet. First numerical experiments tend to show that the choices $\varepsilon = 1$ and $\varepsilon = 0$ produce essentially similar results. Concerning the choice of scaling factor $h_K^{k+2}$ in (2.18), it seems to yield in practice (when $\varepsilon > 0$) the expected orders of convergence. The influence of the $\epsilon$-term will be further investigated in future works.

Let $\underline{\phi}_ \mathcal{D} = \underline{I}_{ \mathcal{D}}^k(\phi)\in \underline{V}_ \mathcal{D}^k$. We can now introduce our nonlinear HHO scheme for Problem (1.1): find

such that, for all $n \in \mathbb{N}$,

For any solution $\left (\underline{\ell}_ \mathcal{D}^n \right)_{n \geq 1}$ to (2.20), we define a sequence of corresponding positive densities $\left (\underline{\mathfrak{u}}_ \mathcal{D}^n \right) _{n \geq 1}$ as follows. To the discrete logarithm potential $\underline{\ell}_ \mathcal{D} \in \underline{V}_ \mathcal{D}^{k}$ we associate the discrete density

defined, consistently with (2.20b), as the collection of positive (non-polynomial) functions

The non-polynomial nature of the components of $\underline{\mathfrak{u}}_ \mathcal{D}$ is, here also, emphasised by the use of Gothic fonts. Finally, to any discrete density $\underline{\mathfrak{u}}_ \mathcal{D}$, we associate the two (positive) piecewise smooth functions $\mathfrak{u}_ \mathcal{M} : \Omega \to \mathbb{R}^\star$ and $\mathfrak{u}_ \mathcal{E} : \partial \mathcal{M} \to \mathbb{R}^\star$ such that

Remark that $\mathfrak{u}_ \mathcal{M} = {\rm{e}}^{\ell_ \mathcal{M}}$ in $\Omega$ and $\mathfrak{u}_ \mathcal{E} = {\rm{e}}^{\ell_ \mathcal{E}}$ on $\partial \mathcal{M}$.

Remark 6 (Initial condition). Remark that we do not define $\mathfrak{u}_ \mathcal{M}^0$ in the same way as $\mathfrak{u}_ \mathcal{M}^{n+1}$. We indeed directly use in the definition (2.20c) the initial datum $u^0$. This strategy allows one to circumvent the definition of some $\ell_ \mathcal{M}^0$, cell interpolate of $\ell^{0} = \log(u^{0})$, the latter quantity being undefined in regions where $u^{0}$ vanishes. The question of defining an initial discrete logarithm potential remains however a major difficulty when it comes to numerical implementation, since it is needed for the initialisation of the Newton method when $n = 0$ (see Section 4.1.2).

Recall that $\int_{\Omega}u^0 = M > 0$. Let $\underline{\ell}_ \mathcal{D}^\infty\in \underline{V}^k_ \mathcal{D}$ be defined as

It can be easily checked that $\underline{\ell}_ \mathcal{D}^\infty$ is the only steady solution to the nonlinear scheme (2.20). Based on (2.21), $\underline{\ell}_ \mathcal{D}^\infty$ is associated to the discrete equilibrium density $\underline{\mathfrak{u}}_ \mathcal{D}^\infty$ such that

In a sense we make clear just below, the discrete equilibrium density $\underline{\mathfrak{u}}_ \mathcal{D}^\infty$ is a reasonable approximation of the thermal equilibrium $u^\infty$ defined by (1.4):

It follows that the reconstructed discrete equilibrium density $\mathfrak{u}^{\infty}_ \mathcal{M}$ satisfies: if $\phi_{\mid K}\in \mathbb{P}^{k+1}(K)$ for all $K\in \mathcal{M}$, then $\mathfrak{u}_ \mathcal{M}^\infty = u^\infty$ in $\Omega$. Remark that the discrete equilibrium density $\underline{\mathfrak{u}}_ \mathcal{D}^\infty$ is not equal to $\underline{I}_{ \mathcal{D}}^k(u^\infty)\in \underline{V}_ \mathcal{D}^k$. This is in contrast with what held true for the low-order nonlinear HFV method of ^[16] (in that case, both cell/face unknowns were constants). This can be explained by the choice of discretisation for $\phi$, which was taken as $-\log(\underline{\omega}_ \mathcal{D})$ (with $\underline{\omega}_ \mathcal{D}$ HFV interpolate of $\omega$) in ^[16] in place of $\underline{\phi}_ \mathcal{D}$ here (the latter choice is inspired by ^[45] in the context of semiconductors). Remark that, in practice, $-\log(\underline{\omega}_ \mathcal{D})$ and $\underline{\phi}_ \mathcal{D}$ may coincide if the integrals are approximated using an evaluation at the barycenter. This was the case in the numerical experiments of ^[16].

3.   Main features of the schemes

We present in this section some theoretical results about the two schemes introduced above. We focus, in particular, on the existence (and stability) of solutions, as well as on questions related to their long-time behaviour. The results presented below generalise those obtained in ^[16] in the low-order HFV context. In particular, the analysis strongly hinges on the entropy structure of both schemes.

Remark 7 (Lowest-order versions of the schemes ($k = 0$)). Note that the lowest-order versions of the two schemes introduced above do not coincide with the exponential fitting and nonlinear HFV schemes of ^[16]. Indeed, the lowest-order versions of the methods (2.12) and (2.20) make use of (enriched) affine cell unknowns, whereas HFV schemes use constants. Also, whereas the nonlinear HFV method is built upon a stable discrete gradient operator (defined on a pyramidal submesh), the present nonlinear scheme is defined following the standard HHO philosophy of splitting consistency and (nonlinear) stabilisation. Therefore, the results presented here are new, even for $k = 0$.

Before presenting individual results for each scheme, let us stress that both schemes exhibit a similar important property: the preservation of the thermal equilibrium.

Proposition 1 (Preservation of the thermal equilibrium). The alternative (fully discrete) exponential fitting scheme of Remark 4 and the nonlinear scheme (2.20) preserve the thermal equilibrium in the sense of Definition 1. More precisely,

● $\log(\underline{\mathfrak{u}}^{\omega, \infty}_ \mathcal{D})+ \underline{\phi}_ \mathcal{D}$ is proportional to $\underline{1}_ \mathcal{D}$;

● $\underline{\ell}^\infty_ \mathcal{D} + \underline{\phi}_ \mathcal{D}$ is proportional to $\underline{1}_ \mathcal{D}$.

For the original exponential fitting scheme (2.12), there holds that $\log(\mathfrak{u}^{\omega, \infty}_ \mathcal{M}) + \phi$ over $\Omega$ and $\log(\mathfrak{u}^{\omega, \infty}_ \mathcal{E}) + \phi_{\mid\partial \mathcal{M}}$ over $\partial \mathcal{M}$ are equal to the same constant.

3.1.   Exponential fitting scheme

We present here the main properties of the exponential fitting HHO scheme (2.12), and give detailed proofs of the results. As a preliminary remark, note that since $\phi\in W^{1, \infty}(\Omega)$,

As a consequence, the tensor $\omega \Lambda$ is a.e. uniformly elliptic. Recalling (2.7), this implies the following stability estimate:

where the multiplicative constant is proportional to $\lambda_\flat$. We first state a well-posedness result, which is mainly a consequence of the previous stability estimate.

Proposition 2 (Well-posedness of the exponential fitting scheme). The exponential fitting scheme (2.12) admits a unique solution $(\underline{\rho}_ \mathcal{D}^n)_{n\geq 1}$. Moreover, the corresponding discrete densities $(\underline{\mathfrak{u}}_ \mathcal{D}^{\omega, n})_{n \geq 1}$ have a mass equal to $M$:

Proof. Let $n \geq 0$, and assume that $\mathfrak{u}_ \mathcal{M}^{\omega, n}$ is defined. We want to show that the system (2.12a)- (2.12b) admits a unique solution. To do so, we first define, for any $\underline{\rho}_ \mathcal{D} \in \underline{V}_ \mathcal{D}^{k}$,

Since $|{\cdot}|_{1, \mathcal{D}}$ is a semi-norm on $\underline{V}_ \mathcal{D}^{k}$ with zero set spanned by $\underline{1}_ \mathcal{D}$, it follows that the map $\|{\cdot}\| _{\omega, \Delta t, \mathcal{D}}$ defines a norm on $\underline{V}_ \mathcal{D}^{k}$. Thus, by (3.2), the bilinear form $A_ \mathcal{D}^\omega: (\underline{\rho}_ \mathcal{D}, \underline{v}_ \mathcal{D}) \mapsto \int_\Omega\omega \rho_ \mathcal{M} v_ \mathcal{M} + \Delta t a_ \mathcal{D}^\omega(\underline{\rho}_ \mathcal{D}, \underline{v}_ \mathcal{D})$ satisfies the following coercivity property:

By the Lax-Milgram lemma, the system (2.12a)-(2.12b) therefore admits a unique solution $\underline{\rho}_ \mathcal{D}^{n+1}$ in $\underline{V}_ \mathcal{D}^{k}$, from which one can uniquely define $\underline{\mathfrak{u}}^{\omega, n+1}_ \mathcal{D}$ by (2.13). To prove mass conservation, we just test (2.12a) by $\underline{1}_ \mathcal{D}$ to get

We conclude by noticing that $\int_\Omega\mathfrak{u}^{\omega, 0}_ \mathcal{M} = \int_\Omega u^{0}$ according to (2.12c). □

We now state our main result about the exponential fitting scheme, which ensures that the solution to (2.12) has similar long-time behaviour as the PDE solution. As usual with the entropy method, the main idea is to get a control of the entropy by its dissipation. Here, such an estimate is a consequence of the discrete Poincaré inequality (2.8).

Proposition 3 (Long-time behaviour of the exponential fitting scheme). Assume that $u^0\in L^2(\Omega)$. Let $(\underline{\rho}_ \mathcal{D}^n)_{n \geq 1}$ be the solution to the exponential fitting scheme (2.12). Then, the following discrete entropy relation holds true:

where the discrete quadratic entropy is defined as

with $\rho_ \mathcal{M}^0$ defined (with a slight abuse in notation, since $\rho_ \mathcal{M}^0$ is not piecewise polynomial) by $\rho_ \mathcal{M}^0 = \frac{u^{0}}{\omega}$, and the discrete dissipation is given by

As a consequence, the reconstructed discrete density converges exponentially fast in time towards the reconstructed discrete equilibrium density: there exists a positive constant $\nu_\omega$, independent of both $h_ \mathcal{D}$ and $\Delta t$, such that

Proof. Let $n \in \mathbb{N}$. By convexity of $x\mapsto x^2$ on $\mathbb{R}$, one has

Therefore, testing (2.12a) against $\underline{\rho}_ \mathcal{D}^{n+1} - \underline{\rho}_ \mathcal{D}^\infty \in \underline{V}_ \mathcal{D}^{k}$, we get

Note that this estimate holds true also for $n = 0$ (using the definition of $\rho_ \mathcal{M}^{0}$). On the other hand, by the expression of $\underline{\rho}_ \mathcal{D}^\infty$ (proportional to $\underline{1}_ \mathcal{D}$), $a_ \mathcal{D}^\omega (\underline{\rho}_ \mathcal{D}^{\infty}, \underline{\rho}_ \mathcal{D}^{n+1} - \underline{\rho}_ \mathcal{D}^\infty) = 0$, hence by bilinearity of $a_ \mathcal{D}^\omega$,

which yields the entropy relation (3.4). To get the exponential decay, one needs to compare $\mathbb{D}_\omega^{n+1}$ with $\mathbb{E}^{n+1}_\omega$. To do so, we let $\underline{v}_ \mathcal{D} = \underline{\rho}_ \mathcal{D}^{n+1} - \underline{\rho}_ \mathcal{D}^\infty \in \underline{V}_ \mathcal{D}^{k}$, and we define the probability measure $\mathrm{d}\mu = \frac{c^M_l}{M}\omega \mathrm{d} x$ on $\Omega$. We define $\langle v_ \mathcal{M}\rangle_\mu$ as the mass of $v_ \mathcal{M}$ for the measure $\mathrm{d}\mu$, i.e.,

The definition (2.13), and the mass preservation identity (3.3), imply that

Therefore, letting $\langle v_ \mathcal{M}\rangle = \frac{1}{|\Omega|_d} \int_\Omega v_ \mathcal{M}$, and applying [11, Lemma 5.2, $q = 2$], we get

Using the definition of $\mathrm{d}\mu$, and the bound (3.1), yields

By definition of $\langle v_ \mathcal{M}\rangle$, one has $\underline{v}_ \mathcal{D} - \langle v_ \mathcal{M}\rangle \underline{1}_ \mathcal{D} \in \underline{V}_{ \mathcal{D}, 0}^{k}$, so we can apply the discrete Poincaré–Wirtinger inequality (2.8) to infer that

Combining the two previous estimates, we get

Now, one can use the stability estimate (3.2) to infer that

Therefore, combining the last two estimates, one infers the existence of $\nu_\omega > 0$, independent of both $h_ \mathcal{D}$ and $\Delta t$, such that the following relation between the entropy and its dissipation holds true:

Plugging this estimate into the entropy relation (3.4), we deduce that

This implies the exponential decay of the entropy:

To conclude, we just use (2.13) and (3.1) to infer that

which, combined with the fact that $\mathfrak{u}^{\omega, 0}_ \mathcal{M} = u^0$, finally yields (3.5). □

Remark 8. (Regularity of the initial datum and topology of the convergence). Notice that in Proposition 3 we have made the extra assumption that $u^{0} \in L^2(\Omega)$. The long-time analysis of the exponentially fitted model indeed relies on the decay of the quadratic entropy (in the Slotboom variable)

In order to guarantee that the initial quadratic entropy is finite, assuming that the initial datum is in $L^2(\Omega)$ is a safe choice. At the end, as a reminiscence of the linearity of the model, the exponential fitting approach gives a convergence (in time) result in the $L^2$-topology (in space). In contrast, the nonlinear approach will yield convergence in a weaker norm (typically $L^1$), but can be used to deal with less regular initial data, which are in $L \log(L)$ only.

3.2.   Nonlinear scheme

We present here some results regarding the analysis of the nonlinear HHO method (2.20). Since we deal with a nonlinear scheme, unlike the exponential fitting scheme, the question of the existence of solutions is the main difficulty here. As often for this type of method, we start by establishing some a priori estimates. For the purpose of analysis, given a discrete logarithm potential $\underline{\ell}_ \mathcal{D} \in \underline{V}_ \mathcal{D}^{k}$, we associate a discrete quasi-Fermi potential $\underline{w}_ \mathcal{D}\in \underline{V}^k_ \mathcal{D}$ defined by

where we recall that $c^M_{nl} = M/\int_\Omega{\rm{e}}^{-\phi_ \mathcal{M}}$. By (2.21) and (2.22), one has

Note that, on the other hand, for any $\underline{\ell}_ \mathcal{D}, \underline{v}_ \mathcal{D}\in \underline{V}_ \mathcal{D}^{k}$, we have

since $(\underline{\ell}_ \mathcal{D} + \underline{\phi}_ \mathcal{D})- \underline{w}_ \mathcal{D}$ is proportional to $\underline{1}_ \mathcal{D}$. Similarly to previous works on nonlinear HFV schemes for semiconductor models ^[45], the discrete quasi-Fermi potentials are the key variables to perform the analysis of the method. As a last remark, notice that since $\phi\in H^1(\Omega)$, by boundedness of the interpolator (cf. [23, Proposition 5.3]),

Let us now state some fundamental a priori relations. As for the exponential fitting scheme, the discrete entropy structure of the nonlinear scheme mainly results from the convexity of the entropy.

Proposition 4 (Fundamental a priori relations). Let $\left (\underline{\ell}_ \mathcal{D}^n \right)_{n \geq 1}$ be a given solution to the nonlinear scheme (2.20), and $\left (\underline{\mathfrak{u}}_ \mathcal{D}^n \right) _{n \geq 1}$ be the corresponding discrete density. Then, the following a priori relations hold true:

(i) the mass is preserved along time:

(ii) a discrete entropy/dissipation relation is satisfied:

where the discrete entropy and dissipation are non-negative quantities defined by

with $\Phi_1:s \mapsto s \log(s) - s +1$ (and $\Phi_1(0) = 1$).

Proof. Let $n\geq 0$. Using $\underline{1}_ \mathcal{D}$ as a test function in (2.20a), we get that the mass is conserved:

Therefore, by (2.20c), we infer (3.9). To establish the entropy relation, we first use the convexity of $\Phi_1$, which yields

Then, since $w_ \mathcal{M}^{n+1} = \log \left (\frac{\mathfrak{u}_ \mathcal{M}^{n+1}}{\mathfrak{u}_ \mathcal{M}^\infty} \right)$ and $\Phi_1 ' = \log$, one has

On the other hand, testing (2.20a) with $\underline{w}_ \mathcal{D}^{n+1} \in \underline{V}_ \mathcal{D}^{k}$, and using (3.7), we get

which finally yields (3.10) by definition of the discrete dissipation. □

Remark that since $\mathfrak{u}_ \mathcal{M}^0 = u^0$, and $u^0\geq 0$ in $\Omega$, $u^0\in L^1(\Omega)$, and $\int_\Omega u^{0} \log(u^{0}) < \infty$, one has $\mathbb{E}^0 < \infty$. Note finally that the previous results hold true for any $\varepsilon \geq 0$ in (2.18).

In the rest of this section, we focus on the existence of solutions and on their long-time behaviour. We henceforth assume that $\varepsilon > 0$. The proofs for both results rely on a discrete a priori estimate, which is obtained by means of a high-order counterpart of [16, Lemma 2]. In order to perform the analysis, we first introduce an inner product $\langle \cdot, \cdot \rangle$ on $\underline{V}_ \mathcal{D}^{k}$:

We denote by $\|{\cdot}\|$ the corresponding Euclidean norm:

Lemma 1 (Discrete boundedness by mass and energy semi-norm). Let $\underline{\ell}_ \mathcal{D} \in \underline{V}_ \mathcal{D}^{k}$, and assume that there exist $C_\sharp > 0$ and $M_\sharp \geq M_\flat > 0$ such that

Then, there exists a positive constant $C$, only depending on $M_\flat$, $M_\sharp$, $C_\sharp$, $\Omega$, $\theta$, $k$ and $h_\flat$ such that

Proof. Let us first remark that

Hence, to estimate $\| \underline{\ell}_ \mathcal{D}\|$, all that remains to bound is $\|\ell_ \mathcal{M}\|_{L^2(\Omega)}$. Recalling the notation $\langle z\rangle = \frac{1}{|\Omega|_d}\int_{\Omega}z$, and applying the discrete Poincaré–Wirtinger inequality (2.8), it holds

with $C_{PW} > 0$ only depending on $\Omega$, $\theta$ and $k$. Thus, by the triangle inequality, we infer

and we are only left with estimating $|\langle\ell_ \mathcal{M}\rangle|$. We proceed in two steps, showing first an upper bound on $\langle\ell_ \mathcal{M}\rangle$, and then a lower bound. Applying Jensen's inequality, the upper bound can be readily obtained:

which yields $\langle\ell_ \mathcal{M}\rangle\leq\log\big(\frac{M_\sharp}{|\Omega|_d}\big)$. To prove the lower bound, we start from (3.13), and we use local reverse Lebesgue embedding (cf. [23, Lemmas 1.25 and 1.12]). This yields

where $C_{RL} > 0$ only depends on $d$, $\theta$ and $k$. Then, remarking that

and integrating over $\Omega$, we get

Now, using the lower bound on $\int_{\Omega}{\rm{e}}^{\ell_ \mathcal{M}}$, and taking the logarithm, we finally infer that

This concludes the proof. □

We now state the existence result, which holds true for positive $\varepsilon$. The proof adopts the methodology developed in ^[16] in the (nonlinear) HFV context.

Theorem 1 (Existence of solutions to the nonlinear scheme (2.20)). Assume that the stabilisation parameter $\varepsilon$ in (2.18) is positive. Then, there exists at least one solution $\left (\underline{\ell}_ \mathcal{D}^n \right)_{n \geq 1}$ to the scheme (2.20). The corresponding discrete densities $\left (\underline{\mathfrak{u}}_ \mathcal{D}^n \right) _{n \geq 1}$, defined by (2.21), are positive.

Proof. The proof proceeds by induction. Let $n \in \mathbb{N}$, and assume that $\mathfrak{u}_ \mathcal{M}^n$ is well defined, following (2.20b) (if $n \geq 1$) or (2.20c) (if $n = 0$). We now prove the existence of a solution $\underline{\ell}_ \mathcal{D}^{n+1} \in \underline{V}_ \mathcal{D}^{k}$ to (2.20a). For convenience, instead of looking for the discrete logarithm potential, we will equivalently seek for the discrete quasi-Fermi potential $\underline{w}_ \mathcal{D}^{n+1} = \underline{\ell}_ \mathcal{D}^{n+1} + \underline{\phi}_ \mathcal{D} - \log(c^M_{nl}) \underline{1}_ \mathcal{D}$ (cf. (3.6)).

First, notice that, given any $\underline{w}_ \mathcal{D} \in \underline{V}_ \mathcal{D}^{k}$, and corresponding discrete logarithm potential $\underline{\ell}_ \mathcal{D}$ (through (3.6)) and discrete density $\underline{\mathfrak{u}}_ \mathcal{D}$ (through (2.21)), the map

is a bounded linear form on $\underline{V}_ \mathcal{D}^{k}$. Therefore, by the Riesz–Fréchet representation theorem, there exists a unique element $\underline{\mathcal{G}}_ \mathcal{D}(\underline{w}_ \mathcal{D}) \in \underline{V}_ \mathcal{D}^{k}$ such that

Remark that $\underline{w}_ \mathcal{D} \mapsto \underline{\mathcal{G}}_ \mathcal{D}(\underline{w}_ \mathcal{D})$ is a continuous (nonlinear) map of $\underline{V}_ \mathcal{D}^{k}$. Note also that, for any discrete quasi-Fermi potential $\underline{w}_ \mathcal{D} \in \underline{V}_ \mathcal{D}^{k}$ such that $\underline{\mathcal{G}}_ \mathcal{D}(\underline{w}_ \mathcal{D}) = \underline{0}_ \mathcal{D}$, by (3.7), the corresponding discrete logarithm potential $\underline{\ell}_ \mathcal{D}$ solves (2.20a). Our aim from now on is thus to show that $\underline{\mathcal{G}}_ \mathcal{D}$ does vanish on $\underline{V}_ \mathcal{D}^{k}$.

To this purpose, we introduce a regularisation of $\underline{\mathcal{G}}_ \mathcal{D}$: given any $\mu > 0$, we let

By definition of $\underline{\mathcal{G}}_ \mathcal{D}$, one has

As already shown in the proof of Proposition 4 (cf. (3.11)), by convexity of $\Phi_1$, one has

where the discrete entropies are defined by

As already mentioned, since $\Phi_1$ is a non-negative function, these two quantities are non-negative. Note that it may occur that $\mathbb{E}^n = 0$ (which is equivalent to $\mathfrak{u}^n_ \mathcal{M} = \mathfrak{u}_ \mathcal{M}^\infty$ in $\Omega$ for $n\geq 1$, or $u^0 = \mathfrak{u}_ \mathcal{M}^\infty$ in $\Omega$), in which case $\underline{\ell}_ \mathcal{D} = \underline{\ell}_ \mathcal{D}^\infty$ is the unique solution to (2.20a) (uniqueness follows from the entropy relation (3.10)). In the following, we therefore assume that $\mathbb{E}^n > 0$. The previous identities, and the non-negativity of the dissipation and entropy, imply that

Letting $r = \sqrt{ \frac{ \mathbb{E}^n}{\mu \Delta t }} > 0$, one has that $\left \langle \underline{\mathcal{G}}_ \mathcal{D}^\mu(\underline{w}_ \mathcal{D}), \underline{w}_ \mathcal{D} \right \rangle \geq 0$ for all $\underline{w}_ \mathcal{D} \in \underline{V}_ \mathcal{D}^{k}$ such that $\| \underline{w}_ \mathcal{D} \| = r$. Therefore, according to [16, Lemma 1] (cf. also [32, Section 9.1]), which is a by-product of Brouwer's fixed-point theorem, there exists $\underline{w}_ \mathcal{D}^\mu \in \underline{V}_ \mathcal{D}^{k}$ such that

Now, plugging $\underline{w}_ \mathcal{D}^\mu$ in (3.14), and using that $\underline{\mathcal{G}}_ \mathcal{D}^\mu(\underline{w}_ \mathcal{D}^\mu) = \underline{0}_ \mathcal{D}$, we get

so that $\mathcal{T}_ \mathcal{D}(\underline{\ell}_ \mathcal{D}^\mu; \underline{w}_ \mathcal{D}^\mu, \underline{w}_ \mathcal{D}^\mu) \leq \frac{ \mathbb{E}^n }{\Delta t }$. Thus, recalling the definition (2.18)-(2.19) of $\mathcal{T}_ \mathcal{D}$, as well as the stability estimate (2.7) for $a_ \mathcal{D}$, we infer that

On the one hand, by (3.6) and the estimate (3.8) on $| \underline{\phi}_ \mathcal{D}|_{1, \mathcal{D}}$, it holds

On the other hand, by definition of $\underline{\mathcal{G}}_ \mathcal{D}^\mu$, one first infers that

Second, using the Cauchy-Schwarz inequality, followed by the bound (3.15) on $\| \underline{w}_ \mathcal{D}^\mu\|$, one gets

where we have also used that $\| \underline{1}_ \mathcal{D}\| = |\Omega|_d^{{}^{1}\!\!\diagup\!\!{}_{2}\;}$. Thus, letting $M^{n} = \int_\Omega \mathfrak{u}_ \mathcal{M}^n > 0$ (recall that $\int_{\Omega}\mathfrak{u}^0_ \mathcal{M} = M > 0$), and $\mu^n = \frac{(M^{n})^2}{4 \Delta t \mathbb{E}^n |\Omega|_d} > 0$, for all $0 < \mu \leq \mu^n$ one has

Leveraging (3.16) and (3.17), one can eventually apply Lemma 1 with $M_\flat = \frac{M^n}{2}$, $M_\sharp = \frac{3M^n}{2}$, and $C_\sharp$ proportional to $\sqrt{\frac{ \mathbb{E}^n }{\varepsilon h_\flat^{k+2} \Delta t }} +1$ (note that these three constants do not depend on $\mu$): there exists a constant $C > 0$, independent of $\mu$, such that

Then, by compactness, there exists $\underline{\ell}_ \mathcal{D}^{n+1} \in \underline{V}_ \mathcal{D}^{k}$ such that, up to extraction (not relabelled), $\underline{\ell}_ \mathcal{D}^\mu \to \underline{\ell}_ \mathcal{D}^{n+1}$ when $\mu \to 0$. On the other hand, $\underline{\mathcal{G}}_ \mathcal{D}^\mu$ tends to $\underline{\mathcal{G}}_ \mathcal{D}$ as $\mu$ tends to $0$. Therefore, letting

we have $\underline{0}_ \mathcal{D} = \underline{\mathcal{G}}_ \mathcal{D}^\mu(\underline{w}_ \mathcal{D}^\mu) \to \underline{\mathcal{G}}_ \mathcal{D}(\underline{w}_ \mathcal{D}^{n+1})$ as $\mu \to 0$, which implies that

It follows that $\underline{\ell}_ \mathcal{D}^{n+1}$ is a solution to (2.20a). □

Remark 9 (Uniqueness of the solution). As for the low-order nonlinear VAG, DDFV and HFV schemes of ^[12,13,16], the uniqueness of the solution to (2.20) is still an open question. A possible approach to show such a result could be to consider the relative discrete entropy of a solution with respect to another solution, and show that this quantity vanishes.

Last, we study the long-time behaviour of the nonlinear HHO scheme.

Proposition 5 (Long-time behaviour of the nonlinear scheme). Assume that the stabilisation parameter $\varepsilon$ in (2.18) is positive, and let $\left (\underline{\ell}_ \mathcal{D}^n \right)_{n \geq 1}$ be a solution to the scheme (2.20). Then, the discrete solution converges in time towards the discrete equilibrium logarithm potential:

Consequently, the corresponding reconstructed discrete density $(\mathfrak{u}^n_ \mathcal{M})_{n \geq 1}$ converges to $\mathfrak{u}_ \mathcal{M}^\infty$ in $L^\infty(\Omega)$.

Proof. First, remark that owing to the entropy relation (3.10), one has

Thus, according to the definition of the discrete dissipation $\mathbb{D}^n$, alongside with the definition (2.18)-(2.19) of $\mathcal{T}_ \mathcal{D}$, and the stability estimate (2.7) for $a_ \mathcal{D}$, we infer that

This implies, in particular, that

Let $n \geq 1$. By (3.6) and (3.8), one has

On the other hand, by the mass preservation (3.9), we have $\int_\Omega {\rm{e}}^{\ell^n_ \mathcal{M}} = M > 0$. Therefore, one can apply Lemma 1, and infer the existence of a positive constant $C$ (which is independent of $n$) such that

It follows, by compactness, that there exists $\underline{\ell}_ \mathcal{D} \in \underline{V}_ \mathcal{D}^{k}$ such that, up to extraction (not relabelled),

By (3.6), (3.19), and continuity of $|{\cdot}|_{1, \mathcal{D}}$ on $\underline{V}_ \mathcal{D}^{k}$, we infer that

This means that there exists $a \in \mathbb{R}$ such that $\underline{\ell}_ \mathcal{D} + \underline{\phi}_ \mathcal{D} = a \underline{1}_ \mathcal{D}$. By mass preservation, we get

so that $a = \log(c^M_{nl})$, which implies that

By uniqueness of the limit, we finally infer the convergence of the whole sequence $\left(\underline{\ell}_ \mathcal{D}^n \right)_{n \geq 1}$ towards $\underline{\ell}_ \mathcal{D}^\infty$ in $\underline{V}^k_ \mathcal{D}$. This implies, in particular, that $\ell^n_ \mathcal{M} \xrightarrow[n\to\infty]{} \ell_ \mathcal{M}^\infty$ in $L^\infty(\Omega)$, by norm equivalence in finite-dimensional vector spaces. Then, by the mean value theorem, we deduce that

which implies, by uniform boundedness (in $n$) of $\left(\ell_ \mathcal{M}^n \right)_{n \geq 1}$, the convergence of the reconstructed discrete density in $L^\infty(\Omega)$. □

Remark 10 (Non-uniformity of the bounds). Note that the estimate (3.20) on the solution to (2.20) is not uniform with respect to the discretisation parameters $h_ \mathcal{D}$ and $\Delta t$, nor with respect to the stabilisation parameter $\varepsilon$. Indeed, having a closer look to the dependencies of the corresponding upper bound (using Lemma 1), one realises that it blows up as soon as either $h_ \mathcal{D}$, $\Delta t$ or $\varepsilon$ tends to zero.

Remark 11 (Convergence to equilibrium). Notice that the time-asymptotic result of Proposition 5 is relatively weaker than the one of Proposition 3 in the exponential fitting context. For the latter result, the convergence to equilibrium is shown to be exponentially fast, and the decay rate $\nu_\omega$ uniform with respect to the discretisation parameters. The numerical results of ^[44] and Section 4.4 indicate that, also for the nonlinear scheme, the convergence is expected to be exponential, with seemingly uniform (and close to the PDE model one) decay rate. At the theoretical level, to prove exponential convergence to equilibrium, one has to establish a control of the discrete entropy by the discrete dissipation. Adapting the arguments from [45, Theorem 3], such a control can actually be established in the present context, but leads to a non-uniform decay rate (also depending on $\varepsilon$), and to a final result still only valid in the case $\varepsilon > 0$. In order to showcase a uniform (and $\varepsilon$-independent) decay rate, and establish a result also valid in the case $\varepsilon = 0$, a (high-order) uniform discrete Logarithmic-Sobolev inequality needs to be available (cf. ^[16] in the low-order HFV context). This is the subject of ongoing research.

4.   Numerical results

In this section, we extensively assess the high-order nonlinear scheme (2.20). We study positivity preservation, convergence, efficiency (accuracy vs. computational cost), and long-time behaviour. We also compare it, in terms of structure preservation, with the linear high-order exponential fitting scheme (2.12). All the test-cases considered below are set in the two-dimensional domain $\Omega = (0, 1)^2$, and are (except for the last one) taken from ^[16], to which we refer for more detailed descriptions. Given a (face) degree $k\geq 0$, the nonlinear scheme (2.20) will be referred to as $\texttt{nlhho_k}$, whereas the exponential fitting one (2.12) as $\texttt{expf_k}$. For the nonlinear scheme, we will always use below the value $\varepsilon = 1$ for the parameter $\varepsilon$ in (2.18). However, in some situations, we will compare the two values $\varepsilon = 1$ and $\varepsilon = 0$. The nonlinear scheme with $\varepsilon = 0$ will then be denoted $\texttt{nlhho_k_0}$.

4.1.   Implementation

All numerical tests presented below have been run on a laptop equipped with an Intel Core i7-9850H processor clocked at 2.60GHz and 32Gb of RAM. Our HHO implementation makes use of monomial basis functions for both the cell and face unknowns. Such a choice is known to introduce numerical instabilities for large values of $k$, we thus restrict our study to $k\leq 3$. The use of orthonormal basis functions, which is expected to improve on this situation (in particular for the convergence of the Newton algorithm in the nonlinear case), shall be studied in future works. We use quadrature formulas based on the Dunavant rules ^[31] (after subtessellation). To cope with non-polynomial integrands, we employ quadrature formulas of order $2k+5$. We performed a few tests (not reported here) with higher-order formulas, and did not observe any significant changes. Last, the local computations are performed sequentially. One could expect a significant gain in terms of performances parallelising the latter. We discuss below some important implementation aspects for both schemes.

4.1.1.   Exponential fitting scheme

For the linear exponential fitting scheme, the implementation follows the classical HHO strategy for linear diffusion problems. We directly solve for the discrete Slotboom variable $(\underline{\rho}_ \mathcal{D}^n)_{n \geq 1}$. As standard for skeletal methods, we do not solve the full linear system, but first perform static condensation, which allows one to locally eliminate the cell unknowns. Since the scheme relies on the same LHS matrix at each time step, we perform once and for all an LU decomposition of the matrix at the beginning of the computation. At each time step, the solution is then inexpensive (the RHS has to be updated, but only through a matrix-vector product).

We do not address in this work the main questions which were highlighted in [16, Section 5.1.2] in the low-order HFV context, about the (harmonic) averaging of $\omega$ (which is related to the choice of quadrature formulas for the high-order scheme), and the preconditioning of the system (which was equivalent, in the simple HFV context, to choose to solve the system in the density variable). These aspects shall be investigated in future works. Nonetheless, in view of the results obtained in ^[16] for the HFV exponential fitting scheme, we believe these potential improvements will have no effect on the positivity violation issues.

4.1.2.   Nonlinear scheme

The numerical scheme (2.20) requires to solve a nonlinear system of equations at each time step. For $n \in \mathbb{N}$, one wants to find $\underline{\ell}_ \mathcal{D}^{n+1} \in \underline{V}_ \mathcal{D}^{k}$ solution to (2.20): this scheme can be written as the equation

with $\underline{\mathcal{G}}_ \mathcal{D}^{n, \Delta t} : \underline{V}_ \mathcal{D}^{k} \to \underline{V}_ \mathcal{D}^{k}$ smooth (nonlinear) vector field. Numerically, to find a zero of $\underline{\mathcal{G}}_ \mathcal{D}^{n, \Delta t}$, we use a Newton method.

In practice, the use of a naive method without any adaptation proves not to be enough to compute a solution in general. In order to get a robust implementation, which can handle various data and meshes, one has to deploy a few techniques. For further use, we let $\| \underline{\ell}_ \mathcal{D}\|_{l^\infty}$ denote the $l^\infty$-norm of the coefficients of $\underline{\ell}_ \mathcal{D}$ in the (cell and face) polynomial bases. The map $\|{\cdot}\|_{l^\infty}$ is a norm on $\underline{V}_ \mathcal{D}^{k}$, which is easily (and at very low cost) computable in practice. To fix the ideas and the notation, the Newton method is defined as follows: given an initialisation $\underline{\ell}_{ \mathcal{D}, 0} \in \underline{V}_ \mathcal{D}^{k}$, and a time step $0 < \delta t \leq \Delta t$, one defines a sequence $(\underline{\ell}_{ \mathcal{D}, i})_{i\geq 0}$ of elements of $\underline{V}_ \mathcal{D}^{k}$ such that

where $\underline{\mathcal{G}}_ \mathcal{D}^{n, \delta t}$ is the vector field associated to the nonlinear scheme (2.20) with time step $\delta t$ instead of $\Delta t$, and $J^{n, \delta t}_{ \underline{\ell}_{ \mathcal{D}, i}}$ is the differential (Jacobian in practice) of $\underline{\mathcal{G}}_ \mathcal{D}^{n, \delta t}$ at $\underline{\ell}_{ \mathcal{D}, i}$. Note that, in practice, we do not solve this linear system, but perform static condensation in order to (locally) eliminate the cell unknowns. The resulting linear system is called "condensed system" in what follows. We discuss below the main tricks deployed to reach robustness in the implementation of the Newton algorithm.

Stopping criterion. We define the relative norm of the residual $r_{i+1}$, and the norm of the objective function $g_i$, as

We consider that the Newton method has converged when either

with $\texttt{tol} = 5. 10^{-10}$, in which case we set $\underline{\ell}_ \mathcal{D}^{n+1} = \underline{\ell}_{ \mathcal{D}, i+1}$. On the other hand, if this criterion is not met at $i = 50$, the method is considered as non-convergent (and we then proceed with a time step reduction, see below). In practice, for the tests collected in this article, we never reached $i = 50$, either because the method converged, or because of a loop break (see below).

Loop break for unreasonably large $\underline{\ell}_ \mathcal{D}$. The computations of $\underline{\mathcal{G}}_ \mathcal{D}^{n, \delta t} \left (\underline{\ell}_{ \mathcal{D}, i} \right)$ and $J^{n, \delta t}_{ \underline{\ell}_{ \mathcal{D}, i}}$ imply punctual evaluations of ${\rm{e}}^{\ell_{K, i}}$ (for $K\in \mathcal{M}$) and ${\rm{e}}^{\ell_{ \sigma, i}}$ (for $\sigma\in \mathcal{E}$) in the quadrature formulas. Such computations can lead to severe numerical issues if the values at the quadrature nodes are too large. Therefore, we declare that $\underline{\ell}_{ \mathcal{D}, i}$ is unreasonably large for the computations if there exists a cell quadrature node $x_{K, q} \in \overline{K}$, or a face quadrature node $x_{ \sigma, q} \in \overline{ \sigma}$, such that

In such a case, the method is immediately considered as non-convergent, and we proceed with a time step reduction (see below). Note that the choice of the value $100$ allows one to compute densities $\underline{\mathfrak{u}}_ \mathcal{D}$ over a range from $10^{-43}$ to $10^{43}$, and hence should not be a significant restriction in practice. In the numerical simulations presented below, the use of this loop-breaking procedure is absolutely necessary in order to avoid the evaluation of too large quantities, leading to some "explosion" of the method and crash of the code. Moreover, we also operate a loop break if the linear solver does not perform a successful resolution of the condensed linear system associated to (4.1), which corresponds to situations for which either $J^{n, \delta t}_{ \underline{\ell}_{ \mathcal{D}, i}}$ or its condensed counterpart are not invertible. Such situations occur in practice, essentially on very coarse meshes.

Adaptive time stepping. The previous strategies can lead to a solution failure for some given time step $\delta t$. If the Newton method did not converge, we try to compute the solution for a smaller time step $\delta t/2$. On the other hand, if the method did converge, we use for the subsequent time step the larger value $2\delta t$. The maximal time step allowed is the initial one, denoted by $\Delta t$. In practice, the scheme may perform numerous time step reductions at the beginning (early times) of the computation.

Initialisation by truncation and filtration. As for any Newton method, the question of the initialisation is fundamental in order to get a robust implementation. It appears that, for $n \geq 1$, the natural initialisation $\underline{\ell}_{ \mathcal{D}, 0} = \underline{\ell}_ \mathcal{D}^n$ is satisfactory when used with the adaptative time stepping strategy. However, for $n = 0$, such a choice is not possible, since $\underline{\ell}_ \mathcal{D}^0$ does not exist in general if $u^{0}$ vanishes locally or is too small (cf. Remark 6). A first way of tackling this problem is to define a truncated initial logarithm potential $\tilde{\ell}^{0}$ as

and to initialise the Newton method with $\tilde{ \underline{\ell}}_ \mathcal{D}^{0} = \underline{I}^k_ \mathcal{D}(\tilde{\ell}^0) \in \underline{V}_ \mathcal{D}^{k}$, provided one can give a sense to the face components. In fact, such a strategy entails another limitation: $\tilde{ \underline{\ell}}_ \mathcal{D}^{0}$ exhibits strong oscillations in the regions where the truncation is performed (this is also true when $u^{0}$ is discontinuous over $\Omega$, as in Section 4.2). These oscillations usually make the method diverge, even with extremely small time step. Therefore, we eventually initialise the method with a "filtered" (non-oscillating) discrete logarithm potential $\underline{\ell}^{0}_ \mathcal{D} \in \underline{V}_ \mathcal{D}^{k}$, which corresponds to a zero-order polynomial projection of $\tilde{\ell}^0$:

still provided one can give a sense to the face components. In practice, using $\underline{\ell}_{ \mathcal{D}, 0} = \underline{\ell}^{0}_ \mathcal{D}$ as the first initialisation (when $n = 0$) yields convergent Newton methods for all tests presented below. The use of filtered initial discrete data seems particularly crucial for high-order schemes ($k \geq 1$). For the lowest-order version of the scheme ($k = 0$), the use of $\tilde{ \underline{\ell}}_ \mathcal{D}^{0}$ as a first initialisation (for $n = 0$) often yields convergent Newton methods (up to time step reduction).

Of course, the chosen values for the stopping criterion and the thresholds are arbitrary and could be modified. However, the set of values advocated here makes the scheme robust enough so as to be capable of computing solutions for all the test-cases in this article.

Remark 12 (Potentials vs. densities). One of the main differences between the present nonlinear scheme and the low-order HFV ones from ^[16,45] lies in the fact that we use here the potential $\ell$ as our (piecewise polynomial) unknown, whereas the density $u$ was used in the low-order schemes. Notice that, in the present context, choosing $u$ as the main variable would require to give a discrete meaning to $\nabla\log(u)$, which is not obvious for the following reason: polynomials of degree $\geq 1$ are not stable by the $\log$ function. As a by-product of seeking for a potential, our stopping criterion only provides information on $\ell$, while we are eventually interested in the corresponding density. Moreover, our criterion only takes into account the coefficients of the polynomials (through the use of the norm $\|{\cdot}\|_{l^\infty}$), but such a measure does not give much information about the effective behaviour of the unknowns. A more relevant stopping criterion could be to consider the residual in terms of densities

(and analogous definition for the face unknowns) in order to ensure a satisfying accuracy on $u$. However, the main drawback of such a criterion is its evaluation cost. In this work, we thus chose to use instead a purely algebraic stopping criterion on $\ell$, whose cost is marginal. The testing of other stopping criteria will be the subject of future investigations. Last, for the HFV schemes of ^[16,45], the following loop-breaking strategy was used: when the computed density had almost-zero (or even non-positive) components, one performed a time step reduction. Here, such a situation cannot occur, since $\ell$ is authorised to take any real value, but this apparent latitude on the potential is in fact pernicious. Indeed, situations in which $\ell$ takes negative values with large magnitude are actually the counterpart of an almost-zero $u$ for the HFV schemes. Like their counterpart, they lead to divergent Newton methods. The main difficulty then lies in the design of a relevant criterion in order to avoid these situations.

4.2.   Positivity

This first section is dedicated to assessing discrete positivity preservation. For the test considered here, we set the advective potential and the anisotropy tensor to

For the initial datum, we take $u^{0} = 10^{-3} \, {\mathbb{1}}_{B} + {\mathbb{1}}_{\Omega \setminus B}$, where $B$ is the Euclidean ball

These data ensure that the solution $u$ is positive on $\mathbb{R}_{+}\times\Omega$. We perform the simulation on the time interval $[0, 5. 10^{-4}]$ with $\Delta t = 10^{-5}$, on a (fine) tilted hexagonal-dominant mesh featuring 4192 cells and 12512 edges. The computed discrete densities are denoted by $(\underline{\mathfrak{u}}^n_ \mathcal{D})_{1 \leq n \leq N_f}$ and $(\underline{\mathfrak{u}}^{\omega, n}_ \mathcal{D})_{1 \leq n \leq 50}$. Remark that the situation $N_f > 50$ may occur if the nonlinear scheme has to perform time step adaptation.

In Table 1, we collect the minimal values reached by the discrete solutions. The values of $\texttt{mincellA}$ (for "average") are defined by

for, respectively, the nonlinear scheme and the exponential fitting scheme. The values of $\texttt{mincellQN}$ are the minimal values taken by the densities at the cell quadrature nodes. Analogous definitions hold for the faces. The values of $\texttt{#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$, so it is not a relevant information to compare the cost of the schemes for different values of $k$. Last, $\texttt{walltime}$ is the total time (in $s$) needed to compute the discrete solution (it includes the pre-computation steps, such as the computation of the matrices representing $G_K$).

Recall that the exponential fitting scheme is linear (with corresponding matrix not depending on time), whence its extremely low cost compared to the nonlinear scheme. Note, however, that when an LU decomposition is unaffordable and an iterative solver has to be used instead, $\texttt{nlhho_k}$ is approximately "only" five times more costly than $\texttt{expf_k}$. The results of Table 1 first indicate that, as expected, all nonlinear schemes preserve the positivity of the discrete solution. On the other hand, none of the linear schemes preserves positivity on the whole domain $\Omega$. In fact, except $\texttt{expf_0}$, all linear schemes studied here do not even preserve the average positivity on each cell, in the sense that there exists $K_0 \in \mathcal{M}$ and $n_0$ integer such that

Moreover, it is interesting to note that the positivity violation peak (which can be approximated by $|\texttt{mincellQN}|$ and $|\texttt{minfaceQN}|$) increases as $k$ increases, whereas in average (values of $|\texttt{mincellA}|$ and $|\texttt{minfaceA}|$) the lack of positivity becomes smaller as the order increases.

At this stage, it is worth pointing out the fact that quantifying the negativity of the solution is much more difficult for high-order schemes, since it is not possible to "count" the number of negative values (which correspond to the degrees of freedom for low-order schemes). While the $\texttt{mincellQN}$ value gives information about the minimum value reached on the whole domain, it does not give any indication about the measure of the set $\lbrace x\in \overline{\Omega} \mid \mathfrak{u}^{\omega, n}_ \mathcal{M}(x) < 0\rbrace$ where the discrete cell unknown takes negative values. The same remark applies to $\texttt{mincellA}$. As an attempt to provide an idea of the size of this set, we display in Table 2 the number of cells with negative average over the whole simulation, defined as the cardinal of the set

These data reveal that, excluding $\texttt{expf_0}$ which performs quite well on this particular test, the higher the order, the smaller the size of the negative-average set.

The previous observations seem to indicate a competition between two phenomena for linear methods. As $k$ increases, the accuracy is improved, and therefore the discrete solution becomes closer to the exact one. Hence, in average, high-order schemes compute solutions with smaller area of negativity, and lesser positivity violation. However, high values of $k$ induce larger oscillations for the polynomial solution: {t}he computed solution takes negative values on smaller sets, but the (pointwise) undershoots become bigger as $k$ increases. At the end, it seems that there is no hope to get a positive discrete solution on the whole domain $\Omega$ with a linear method.

Remark 13 (An accuracy criterion taking into account positivity). The previous observations suggest that, for applications in which preserving the positivity of the solution is an essential feature, the accuracy of the scheme should not simply be defined as an $L^p$-distance between the reconstructed discrete solution $\mathfrak{u}_ \mathcal{M}$ and the exact one $u$. We believe that a relevant criterion in order to take into account both "classical accuracy" (distance between $\mathfrak{u}_ \mathcal{M}$ and $u$) and positivity is to look at the relative Boltzmann entropy (or other kinds of relative $\Phi$-entropies as defined in ^[6]) with respect to the exact solution, that is

where $\Phi_1(s) = s\log(s) -s +1$ for $s > 0$, $\Phi_1(0) = 1$, and $\Phi_1(s)$ takes large values for $s < 0$. The interest of such a definition is twofold. First, the negativity of $\mathfrak{u}_ \mathcal{M}$ is penalised. Second, if $\mathfrak{u}_ \mathcal{M}$ is positive and $\int_\Omega \mathfrak{u}_ \mathcal{M} = \int_\Omega u$ (which is the case in practice for problems with homogeneous Neumann boundary conditions), by Csiszár–Kullback inequality (see e.g., [11, Lemma 5.6]), one has

4.3.   Convergence and efficiency

We here study the convergence as $(h_ \mathcal{D}, \Delta t) \to (0, 0)$ of the nonlinear scheme for different values of the polynomial degree $k$. We consider a test-case with known exact solution. We set the advective potential and anisotropy tensor to

for $l_{x_1} > 0$. The exact solution is then given by

where $C_1 > 0$ and $\alpha = l_{x_1} \left (\frac{1}{4} + \pi^2 \right)$. Note that $u^{0}$ vanishes on $\{ x_1 = 1 \}$, but for any $t > 0$, $u(t, \cdot) > 0$. Here, our experiments are performed using $l_{x_1} = 1$ and $C_1 = 10^{-1}$.

We compute the discrete solutions on the time interval $[0, 0.1]$, and we denote by $(\underline{\mathfrak{u}}^n_ \mathcal{D})_{1 \leq n \leq N_f}$ the corresponding discrete densities. We monitor the relative $L^2_t(L^2_x)$-error and $L^2_t(H^1_x)$-error on the solution, respectively defined by

and

where $\delta t^n = t^{n}-t^{n-1}$, and $\sum_{1 \leq n \leq N_f} \delta t^n = 0.1$. The discrete gradient of the densities $\mathfrak{G}_ \mathcal{M}(\underline{\mathfrak{u}}^n_ \mathcal{D})$ is defined as follows. For all $K\in \mathcal{M}$, $(\mathfrak{G}_ \mathcal{M}(\underline{\mathfrak{u}}^n_ \mathcal{D}))_{\mid K} = \mathfrak{G}_K(\underline{\mathfrak{u}}^n_K)$, where $\mathfrak{G}_K(\underline{\mathfrak{u}}^n_K)$ is a smooth vector field on $K$ defined by mimicking at the discrete level the relation $\nabla u = {\rm{e}}^\ell \nabla \ell$:

Notice that, with the chosen error measures, we do not take into account the time $t = 0$. We perform our simulations on a triangular mesh family $(\mathcal{D}_i)_{1\leq i \leq 5}$, such that $h_{ \mathcal{D}_{i}} /h_{ \mathcal{D}_{i+1}} = 2$. Since the time discretisation is of order one, in $L^\infty_t(L^2_x)$-norm, we expect the error to decrease as

where $C_T, C_S(k) > 0$ are multiplicative constants respectively related to time and space discretisations, with $k\mapsto C_S(k)$ decreasing. We have $h_{ \mathcal{D}_i} = h_{ \mathcal{D}_1}/2^{i-1}$, so to balance the time and space contributions of the error upper bound, we need to take

where $\Delta t(1, k) = \frac{C_S(k)}{C_T}h_{ \mathcal{D}_1}^{k+2}$. For the values of $k$ we consider, we assume (and we verify in practice that it is relevant) that $\Delta t(1, k)\geq 0.05/2^{k+2}$ (see ^[1] for a theoretical study of $k\mapsto C_S(k)$ in the HHO context). Thus, for given $i$ and $k$, we define our (maximal) time step as

On Figure 1, we plot the relative errors as functions of the mesh size $h_ \mathcal{D}$ for $k\in\{0, 1, 2, 3\}$. For completeness, we also include the scheme $\texttt{nlhho_2_0}$ (i.e., with $\varepsilon = 0$ in (2.18) for $k = 2$) in our comparison. First, we observe that $\texttt{nlhho_2}$ and $\texttt{nlhho_2_0}$ have the same behaviour (the two plots are superimposed). Tests with other values of $k$, not shown here, indicate that the influence of $\varepsilon$ ($0$ or $1$) on the accuracy of the scheme is not noticeable. Second, we see that, as one could expect, the method $\texttt{nlhho_k}$ converges at order $k+2$ in $L^2_t(L^2_x)$-norm. In the $L^2_t(H^1_x)$-norm, if the expected convergence order of $k+1$ is attained for $k = 0$ and $k = 1$, then some sort of saturation appears for $k = 2$ and $k = 3$. Since this saturation does not show up in $L^2_t(L^2_x)$-norm, we suspect this might be due to our definition (4.4) of the discrete density gradient. Indeed, remark that, at the discrete level, the chain rule is violated, thus (4.4) is not exactly a discrete version of $\nabla{\rm{e}}^\ell$.

We now study efficiency, that is to say accuracy for a given computational cost. On Figure 2, we plot the relative errors as functions of the simulation walltime (in $s$). Here again, the results for the schemes $\texttt{nlhho_2}$ and $\texttt{nlhho_2_0}$ are superimposed. It is quite remarkable to observe that, even with a low-order discretisation in time, significant efficiency gains can be obtained using a larger value of $k$, at least for values of $k\leq 2$. The gain is expected to be even larger after parallelising the local computations. Of course, the use of higher-order time-stepping methods should also lead to significant gains of efficiency. This will be investigated in future works.

Remark 14 (High-order schemes in time and space). The extension of the nonlinear scheme (2.20) to arbitrary orders in time and space is a rather natural goal in order to achieve optimal efficiency. However, even with a time discretisation of order $2$ (like for example Crank–Nicolson, which is perhaps the most natural extension to backward Euler), there is currently no successful approach retaining the discrete entropy structure. Since this structure is the cornerstone of the analysis (including the existence of solutions), it is of utmost importance to preserve it. Some numerical investigations on nonlinear entropic TPFA schemes for diffusive problems with BDF2 time discretisation have been performed in [21, Chapter 3], and indicate that such a time discretisation could also be a good candidate, even in regard of long-time behaviour (see [21, Chapter 3.4.4]). An alternative approach is to consider space-time methods, as in ^[8] in the context of conforming Galerkin discretisations of cross-diffusion systems. The extension of space-time techniques to polytopal grids is currently an active research area.

For completeness, we finally perform simulations on distorted quadrangular meshes, and display the relative $L^2_t(L^2_x)$-errors on Figure 3. Note that we use the same time step definition as for the previous simulations, whereas the initial mesh is coarser. As expected, the behaviour of the schemes is not strongly impacted by the mesh geometry, and $\texttt{nlhho_k}$ converges at order $k+2$ in $L^2_t(L^2_x)$-norm. When it comes to efficiency, increasing the value of $k$ leads to better accuracy for fixed computational cost, but the efficiency gain saturates for $k\geq 2$. It is also worth noting that on the coarsest mesh, $\texttt{nlhho_3}$ has to perform more time step reductions than the other schemes, because at some iterations the linear solver is unable to perform LU decomposition. These time step reductions do not only occur at the beginning of the simulation, and are probably related to the bad conditioning of the system for high-order polynomials (we use here monomial bases). Based on these observations, using $\texttt{nlhho_2}$ seems to be a sound choice to optimise efficiency while ensuring a good numerical stability.

4.4.   Discrete long-time behaviour

We are now interested in the long-time behaviour of discrete solutions. We first use the same test-case as in Section 4.3, but this time with an anisotropic tensor: we set $l_{x_1} = 10^{-2}$. The corresponding steady-state is

We compute the discrete solutions on the time interval $[0,350]$, with $\Delta t = 10^{-1}$, on two Kershaw meshes of sizes $0.02$ and $0.006$. On Figure 4, we display the evolution along time of the $L^1$-distance between the reconstructed discrete densities and $u^\infty$, computed as

for, respectively, the exponential fitting scheme, and the nonlinear scheme. We here focus on $\texttt{expf_1}$, and on $\texttt{nlhho_k}$ for $k \in \lbrace 0, 1, 2 \rbrace$ (as well as on $\texttt{nlhho_1_0}$). For all schemes, we observe the exponential convergence towards the thermal equilibrium, until machine precision is reached. Remark that, for the test-case considered here, $\phi \in \mathbb{P}^1(\Omega)$, therefore $\phi_ \mathcal{M} = \phi$ for all $k\geq 0$. It follows that $\mathfrak{u}_ \mathcal{M}^\infty = u^\infty$ (recall that we always have $\mathfrak{u}^{\omega, \infty}_ \mathcal{M} = u^\infty$). This is exactly what we observe in the numerical experiments. As previously, $\texttt{nlhho_1}$ and $\texttt{nlhho_1_0}$ exhibit an extremely similar behaviour. Also, for $k \geq 1$, we observe that the decay rates are similar to the exact one $\alpha$, and do not seem to depend on the size of the mesh. For $k = 0$, the decay rate differs a bit from $\alpha$ on the coarsest mesh, but these two rates seem to coincide on a sufficiently refined mesh.

As a last test-case, we consider an advective potential and an anisotropy tensor set to

Our initial datum is

The corresponding thermal equilibrium therefore reads

Remark that the potential $\phi$ is not (piecewise) polynomial. As previously, we investigate the long-time behaviour of the schemes. We compute the discrete solutions on the time interval $[0, 5]$, with $\Delta t = 0.2$, on two distorted quadrangular meshes featuring, respectively, 64 and 1024 cells. On Figure 5, we display the evolution of the $L^1$-distance to equilibrium, as defined in (4.5), for both the $\texttt{expf_k}$ and $\texttt{nlhho_k}$ schemes, for $k\in\{0, 1, 2, 3\}$. For all schemes, we observe the exponential convergence towards the thermal equilibrium, until some precision is reached. For the exponential fitting schemes, machine precision is attained (which is expected since $\mathfrak{u}^{\omega, \infty}_ \mathcal{M} = u^\infty$), whereas for the nonlinear schemes (for which $\mathfrak{u}^\infty_ \mathcal{M}$ is an approximation of $u^\infty$), the precision increases, as expected, with the polynomial degree and as the mesh is refined. Also, all schemes with $k\geq 1$ seem to exhibit a similar, meshsize-independent decay rate. For $k = 0$, the decay rate seems slightly sensitive to the mesh size, but tends to reach the expected value on a sufficiently refined mesh.

5.   Conclusions

In this paper, we have studied two arbitrary-order hybrid methods for the approximation of linear, anisotropic, potential-driven advection-diffusion equations on general polytopal meshes. The first one is a linear scheme, which is based on the exponential fitting strategy, whereas the second is a nonlinear scheme, whose building principles are adapted from the low-order constructions of ^[12,13,16]. We proved that both schemes admit solutions, possess a discrete entropy structure, and preserve the mass, the thermal equilibrium, and the long-time asymptotics. Moreover, the solutions to the nonlinear scheme are positive by construction. We have validated these theoretical results on a set of numerical test-cases. We have unraveled the positivity violation of the linear methods, which justifies the use of (more costly) nonlinear methods. In the meantime, the use of nonlinear schemes with polynomial unknowns of higher degree results in an important gain of efficiency (accuracy vs. computational cost). These results confirm the benefits of using high-order nonlinear schemes in order to get reliable approximations of dissipative problems. Future research directions include a full analysis of the nonlinear scheme, in particular of its convergence (with respect to the discretisation parameters) and time-asymptotic properties, as well as the development of similar schemes for more complex, nonlinear problems, like semiconductor models (based on ^[45]).

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

The authors would like to thank the anonymous reviewers for their remarks and suggestions which helped improving the quality of the presentation. The authors also thank J. Droniou for his insightful comments about this work, and for pointing out a simplification of the proof of Lemma 1. The authors finally thank C. Chainais-Hillairet and M. Herda for fruitful discussions on the topic. This research was funded in part by the Austrian Science Fund (FWF) project 10.55776/F65. The authors also acknowledge support from the LabEx CEMPI (ANR-11-LABX-0007).

Conflict of interest

The authors declare no conflicts of interest.