Numerical analysis of a nonlinear free-energy diminishing Discrete Duality Finite Volume scheme for convection diffusion equations

TL;DR

This paper introduces a nonlinear Discrete Duality Finite Volume scheme for convection diffusion equations that preserves energy dissipation at the discrete level, ensuring accurate long-term behavior and convergence to weak solutions.

Contribution

The paper develops a novel nonlinear scheme that maintains energy relations on distorted meshes and proves its convergence and positivity properties.

Findings

Scheme preserves energy dissipation even on distorted meshes

Convergence of the scheme to weak solutions is established

Numerical experiments confirm good long-time behavior

Abstract

We propose a nonlinear Discrete Duality Finite Volume scheme to approximate the solutions of drift diffusion equations. The scheme is built to preserve at the discrete level even on severely distorted meshes the energy / energy dissipation relation. This relation is of paramount importance to capture the long-time behavior of the problem in an accurate way. To enforce it, the linear convection diffusion equation is rewritten in a nonlinear form before being discretized. We establish the existence of positive solutions to the scheme. Based on compactness arguments, the convergence of the approximate solution towards a weak solution is established. Finally, we provide numerical evidences of the good behavior of the scheme when the discretization parameters tend to 0 and when time goes to infinity.

Tables

Table 1. Table 1 . Numerical results on the Quadrangle mesh family, final time T=0.25.

M	dt	$\kappa =0$		$\kappa ={10}^{-1}$
		normU	ordU	normU	ordU
1	4.032E-03	1.798E-01	—	1.796E-01	—
2	1.008E-03	9.316E-02	0.95	9.313E-02	0.94
3	2.520E-04	4.717E-02	1.03	4.716E-02	1.03
4	6.300E-05	2.361E-02	1.05	2.361E-02	1.05
5	1.575E-05	1.135E-02	0.87	1.135E-02	0.87

Table 2. Table 2 . Numerical results on the Kershaw mesh family, final time T=0.25.

M	dt	errgu	ordgu	erru	ordu	$N_{\text{max}}$	$N_{\text{mean}}$	Min $u^n$
1	2.0E-03	6.693E-02	—	7.254E-03	—	9	2.15	1.010E-01
2	5.0E-04	2.353E-02	1.54	1.751E-03	2.09	8	2.02	2.582E-02
3	1.25E-04	1.235E-02	1.61	7.237E-04	2.20	7	1.49	6.488E-03
4	3.125E-05	7.819E-03	1.60	3.962E-04	2.11	7	1.07	1.628E-03
5	3.125E-05	5.507E-03	1.58	2.556E-04	1.98	7	1.04	1.628E-03

Table 3. Table 3 . Numerical results on the Quadrangle mesh family, final time T=0.25.

M	dt	errgu	ordgu	erru	ordu	$N_{\text{max}}$	$N_{\text{mean}}$	Min $u^n$
1	4.032E-03	1.754E-01	—	2.149E-02	—	9	2.26	1.803E-01
2	1.008E-03	5.933E-02	1.56	5.055E-03	2.08	9	2.04	5.079E-02
3	2.520E-04	2.294E-02	1.44	1.299E-03	2.06	8	1.96	1.352E-02
4	6.300E-05	8.631E-03	1.48	3.256E-04	2.09	8	1.22	3.349E-03
5	1.250E-05	2.715E-03	1.37	7.702E-05	1.70	7	1.01	8.695E-04

Equations

$${\partial }_{t}u+\mathrm{div}\mathbf{J}=0,\text{mbox}in{Q}_T=\Omega \times (0,T),$$

$$\mathbf{J}=-\mathbf{\Lambda }\nabla u-u\mathbf{\Lambda }\nabla V,\text{mbox}in{Q}_T,$$

$$\mathbf{J}\cdot \mathbf{n}=0,\text{mbox}on\partial \Omega \times (0,T),$$

$$u(\cdot ,0)=u_0,\text{mbox}in\Omega ,$$

$${\int }_{\Omega }{u}_0\mathrm{d}\mathbf{x}>0\text{mbox}and{\int }_{\Omega }H(u_0)\mathrm{d}\mathbf{x}<\infty ,$$

$${\lambda }_m|\mathbf{v}{|}^2\le \mathbf{\Lambda }(\mathbf{x})\mathbf{v}\cdot \mathbf{v}\le {\lambda }^M|\mathbf{v}{|}^2,\text{for all v∈R2 and almost all x∈Ω}.$$

$$\mathbf{J}=-u\mathbf{\Lambda }\nabla (\log u+V).$$

$$\frac{d\mathbb{E}}{dt}+\mathbb{I}=0,$$

$$\mathbb{E}(t)=$$

$$\mathbb{I}(t)=$$

$$u_{\infty }=\left({\int }_{\Omega }{u}_{0}dx/{\int }_{\Omega }{e}^{-V}dx\right)e^{-V}$$

$${\iint }_{Q_T}u|\nabla \log (u){|}^2\mathrm{d}\mathbf{x}\mathrm{d}t=4{\iint }_{Q_T}|\nabla \sqrt{u}{|}^2\mathrm{d}\mathbf{x}\mathrm{d}t\le C.$$

$$\left\{u:Q_T\to {\mathbb{R}}_{+}|{\int }_{\Omega }H(u(\mathbf{x},\cdot ))\mathrm{d}\mathbf{x}\in L^{\infty }(0,T)\text{and}\sqrt{u}\in L^2(0,T;H^1(\Omega ))\right\}.$$

$${\iint }_{Q_T}u{\partial }_t\phi \mathrm{d}\mathbf{x}\mathrm{d}t+{\int }_{\Omega }{u}_0\phi (\cdot ,0)\mathrm{d}\mathbf{x}-{\iint }_{Q_T}\left(u\nabla V+\nabla u\right)\cdot \mathbf{\Lambda }\nabla \phi \mathrm{d}\mathbf{x}\mathrm{d}t=0.$$

$${\mathrm{m}}_{\mathcal{D}}=\frac{1}{2}{\mathrm{m}}_{\sigma }{\mathrm{m}}_{{\sigma }^{\ast }}\sin ({\alpha }_{\mathcal{D}}),\forall \mathcal{D}={\mathcal{D}}_{\sigma ,{\sigma }^{\ast }}\in \mathfrak{D}.$$

$${\theta }_{\mathcal{D}}=\frac{1}{2\sin ({\alpha }_{\mathcal{D}})}\left(\frac{{\mathrm{m}}_{\sigma }}{{\mathrm{m}}_{{\sigma }^{\ast }}}+\frac{{\mathrm{m}}_{{\sigma }^{\ast }}}{{\mathrm{m}}_{\sigma }}\right)\ge 1,{\tilde{\theta }}_{\mathcal{D}}=\max \left(\underset{K\in {\mathfrak{M}}_{\mathcal{D}}}{\max }\frac{{\mathrm{m}}_{\mathcal{D}}}{{\mathrm{m}}_{\mathcal{D}\cap K}};\underset{K^{\ast }\in {\mathfrak{M}}_{\mathcal{D}}^{\ast }}{\max }\frac{{\mathrm{m}}_{\mathcal{D}}}{{\mathrm{m}}_{\mathcal{D}\cap K^{\ast }}}\right).$$

$$1\le {\theta }_{\mathcal{D}},{\tilde{\theta }}_{\mathcal{D}}\le {\theta }^{\star },\forall \mathcal{D}\in \mathfrak{D}.$$

$$\sin ({\alpha }_{\mathcal{D}})\ge \frac{1}{{\theta }^{\star }},\forall \mathcal{D}\in \mathfrak{D}.$$

$$\sum _{\mathcal{D}\in {\mathfrak{D}}_K}{\mathrm{m}}_{\mathcal{D}}\le {\theta }^{\star }{\mathrm{m}}_K\text{and}\sum _{\mathcal{D}\in {\mathfrak{D}}_{K^{\ast }}}{\mathrm{m}}_{\mathcal{D}}\le {\theta }^{\star }{\mathrm{m}}_{K^{\ast }}$$

$$u_{\mathcal{T}}\in {\mathbb{R}}^{\mathcal{T}}⟺u_{\mathcal{T}}=\left({\left(u_K\right)}_{\mathcal{K}\in \overline{\mathfrak{M}}},{\left(u_{K^{\ast }}\right)}_{{\mathcal{K}}^{\ast }\in {\overline{\mathfrak{M}}}^{\ast }}\right)$$

$${\mathbf{\xi }}_{\mathfrak{D}}\in {\left({\mathbb{R}}^2\right)}^{\mathfrak{D}}⟺{\mathbf{\xi }}_{\mathfrak{D}}={\left({\mathbf{\xi }}_{\mathcal{D}}\right)}_{\mathcal{D}\in \mathfrak{D}}.$$

$$\displaystyle\left\llbracket v_{\mathcal{T}},u_{\mathcal{T}}\right\rrbracket_{\mathcal{T}}$$

$${\left({\mathbf{\xi }}_{\mathfrak{D}},{\mathbf{\phi }}_{\mathfrak{D}}\right)}_{\mathbf{\Lambda },\mathfrak{D}}$$

$${\mathbf{\Lambda }}^{\mathcal{D}}=\frac{1}{{\mathrm{m}}_{\mathcal{D}}}{\int }_{\mathcal{D}}\mathbf{\Lambda }(\mathbf{x})\mathrm{d}\mathbf{x},\forall \mathcal{D}\in \mathfrak{D}.$$

$${\Vert {\mathbf{\xi }}_{\mathfrak{D}}\Vert }_{\mathbf{\Lambda },\mathfrak{D}}^2={\left({\mathbf{\xi }}_{\mathfrak{D}},{\mathbf{\xi }}_{\mathfrak{D}}\right)}_{\mathbf{\Lambda },\mathfrak{D}},\forall {\mathbf{\xi }}_{\mathfrak{D}}\in {\left({\mathbb{R}}^2\right)}^{\mathfrak{D}}.$$

$${\nabla }^{\mathcal{D}}{u}_{\mathcal{T}}=\frac{1}{\sin ({\alpha }_{\mathcal{D}})}\left(\frac{u_L-u_K}{{\mathrm{m}}_{{\sigma }^{\ast }}}{\mathbf{n}}_{\sigma K}+\frac{u_{L^{\ast }}-u_{K^{\ast }}}{{\mathrm{m}}_{\sigma }}{\mathbf{n}}_{{\sigma }^{\ast }{K}^{\ast }}\right),\forall \mathcal{D}\in \mathfrak{D}.$$

$${\nabla }^{\mathcal{D}}{u}_{\mathcal{T}}=\frac{1}{2{\mathrm{m}}_{\mathcal{D}}}\left({\mathrm{m}}_{\sigma }(u_L-u_K){\mathbf{n}}_{\sigma K}+{\mathrm{m}}_{{\sigma }^{\ast }}(u_{L^{\ast }}-u_{K^{\ast }}){\mathbf{n}}_{{\sigma }^{\ast }{K}^{\ast }}\right),\forall \mathcal{D}\in \mathfrak{D}.$$

$${\mathrm{div}}^{\mathcal{T}}{\mathbf{\xi }}_{\mathfrak{D}}=\left({\mathrm{div}}^{\mathfrak{M}}{\mathbf{\xi }}_{\mathfrak{D}},{\mathrm{div}}^{\partial \mathfrak{M}}{\mathbf{\xi }}_{\mathfrak{D}},{\mathrm{div}}^{{\mathfrak{M}}^{\ast }}{\mathbf{\xi }}_{\mathfrak{D}},{\mathrm{div}}^{\partial {\mathfrak{M}}^{\ast }}{\mathbf{\xi }}_{\mathfrak{D}}\right),$$

$$\forall K\in \mathfrak{M},{\mathrm{div}}_K{\mathbf{\xi }}_{\mathfrak{D}}=\frac{1}{{\mathrm{m}}_K}{\sum _{\mathcal{D}\in {\mathfrak{D}}_K}}_{\mathcal{D}={\mathcal{D}}_{\sigma ,{\sigma }^{\ast }}}{\mathrm{m}}_{\sigma }{\mathbf{\xi }}_{\mathcal{D}}\cdot {\mathbf{n}}_{\sigma K},$$

$${\mathcal{P}}^{\mathcal{T}}{u}_{\mathcal{T}}=\left({\mathcal{P}}^{\mathfrak{M}}{u}_{\mathcal{T}},{\mathcal{P}}^{\partial \mathfrak{M}}{u}_{\mathcal{T}},{\mathcal{P}}^{{\mathfrak{M}}^{\ast }}{u}_{\mathcal{T}},{\mathcal{P}}^{\partial {\mathfrak{M}}^{\ast }}{u}_{\mathcal{T}}\right),$$

$$[[{\mathcal{P}}^{\mathcal{T}}{u}_{\mathcal{T}},v_{\mathcal{T}}{]]}_{\mathcal{T}}=\frac{1}{2}\frac{1}{\mathrm{Size}(\mathcal{T}{)}^{\beta }}\sum _{K^{\ast }\in {\overline{\mathfrak{M}}}^{\ast }}\sum _{K\in \mathfrak{M}}{\mathrm{m}}_{K\cap K^{\ast }}(u_K-u_{K^{\ast }})(v_K-v_{K^{\ast }})\text{mbox}with\beta \in (0,2).$$

$$r^{\mathcal{D}}(u_{\mathcal{T}})=\frac{1}{4}(u_K+u_L+u_{K^{\ast }}+u_{L^{\ast }}).$$

$$\delta^{\mathcal{D}}u_{\mathcal{T}}=\left(\begin{array}[]{c}u_{K}-u_{L}\\ u_{{{K}^{*}}}-u_{{{L}^{*}}}\end{array}\right).$$

Full text

The authors are supported by the Inria teams RAPSODI and COFFEE, the LabEx CEMPI (ANR-11-LABX-0007-01), the GEOPOR project (ANR-13-JS01-0007-01) and the MOONRISE project (ANR-14-CE23-0007).

Keywords. Convection diffusion equation, Discrete Duality Finite Volumes, convergence analysis, discrete entropy method

**AMS subjects classification. ** 65M08, 65M12, 35K20

1. Introduction

1.1. Motivation

The modeling of systems of interacting particles, like electrons in electronic devices, ions in plasmas, chemical species in biological membranes for instance, leads to systems of evolutive partial differential equations. The knowledge of the large time behavior of such systems is crucial for the understanding of the underlying physical phenomena. In many cases, the relaxation to an equilibrium configuration is based on the second law of thermodynamics and on the dissipation of some entropies.

Based on works in kinetic theory, mathematicians have intensively developed the entropy method for the study of the large time behavior of different systems of PDEs. Let us mention works on Boltzmann and Landau equations [43], on linear Fokker-Planck equations [15], on porous media equations [16], on reaction-diffusion systems [23, 24, 34], on drift-diffusion systems for semiconductor devices [31, 32, 30]. We also refer to the survey paper [5] and to the reference book [37]. Similar results were obtained based on the interpretation of PDE models as Wasserstein gradient flows [1]. We refer for instance to [36, 10] for linear Fokker Planck equations, to [40] for the porous medium equation, to [11] for granular media. This list is far from being exhaustive.

The knowledge of the large time behavior of such evolution equations, the existence of some entropies which are dissipated along time are structural features, as positivity of densities or conservation of mass, that should be preserved at the discrete level by numerical schemes. The question of the large time behavior of numerical schemes has been investigated for instance for coagulation-fragmentation models [29], for nonlinear diffusion equations [17, 38], for reaction-diffusion systems [34, 35], for drift-diffusion systems [19, 7]. These last works show that the Scharfetter-Gummel numerical fluxes, first introduced in [42] for the approximation of convection-diffusion fluxes and widely used later for the simulation of semiconductor devices, preserve the thermal equilibrium. Their use in the numerical approximation of drift-diffusion systems ensure the exponential decay towards equilibrium of the numerical scheme. Unfortunately, such numerical fluxes can only be applied in two-points flux approximation finite volume schemes and therefore on restricted meshes. Moreover, they do not extend to anisotropic convection-diffusion equations.

Therefore, it seems crucial to propose new finite volume schemes which preserve the large-time behavior of anisotropic convection-diffusion equations and which apply on almost general meshes. In [14], the authors proposed and analyzed a VAG scheme satisfying these prescribed properties. In this work, we propose and study the convergence analysis of a nonlinear free-energy diminishing discrete duality finite volume scheme [12].

1.2. Presentation of the continuous problem

We focus on a very basic drift-diffusion equation with potential convection and anisotropy. Let $\Omega$ be a polygonal connected open bounded subset of $\mathbb{R}^{2}$ and let $T>0$ be a finite time horizon. The problem writes:

[TABLE]

with ${\mathbf{n}}$ the outward unit normal to ${\partial}\Omega$ and the following assumptions on the data:

1. (A1)

The initial data $u_{0}$ is measurable, nonnegative and satisfies

[TABLE]

where $H(s)=s\log s-s+1$ for all $s\geq 0$. 2. (A2)

The exterior potential $V$ belongs to $C^{1}({\overline{\Omega}},\mathbb{R})$. Without loss of generality, we assume that $V\geq 0$ in $\overline{\Omega}$. 3. (A3)

The anisotropy tensor ${\boldsymbol{\Lambda}}$ is supposed to be bounded (i.e., ${\boldsymbol{\Lambda}}\in L^{\infty}(\Omega)^{2\times 2}$), symmetric (i.e., ${\boldsymbol{\Lambda}}={\boldsymbol{\Lambda}}^{T}$ a.e. in $\Omega$), and uniformly elliptic: there exist $\lambda_{m}>0$ and $\lambda^{M}>0$ such that

[TABLE]

The flux ${\mathbf{J}}$ can be reformulated in the nonlinear form

[TABLE]

Testing equation (1a) by $\log(u)+V$ leads to the so-called energy/energy dissipation relation (energy/dissipation for short)

[TABLE]

where the free energy $\mathbb{E}$ and the dissipation $\mathbb{I}$ for (1) are respectively defined by

[TABLE]

Since $u$ is nonnegative, so does ${\mathbb{I}}$ and the free energy $\mathbb{E}$ is decaying with time. As highlighted for instance in [6] and [10], the solution $u$ to (1) converges towards the steady-state

[TABLE]

when time goes to infinity. In the case where ${\boldsymbol{\Lambda}}$ does not depend on ${\boldsymbol{x}}$ and where both $\Omega$ and $V$ are convex, this convergence is exponentially fast.

The energy/energy dissipation relation (4) provides a control on the Fisher information

[TABLE]

Thus it is natural to seek the solution in the space

[TABLE]

This motivates the following notion of weak solution.

Definition 1.1**.**

A function $u:Q_{T}\to\mathbb{R}_{+}$ is said to be a weak solution to the problem (1) if $H(u)\in L^{\infty}(0,T,L^{1}(\Omega))$, $\sqrt{u}\in L^{2}(0,T,H^{1}(\Omega))$, and (1) is satisfied in the distributional sense, i.e., for all $\varphi\in C^{\infty}_{c}(\overline{\Omega}\times[0,T))$, there holds

[TABLE]

1.3. Outline of the paper

In Section 2, we introduce the numerical scheme and state the main results of the paper: existence of a positive solution to the scheme and convergence of a sequence of approximate solutions towards a weak solution. The existence of a solution to the scheme is established in Section 3. It strongly relies on the conservation of mass at the discrete level and on a discrete counterpart of an energy/dissipation estimate. Section 4 is devoted to the proof of convergence of the scheme. The effective behavior of the numerical method is eventually discussed in Section 5. It is shown that the method is second order accurate w.r.t. space in $L^{2}$ norm, whereas the approximate gradient super-converges with observed order 3/2. Moreover, the method exhibit a very accurate long-time behavior.

2. Presentation of the scheme and main results

2.1. Meshes and notations

In order to define a DDFV scheme, as for instance in [25, 3], we need to introduce three different meshes – the primal mesh, the dual mesh and the diamond mesh – and some associated notations.

The primal mesh denoted $\overline{\mathfrak{M}}$ is composed of the interior primal mesh $\mathfrak{M}$ (a partition of $\Omega$ with polygonal control volumes) and the set $\partial\mathfrak{M}$ of boundary edges seen as degenerate control volumes. For all $K\in\overline{\mathfrak{M}}$, we define $x_{K}$ the center of $K$. The family of centers is denoted by ${\mathfrak{X}}=\{x_{{{K}}},{{K}}\in\overline{\mathfrak{M}}\}$.

Let ${\mathfrak{X}}^{*}$ denote the set of the vertices of the primal control volumes in $\overline{\mathfrak{M}}$. Distinguishing the interior vertices from the vertices lying on the boundary, we split ${\mathfrak{X}}^{*}$ into ${\mathfrak{X}}^{*}={\mathfrak{X}}_{int}^{*}\cup{\mathfrak{X}}_{ext}^{*}$. To any point $x_{{{K}^{*}}}\in{\mathfrak{X}}_{int}^{*}$, we associate the polygon ${{K}^{*}}$, whose vertices are $\{x_{{{K}}}\in{\mathfrak{X}}/x_{{{K}^{*}}}\in{\overline{{K}}},{{K}}\in\mathfrak{M}\}$. The set of these polygons defines the interior dual mesh denoted by $\mathfrak{M}^{*}$. To any point $x_{{{K}^{*}}}\in{\mathfrak{X}}_{ext}^{*}$, we then associate the polygon ${{K}^{*}}$, whose vertices are $\{x_{{{K}^{*}}}\}\cup\{x_{{{K}}}\in X/x_{{{K}^{*}}}\in{\bar{{K}}},{{K}}\in\overline{\mathfrak{M}}\}$. The set of these polygons is denoted by $\partial\mathfrak{M}^{*}$ called the boundary dual mesh and the dual mesh is $\mathfrak{M}^{*}\cup\partial\mathfrak{M}^{*}$, denoted by $\overline{\mathfrak{M}^{*}}$.

For all neighboring primal cells ${{K}}$ and ${{L}}$, we assume that $\partial{{K}}\cap\partial{{L}}$ is a segment, corresponding to an edge of the mesh $\mathfrak{M}$, denoted by $\sigma={{K}}|{{L}}$. Let $\mathcal{E}$ be the set of such edges. We similarly define the set $\mathcal{E}^{*}$ of the edges of the dual mesh. For each couple $(\sigma,\sigma^{*})\in\mathcal{E}\times\mathcal{E}^{*}$ such that $\sigma={{K}}|{{L}}=(x_{{{K}^{*}}},x_{{{L}^{*}}})$ and $\sigma^{*}={{K}^{*}}|{{L}^{*}}=(x_{{{K}}},x_{{{L}}})$, we define the quadrilateral diamond cell ${\mathcal{D}_{\sigma,\sigma^{*}}}$ whose diagonals are $\sigma$ and $\sigma^{*}$, as shown on Figure 1. If $\sigma\in\mathcal{E}\cap\partial\Omega$, we note that the diamond degenerates into a triangle. The set of the diamond cells defines the diamond mesh $\mathfrak{D}$. It is a partition of $\Omega$. We can rewrite $\mathfrak{D}=\mathfrak{D}^{ext}\cup\mathfrak{D}^{int}$ where $\mathfrak{D}^{ext}$ is the set of all the boundary diamonds and $\mathfrak{D}^{int}$ the set of all the interior diamonds.

Finally, the DDFV mesh is made of $\mathcal{T}=(\overline{\mathfrak{M}},\overline{\mathfrak{M}^{*}})$ and $\mathfrak{D}$. For each primal or dual cell $M$ ($M\in{\overline{\mathfrak{M}}}$ or $M\in{\overline{\mathfrak{M}}^{*}}$), we define ${\rm m}_{M}$ the measure of $M$, $\mathcal{E}_{M}$ the set of the edges of $M$ (it coincides with the edge $\sigma=M$ if $M\in\partial\mathfrak{M}$), $\mathfrak{D}_{M}$ the set of diamonds ${\mathcal{D}_{\sigma,\sigma^{*}}}\in\mathfrak{D}$ such that ${\rm m}({\mathcal{D}_{\sigma,\sigma^{*}}}\cap M)>0$, and $d_{M}$ the diameter of $M$.

For a diamond ${\mathcal{D}_{\sigma,\sigma^{*}}}$, whose vertices are $(x_{{{K}}},x_{{{K}^{*}}},x_{{{L}}},x_{{{L}^{*}}})$, we define: $x_{\mathcal{D}}$ the center of the diamond cell ${\mathcal{D}}$: $\{x_{{\scriptscriptstyle\mathcal{D}}}\}=\sigma\cap\sigma^{*}$, ${\rm m}_{\sigma}$ the length of the primal edge $\sigma$, ${\rm m}_{\sigma^{*}}$ the length of the dual edge $\sigma^{*}$, $d_{\mathcal{D}}$ the diameter of ${\mathcal{D}}$, $\alpha_{\mathcal{D}}$ the angle between $(x_{{{K}}},x_{{{L}}})$ and $(x_{{{K}^{*}}},x_{{{L}^{*}}})$. We will also use two direct basis $(\boldsymbol{{\boldsymbol{\tau}}}_{{{{K}^{*}}\!\!\!,{{L}^{*}}}},{\mathbf{n}}_{\sigma K})$ and $(\boldsymbol{{\mathbf{n}}}_{\sigma^{*}{{K}^{*}}},\boldsymbol{{\boldsymbol{\tau}}}_{{K\!,L}})$, where ${\mathbf{n}}_{\sigma K}$ is the unit normal to $\sigma$, outward ${{K}}$, $\boldsymbol{{\mathbf{n}}}_{\sigma^{*}{{K}^{*}}}$ is the unit normal to $\sigma^{*}$, outward ${{K}^{*}}$, $\boldsymbol{{\boldsymbol{\tau}}}_{{{{K}^{*}}\!\!\!,{{L}^{*}}}}$ is the unit tangent vector to $\sigma$, oriented from ${{K}^{*}}$ to ${{L}^{*}}$, $\boldsymbol{{\boldsymbol{\tau}}}_{{K\!,L}}$ is the unit tangent vector to $\sigma^{*}$, oriented from ${{K}}$ to ${{L}}$. Denoting by ${\rm m}_{{\mathcal{D}}}$ the $2$-dimensional Lebesgue measure of ${\mathcal{D}}$, one has

[TABLE]

We define two local regularity factors $\theta_{\mathcal{D}},\tilde{\theta}_{\mathcal{D}}$ of the diamond cell ${\mathcal{D}}={\mathcal{D}}_{\sigma,\sigma^{\ast}} \in\mathfrak{D}$ by

[TABLE]

In what follows, we assume that there exists $\theta^{\star}\geq 1$ such that

[TABLE]

In particular, this implies that

[TABLE]

Moreover, owing to the definition of $\tilde{\theta}_{\mathcal{D}}$ and to (11), one has

[TABLE]

Finally, we define the size of the mesh: ${\rm size}(\mathcal{T})={\max_{{\mathcal{D}}\in\mathfrak{D}}}{\rm d}_{{\mathcal{D}}}.$

2.2. Discrete unknowns and discrete operators

We first define the different sets of discrete unknowns. As it is usual for DDFV methods, we need several types of degrees of freedom to represent scalar and vector fields in the discrete setting. We introduce $\mathbb{R}^{\mathcal{T}}$ the linear space of scalar fields constant on the cells of ${\overline{\mathfrak{M}}}$ and ${\overline{\mathfrak{M}}^{*}}$:

[TABLE]

and $\left(\mathbb{R}^{2}\right)^{\mathfrak{D}}$ the linear space of vector fields constant on the diamonds:

[TABLE]

Let us mention that we similarly denote by $\mathbb{R}^{\mathfrak{D}}$ the set of scalar fields constant on the diamonds.

Then,we define the positive semi-definite bilinear form111Although it mimics the continuous $L^{2}(\Omega)$ scalar product, the bilinear form $\llbracket\cdot,\cdot\rrbracket_{\mathcal{T}}$ is not a scalar product since it does not involve the primal boundary edges ${\partial}\mathfrak{M}$. It is therefore not definite. $\llbracket\cdot,\cdot\rrbracket_{\mathcal{T}}$ on $\mathbb{R}^{\mathcal{T}}$ and the scalar product $(\cdot,\cdot)_{{\boldsymbol{\Lambda}},\mathfrak{D}}$ on $\left(\mathbb{R}^{2}\right)^{\mathfrak{D}}$ by

[TABLE]

where

[TABLE]

We denote by $\|\cdot\|_{{\boldsymbol{\Lambda}},\mathfrak{D}}$ the Euclidian norm associated to the scalar product $\left(\cdot,\cdot\right)_{{\boldsymbol{\Lambda}},\mathfrak{D}}$, i.e.,

[TABLE]

The DDFV method is based on the definitions of a discrete gradient and of a discrete divergence, which are linked by duality formula as shown in [25]. The discrete gradient has been introduced in [20] and developed in [25]. It is a mapping from $\mathbb{R}^{\mathcal{T}}$ to $\left(\mathbb{R}^{2}\right)^{\mathfrak{D}}$ defined by $\nabla^{\mathfrak{D}}u_{\mathcal{T}}=\displaystyle\left(\nabla^{{\mathcal{D}}}u_{\mathcal{T}}\right)_{{\mathcal{D}}\in\mathfrak{D}}$ for all $u_{\mathcal{T}}\in\mathbb{R}^{\mathcal{T}}$, where

[TABLE]

Using (9), the discrete gradient can be equivalently written:

[TABLE]

The discrete divergence has been introduced in [25]. It is a mapping $\mathrm{div}^{\mathcal{T}}$ from $\left(\mathbb{R}^{2}\right)^{\mathfrak{D}}$ to $\mathbb{R}^{\mathcal{T}}$ defined for all ${\boldsymbol{\xi}}_{\mathfrak{D}}\in\left(\mathbb{R}^{2}\right)^{\mathfrak{D}}$ by

[TABLE]

with $\mathrm{div}^{\mathfrak{M}}{\boldsymbol{\xi}}_{\mathfrak{D}}=\left(\mathrm{div}_{{{K}}}{\boldsymbol{\xi}}_{\mathfrak{D}}\right)_{{\scriptscriptstyle\mathcal{K}}\in\mathfrak{M}}$, $\mathrm{div}^{\partial\mathfrak{M}}{\boldsymbol{\xi}}_{\mathfrak{D}}=0$, $\mathrm{div}^{\mathfrak{M}^{*}}{\boldsymbol{\xi}}_{\mathfrak{D}}=\left(\mathrm{div}_{{{K}^{*}}}{\boldsymbol{\xi}}_{\mathfrak{D}}\right)_{{\scriptscriptstyle\mathcal{K}^{*}}\in\mathfrak{M}^{*}}$ and $\mathrm{div}^{\partial\mathfrak{M}^{*}}{\boldsymbol{\xi}}_{\mathfrak{D}}=\left(\mathrm{div}_{{{K}^{*}}}{\boldsymbol{\xi}}_{\mathfrak{D}}\right)_{{\scriptscriptstyle\mathcal{K}^{*}}\in\partial\mathfrak{M}^{*}}$ such that:

[TABLE]

and analogous definitions for $\mathrm{div}_{{{K}^{*}}}{\boldsymbol{\xi}}_{\mathfrak{D}}$ for ${{K}^{*}}\in\overline{\mathfrak{M}}^{*}$.

In [2], the authors study the convergence of DDFV schemes for degenerate hyperbolic-parabolic problems. They show that a penalization operator is needed in order to establish the convergence proof. Indeed, this penalization operator ensures that the two components of a discrete function (reconstructions on the primal and dual meshes) converge to the same limit. For similar reasons (see Section 4), we consider the same penalization operator ${\mathcal{P}}^{\mathcal{T}}:\mathbb{R}^{\mathcal{T}}\to\mathbb{R}^{\mathcal{T}}$ as in [2] and [18]. It is defined for all $u_{\mathcal{T}}\in\mathbb{R}^{\mathcal{T}}$ by

[TABLE]

with $\mathcal{P}^{\mathfrak{M}}u_{\mathcal{T}}=\left(\mathcal{P}_{K}u_{\mathcal{T}}\right)_{K\in\mathfrak{M}}$, $\mathcal{P}^{\partial\mathfrak{M}}u_{\mathcal{T}}=0$, $\mathcal{P}^{\mathfrak{M}^{*}}u_{\mathcal{T}}=\left(\mathcal{P}_{{{K}^{*}}}u_{\mathcal{T}}\right)_{{{K}^{*}}\in\mathfrak{M}^{*}}$ and $\mathcal{P}^{\partial\mathfrak{M}^{*}}u_{\mathcal{T}}=\left(\mathcal{P}_{{{K}^{*}}}u_{\mathcal{T}}\right)_{{{K}^{*}}\in\partial\mathfrak{M}^{*}}$ such that, for a given parameter $\beta\in(0,2)$,

$\begin{aligned} \forall\ {{K}}\in\mathfrak{M},&\ \ \mathcal{P}_{K}u_{\mathcal{T}}=\displaystyle\frac{1}{{\rm m}_{{{K}}}}\frac{1}{\mathrm{size}(\mathcal{T})^{\beta}}\operatorname*{\sum}_{{{K}^{*}}\in{\overline{\mathfrak{M}}^{*}}}{\rm m}_{K\cap{{K}^{*}}}(u_{{{K}}}-u_{{{K}^{*}}}),\\ \forall\ {{K}^{*}}\in{\overline{\mathfrak{M}}^{*}},&\ \ \mathcal{P}_{{{K}^{*}}}u_{\mathcal{T}}=\displaystyle\frac{1}{{\rm m}_{{{K}^{*}}}}\frac{1}{\mathrm{size}(\mathcal{T})^{\beta}}\operatorname*{\sum}_{{{K}}\in\mathfrak{M}}{\rm m}_{K\cap{{K}^{*}}}(u_{{{K}^{*}}}-u_{{{K}}}).\end{aligned}$

It clearly satisfies : for all $u_{\mathcal{T}},v_{\mathcal{T}}\in\mathbb{R}^{\mathcal{T}}$,

[TABLE]

Finally, we introduce a reconstruction operator on diamonds $r^{\mathfrak{D}}$. It is a mapping from $\mathbb{R}^{\mathcal{T}}$ to $\mathbb{R}^{\mathfrak{D}}$ defined for all $u_{\mathcal{T}}\in\mathbb{R}^{\mathcal{T}}$ by $r^{\mathfrak{D}}[u_{\mathcal{T}}]=\left(r^{\mathcal{D}}(u_{\mathcal{T}})\right)_{{\mathcal{D}}\in\mathfrak{D}}$, where for ${\mathcal{D}}\in\mathfrak{D}$, whose vertices are $x_{K}$, $x_{L}$, $x_{{K}^{*}}$, $x_{{L}^{*}}$,

[TABLE]

We conclude this section with a remark on the particular structure of the scalar product of two discrete gradients $(\nabla^{\mathfrak{D}}u_{\mathcal{T}},\nabla^{\mathfrak{D}}v_{\mathcal{T}})_{{\boldsymbol{\Lambda}},\mathfrak{D}}$ for $u_{\mathcal{T}},v_{\mathcal{T}}\in\mathbb{R}^{\mathcal{T}}$. Indeed, for $u_{\mathcal{T}}\in\mathbb{R}^{\mathcal{T}}$ and ${\mathcal{D}}\in\mathfrak{D}$, we define $\delta^{\mathcal{D}}u_{\mathcal{T}}$ by

[TABLE]

Then, we can write

[TABLE]

where the local matrices ${\mathbb{A}}^{{\mathcal{D}}}$ are defined by

[TABLE]

It follows from elementary calculations left to the reader that the condition number of ${\mathbb{A}}^{\mathcal{D}}$ with respect to the 2-norm can be bounded by

[TABLE]

2.3. The nonlinear DDFV scheme

Let $N_{T}$ be a positive integer, we consider for simplicity the constant time step is given by $\Delta t=T/{N_{T}}$. For $n\in\{0,\dots,N_{T}\}$, we denote by $t^{n}=n\Delta t$. We first discretize the initial condition by taking the mean values of $u_{0}$, i.e.,

[TABLE]

and the exterior potential $V$ by taking its nodal values on the primal and dual cells, i.e.,

[TABLE]

It defines in particular $u_{\mathcal{T}}^{0}$ and $V_{\mathcal{T}}$.

The scheme requires a stabilization parameter denoted by $\kappa>0$. It is a fixed parameter. Then, for all $n\geq 0$, we look for $u_{\mathcal{T}}^{n+1}\in(\mathbb{R}_{+}^{\ast})^{\mathcal{T}}$ solution to the following variational formulation:

[TABLE]

Let us mention that, in view of its implementation, the scheme can be rewritten on each mesh as follows:

[TABLE]

2.4. Functional spaces

For a given vector $u_{\mathcal{T}}$ defined on a DDFV mesh $\mathcal{T}$ of size $h$, one usually reconstructs three different approximate solutions : $u_{h,\mathfrak{M}}$ is a piecewise constant reconstruction on the primal mesh, $u_{h,\overline{\mathfrak{M}^{*}}}$ is a piecewise constant reconstruction on the dual mesh and $u_{h}$ is the mean value of $u_{h,\mathfrak{M}}$ and $u_{h,\overline{\mathfrak{M}^{*}}}$. They are defined by

[TABLE]

Then, the set of the approximate solutions is denoted by $H_{\mathcal{T}}$:

[TABLE]

In the sequel, we will also need some reconstruction of the approximate solutions on the diamond cells. Thanks to the reconstruction operator on diamonds $r^{\mathfrak{D}}$, we can define $u_{\mathcal{D}}=r^{\mathcal{D}}(u_{\mathcal{T}})$ for all ${\mathcal{D}}\in\mathfrak{D}$ for instance. Therefore, we can define a piecewise constant function on diamond cells $u_{h,\mathfrak{D}}$ by $u_{h,\mathfrak{D}}=\operatorname*{\sum}_{{\mathcal{D}}\in\mathfrak{D}}u_{\mathcal{D}}\boldsymbol{1}_{\mathcal{D}}$. The set of such functions is denoted $H_{\mathfrak{D}}$.

For a function $u_{h}\in H_{\mathcal{T}}$, we define its approximate gradient $\nabla^{h}u_{h}\in(H_{\mathfrak{D}})^{2}$ by

[TABLE]

As the problem (1) is an evolutive problem, the numerical scheme (20) defines $u_{\mathcal{T}}^{n}\in\mathbb{R}^{\mathcal{T}}$ for all $n\in\{0,\ldots,N_{T}\}$. We consider approximate solutions which are piecewise constant in time. Therefore, we define the space-time approximation spaces $H_{\mathcal{T},{\Delta t}}$ and $H_{\mathfrak{D},{\Delta t}}$ based respectively on $H_{\mathcal{T}}$ and $H_{\mathfrak{D}}$:

[TABLE]

We still keep the notation $\nabla^{h}$ to define the approximate gradient of $u_{h,{\Delta t}}\in H_{\mathcal{T},{\Delta t}}$:

[TABLE]

Therefore, for all $u_{h,{\Delta t}}\in H_{\mathcal{T},{\Delta t}}$, we have $\nabla^{h}u_{h,{\Delta t}}\in(H_{\mathfrak{D},{\Delta t}})^{2}$. Furthermore, we introduce the following reconstructions

[TABLE]

We may now introduce some norms on the functional spaces $H_{\mathcal{T}}$ and $H_{\mathcal{T},{\Delta t}}$. For a discrete solution $u_{\mathcal{T}}\in\mathbb{R}^{\mathcal{T}}$, we define $\left|u_{\mathcal{T}}\right|_{p,\mathcal{T}}$ for $1\leq p\leq\infty$ by

[TABLE]

It permits to define discrete $W^{1,p}$-norms ($1\leq p\leq+\infty$) and a discrete $W^{-1,1}$-norm on $H_{\mathcal{T}}$. For all $u_{h}\in H_{\mathcal{T}}$, we set

[TABLE]

Let us just remark that, as $\nabla^{h}u_{h}$ is a piecewise constant function on diamonds, we have:

[TABLE]

Then, we define some discrete $L^{q}(0,T;W^{1,p}(\Omega))$ ($1\leq p,q<+\infty$), $L^{\infty}(0,T;W^{1,\infty}(\Omega))$ and $L^{\infty}(0,T;L^{p}(\Omega))$-norms on $H_{\mathcal{T},{\Delta t}}$. For all $u_{h,{\Delta t}}\in H_{\mathcal{T},{\Delta t}}$, we set:

[TABLE]

2.5. Main results

The numerical analysis of the scheme strongly relies on a discrete version of the energy/energy dissipation relation (4). In order to make it explicit, let us introduce the discrete counterpart $({\mathbb{E}}_{\mathcal{T}}^{n})_{n\geq 0}$ of the free energy ${\mathbb{E}}$ defined by (5):

[TABLE]

and the discrete counterpart $({\mathbb{I}}^{n}_{\mathcal{T}})_{n\geq 1}$ of the dissipation ${\mathbb{I}}$ defined by (6):

[TABLE]

The first main result of our paper is the existence of a positive solution to the nonlinear scheme (20); it is stated in Theorem 2.1. The mesh is given and fulfills the very permissive requirements of Section 2.1. Our nonlinear scheme (20) yields a nonlinear system of algebraic equations. The fact that this system admits a solution is not obvious and is ensured by Theorem 2.1. The proof strongly relies on the fact that the scheme fulfills a discrete entropy/dissipation relation.

Theorem 2.1** (Existence of a discrete solution).**

For all $n\geq 0$, there exists a solution $u_{\mathcal{T}}^{n+1}\in\left(\mathbb{R}_{+}^{\ast}\right)^{\mathcal{T}}$ to the nonlinear system (20) that satisfies the discrete entropy/entropy dissipation estimate

[TABLE]

Once the existence of $u_{\mathcal{T}}^{n+1}$ at hand for all $n\geq 0$, we can reconstruct the approximate solutions $u_{h,{\Delta t}}$, $u_{h,{\Delta t},\mathfrak{M}}$, and $u_{h,{\Delta t},\mathfrak{M}^{*}}$. The convergence of these approximate solutions towards a weak solution when the mesh size and the time step tend to 0 is then a very natural question. This question is addressed in Theorem 2.2.

In what follows, $\left(\mathcal{T}_{m}\right)_{m\geq 1}=\left(\overline{\mathfrak{M}}_{m},\overline{\mathfrak{M}_{m}^{*}}\right)_{m\geq 1}$ denotes a sequence of admissible discretization of $\Omega$ and $\left(\mathfrak{D}_{m}\right)_{m\geq 1}$ denotes the corresponding diamond mesh. We assume that the

[TABLE]

the regularity factors $\theta_{\mathcal{D}}$ and $\widetilde{\theta}_{\mathcal{D}}$ being defined by (10).

Concerning the time discretization, we consider a sequence $\left({N_{T,m}}\right)_{m\geq 1}$ of positive integers tending to $+\infty$, and we denote by $\left({\Delta t}_{m}\right)_{m\geq 1}=\left(\frac{T}{N_{T,m}}\right)_{m\geq 1}$ the corresponding sequence of time steps. For technical reasons that will appear later on, and even though this condition does not seem to be mandatory from a practical point of view, we have to make the assumption that there exists some constant $C_{1}>0$ such that

[TABLE]

The existence of a discrete solution $u_{\mathcal{T}_{m}}^{n+1}$ to the scheme (20) for all $n\in\{0,\dots,N_{T,m}-1\}$ and all $m\geq 1$ stated in Theorem 2.1 allows us to define the approximate solutions $u_{h_{m},{\Delta t}_{m}}$, $u_{h_{m},{\Delta t}_{m},\mathfrak{M}_{m}}$, and $u_{h_{m},{\Delta t}_{m},\overline{\mathfrak{M}^{*}_{m}}}$ for all $m\geq 1$. The next theorem ensures that, up to a subsequence, the sequences of approximate solution converge towards a weak solution of the problem (1).

Theorem 2.2** (Convergence towards a weak solution).**

Assume that (26) and (27) holds. Then there exists a weak solution $u$ in the sense of Definition 1.1 such that, up to a subsequence,

[TABLE]

for all $p\in[1,\infty)$.

Further convergence properties are established during the proof of Theorem 2.2. We don’t make them explicit here in order to minimize the notations and to improve the readability of the paper. We refer to Section 4.1 for refined statements.

Remark 2.3**.**

The convergence of the scheme is only assessed up to a subsequence in Theorem 2.2. This comes from the fact that the uniqueness of weak solutions in the sense of Definition 1.1 is still an open problem even for initial data $u_{0}$ belonging to $L^{2}(\Omega)$. However, we conjecture that for $u_{0}$ being such that $H(u_{0})$ belongs to $L^{1}(\Omega)$, the weak solutions in the sense of Definition 1.1 are renormalized solutions (see for instance [9]). Uniqueness should follow, implying the convergence of the whole sequence.

The main goal of the paper is to prove Theorems 2.1 and 2.2. The proof is articulated as follows: In Section 3, we derive some estimates on the discrete solution. These a priori estimates allow us to show that the nonlinear system originating from the (20) admits (at least) one solution, as claimed in Theorem 2.1. Most of the estimates derived in Section 3 are uniform w.r.t. $m$. This provides enough compactness on the approximate solutions to pass to the limit $m\to\infty$ in Section 4.

3. Energy and dissipation estimates, existence of a solution to the scheme

3.1. a priori estimates

The first statement of this section is devoted to what we call the fundamental estimates, that are discrete counterparts of the conservation of mass and of the energy/dissipation relation (4). All the further a priori estimates on the discrete solution are based on these two estimates.

Proposition 3.1** (fundamental estimates).**

Let $\left(u_{\mathcal{T}}^{n}\right)_{n\geq 1}$, with $u_{\mathcal{T}}^{n}\in\left(\mathbb{R}_{+}^{*}\right)^{\mathcal{T}}$ for all $n\geq 0$, be a solution to the scheme (20) corresponding to the initial data $u_{0}$. Then,

1. (i)

the mass is conserved along time, i.e.,

[TABLE] 2. (ii)

the discrete free energy is dissipated along time, i.e.,

[TABLE]

Moreover, the discrete free energy is decaying along time and is bounded:

[TABLE]

and the “integrated over time” dissipation is also bounded:

[TABLE]

Proof.

Equation (28) is obtained directly by choosing $\psi_{\mathcal{T}}=1_{\mathcal{T}}$ in (20a). In order to get Estimate (29), it suffices to take $\psi_{\mathcal{T}}=g_{\mathcal{T}}^{n+1}$ in (20a) and to remark that, because of the convexity of $u\mapsto H(u)+uV$, one has $\llbracket u_{\mathcal{T}}^{n+1}-u_{\mathcal{T}}^{n},g_{\mathcal{T}}^{n+1}\rrbracket_{\mathcal{T}}\geq{\mathbb{E}}_{\mathcal{T}}^{n+1}-{\mathbb{E}}^{n}_{\mathcal{T}}.$ Inequality (30) is just a consequence of (29), of the nonnegativity of the dissipation and of the penalization term and of Jensen’s inequality. By summation of (29) over $n$, we deduce (31). ∎

The goal of the remaining part of this section is to take advantage of the fundamental estimates of Proposition 3.1 to derive some further estimates to be used in the numerical analysis. As in the continuous framework, the energy/dissipation estimate (29) is used in order to estimate the discrete counterpart of the Fisher information. But in the discrete framework, the chain rule appearing in (7) does not longer hold, and we have to manipulate several objects related to the Fisher information. The last goal of this section is to prove that, for a fixed grid and a fixed time step, the discrete solutions $\left(u_{\mathcal{T}}^{n}\right)_{n\geq 1}$ is uniformly bounded away from 0, cf. Lemma 3.5.

Define the discrete fields $\xi_{\mathcal{T}}^{n}=\sqrt{u_{\mathcal{T}}^{n}}$ which play a key role as in the continuous level. In order to relate different discrete counterparts of the Fisher information, we first have to derive some properties on the local diffusion matrices ${\mathbb{A}}^{\mathcal{D}}$, ${\mathcal{D}}\in\mathfrak{D}$.

Let ${\mathcal{D}}\in\mathfrak{D}$, then we define the diagonal matrix ${\mathbb{B}}^{\mathcal{D}}$ by

[TABLE]

For all ${\boldsymbol{w}}=(w_{\sigma},w_{\sigma^{*}})^{T}\in\mathbb{R}^{2}$ and all ${\mathcal{D}}\in\mathfrak{D}$, there holds

[TABLE]

where $\|\cdot\|_{q}$ is the usual matrix $q$-norm and $|\cdot|$ is the Euclidian norm on $\mathbb{R}^{2}$. It follows from the equivalence of the matrix $1$- and $2$- norms on the finite dimensional space $\mathbb{R}^{2\times 2}$ that

[TABLE]

for some $\gamma\geq 1$. Therefore, in view of (17), we get the existence of $C_{2}$ depending only on $\theta^{\star}$, $\lambda_{m}$ and $\lambda^{M}$ such that

[TABLE]

We introduce now a discrete counterpart of $\int_{\Omega}u\Lambda\nabla\log u\cdot\nabla\log u$, it is the quantity $\widehat{\mathbb{I}}_{\mathcal{T}}^{n}$ defined by

[TABLE]

Let us first relate this quantity to a discrete Fisher information.

Lemma 3.2**.**

For all $n\geq 0$, there holds

[TABLE]

Proof.

Thanks to the first inequality of (33), one has

[TABLE]

It results from the elementary inequality

[TABLE]

that for all ${\mathcal{D}}\in\mathfrak{D}$, one has

[TABLE]

This yields

[TABLE]

and since $\max({u_{K}^{n+1}},{u_{L}^{n+1}},{u_{K^{*}}^{n+1}},{u_{L^{*}}^{n+1}})\leq 4r^{\mathcal{D}}(u_{\mathcal{T}}^{n+1})$, one gets

[TABLE]

In order to conclude the proof of Lemma 3.2, it only remains to combine (35) and (36). ∎

We now want to get a bound on $\widehat{\mathbb{I}}^{n+1}_{\mathcal{T}}$ in order to deduce some bound on $\|\nabla^{\mathfrak{D}}\xi_{\mathcal{T}}^{n+1}\|_{{\mathbf{\Lambda}},\mathfrak{D}}$. Therefore, we first need to establish an estimate on the discrete reconstruction by diamond $r^{\mathfrak{D}}[u_{\mathcal{T}}^{n+1}]$.

Lemma 3.3**.**

Let $r^{\mathfrak{D}}[u_{\mathcal{T}}^{n+1}]\in\mathbb{R}^{\mathfrak{D}}$ be defined by (15). There exists $C>0$, depending only on $\Omega$, $\lambda_{m}$, $\lambda^{M}$ and $\theta^{\star}$ such that

[TABLE]

Proof.

The definition (15) implies that

[TABLE]

where we have set

[TABLE]

The terms $T^{n+1}$ and $T^{{n+1,*}}$ can be rewritten

[TABLE]

where ${{\mathcal{D}}_{L}}$ denotes the unique diamond cell associated to the primal boundary cell $L\in{\partial}\mathfrak{M}$. The terms $T_{\mathfrak{M}}^{n+1}$ and $T^{{n+1,*}}$ can be estimated thanks to the regularity of the mesh (13) by

[TABLE]

We deduce from (28) that

[TABLE]

Let us now focus on the term $T_{{\partial}\mathfrak{M}}^{n+1}$. The area of the diamond cell ${\mathcal{D}}$ corresponding to a boundary edge $\sigma\subset{\partial}\Omega$ (or equivalently to a primal boundary cell $L\in{\partial}\mathfrak{M}$) can be estimated by

[TABLE]

Therefore, we get that

[TABLE]

where $\gamma_{{\partial}\mathfrak{M}}u_{\mathcal{T}}({\boldsymbol{x}})=\operatorname*{\sum}_{{{L}}\in\partial\mathfrak{M}}u_{{{L}}}\boldsymbol{1}_{{{L}}}({\boldsymbol{x}}),$ $\forall{\boldsymbol{x}}\in{\partial}\Omega.$ The trace inequality stated in Theorem A.1 gives with $v_{\mathcal{T}}=\xi_{\mathcal{T}}^{n+1}$

[TABLE]

Thanks to (3) and the regularity of the mesh (10) and (11), there exists $C>0$ depending only on $\lambda_{m}$, $\lambda^{M}$ and $\theta^{\star}$ such that

[TABLE]

Since $|\xi_{\mathcal{T}}^{n+1}|_{2,\mathcal{T}}^{2}=|u_{\mathcal{T}}^{n+1}|_{1,\mathcal{T}}$, Proposition 3.1 and Lemma 3.2 provide that

[TABLE]

∎

Thanks to (37), it is now possible to relate $\widehat{\mathbb{I}}^{n+1}_{\mathcal{T}}$ to ${\mathbb{I}}^{n+1}_{\mathcal{T}}$.

Lemma 3.4**.**

There exist $C>0$ and $h^{\star}>0$ depending only on $u_{0}$, $V$, $\lambda_{m}$, $\lambda^{M}$ and $\theta^{\star}$ such that

[TABLE]

Proof.

Bearing in mind the definition (24) of the dissipation ${\mathbb{I}}_{\mathcal{T}}^{n+1}$, we deduce thanks to (33) that

[TABLE]

But, as $g_{\mathcal{T}}^{n+1}=\log u_{\mathcal{T}}^{n+1}+V_{\mathcal{T}}$, the elementary inequality $(a+b)^{2}\leq 2(a^{2}+b^{2})$ implies that

[TABLE]

It follows from the regularity of $V$ and from the regularity of the mesh (11) that

[TABLE]

for some $C$ depending only on $\lambda^{M}$ and $\theta^{\star}$. Therefore, we deduce from Lemma 3.3 that

[TABLE]

We infer from (41) and (42) that if ${\rm size}(\mathcal{T})$ is small enough, (40) holds. ∎

We have at hand the necessary tool to address the uniform positivity of the solutions. The next lemma states that the discrete solutions remain bounded away from [math] by a small quantity $\epsilon>0$ depending on the data of the continuous problem and on the discretization parameters. This information is of great importance since, because of the singularity of the $\log$, the nonlinear functional corresponding to the scheme is not continuous on the boundary of $\left(\mathbb{R}_{+}^{*}\right)^{\mathcal{T}}$. The proof is inspired from the ones of [13, Lemma 3.10] and [14, Lemma 3.7]. We sketch it here to highlight how we overpass the difficulties related to the fact that there is a limited communication between the primal and dual meshes.

Lemma 3.5**.**

There exists $\epsilon>0$ depending on the data $u_{0}$ and $V$, on the mesh $\mathcal{T}$, on the time step ${\Delta t}$, and on the stabilization parameter $\kappa$, such that

[TABLE]

Proof.

In this proof as elsewhere in the paper, the generic constants $C$ only depend on the data of the continuous problem and on the regularity bound $\theta^{\star}$ for the mesh. In order to highlight the dependency of a quantity with respect to the mesh or to the time step, we use subscripts. For instance $C_{\Delta t}$ may depend on ${\Delta t}$, whereas $C_{\mathcal{T}}$ may depend on the mesh and $C_{\mathcal{T},{\Delta t}}$ may depend on the mesh and on the time step.

Owing to (28), there exists $M_{0}\in\mathfrak{M}\cap\overline{\mathfrak{M}^{*}}$ such that

[TABLE]

implying in particular that

[TABLE]

On the other hand, it follows from (31) that ${\mathbb{I}}_{\mathcal{T}}^{n+1} \leq C_{\Delta t}$. Together with Estimate (40), this provides that

[TABLE]

Assume for instance that $M_{0}=K_{0}\in\mathfrak{M}$ (the case $M_{0}\in\overline{\mathfrak{M}^{*}}$ is similar). Since $B^{\mathcal{D}}_{\sigma}>\frac{1}{C}$ for all ${\mathcal{D}}\in\mathfrak{D}$, we deduce from (44)–(45) that $\log(u_{K_{1}}^{n+1})\geq-C_{\Delta t}$ for all neighboring cell ${K_{1}}\in\mathfrak{M}$ such that $K_{0}|K_{1}\in\mathcal{E}$ and for all $L\in{\partial}\mathfrak{M}\cap{\partial}K_{0}$. It follows from a simple induction based on the reproduction of this argument (see [13, Lemma 3.10] or [14, Lemma 3.7] for details) that

[TABLE]

Let $K^{*}\in\overline{\mathfrak{M}^{*}}$, then there exists $K\in\mathfrak{M}$ such that m$({{K}}\cap{{K}^{*}})>0$. Thanks to the penalization term in Estimate (29), one has that

[TABLE]

The regularity of $V$ implies that

[TABLE]

hence, using (46), one gets that

[TABLE]

The relation (43) follows from (46) and (47). ∎

3.2. Existence of a solution to the scheme

The numerical scheme (20) amounts at each time step $n\geq 0$ to solve a nonlinear system $\boldsymbol{\mathcal{F}}^{n}(u_{\mathcal{T}}^{n+1})={\bf 0}$. The existence of a solution $u_{\mathcal{T}}^{n+1}$ is therefore non trivial. It is established in the following proposition.

Proposition 3.6**.**

For all $n\geq 0$, there exists (at least) one solution $u_{\mathcal{T}}^{n+1}\in(\mathbb{R}_{+}^{\ast})^{\mathcal{T}}$ to the nonlinear system (20).

The proof of Proposition 3.6 relies on a topological degree argument [39, 21, 27]. The key point is that, owing to Lemma 3.5, one can restrain our search for the solution on a compact subset of $\left(\mathbb{R}_{+}^{*}\right)^{\mathcal{T}}$ on which the functional $\boldsymbol{\mathcal{F}}^{n}$ is (uniformly) continuous, making Leray-Schauder’s theorem applicable. We do not detail the proof here since it is very close to the one of [14, Proposition 3.8].

4. Convergence w.r.t. discretization parameters

4.1. Compactness of the approximate solutions

Thanks to Proposition 3.6, we have at hand discrete solutions $\left(u_{\mathcal{T}}^{n+1}\right)_{n\geq 0}$ corresponding to all the time steps, and thus the corresponding reconstructions $u_{h,\Delta t}\in H_{\mathcal{T},{\Delta t}}$, $u_{h,\mathfrak{D},\Delta t}$ as defined in Section 2.4. We also define $\xi_{h,\Delta t}\in H_{\mathcal{T},{\Delta t}}$ based on $\xi_{\mathcal{T}}^{n}=\sqrt{u_{\mathcal{T}}^{n}}$ for all $n\in\{0,\dots,N_{T}\}$.

Thanks to the estimates established in Section 3.1, we can obtain some further estimates satisfied by the discrete reconstructions. These estimates, stated in Lemma 4.1, will then be used to deduce some compactness properties of sequences of approximate solutions.

Lemma 4.1**.**

1. (i)

There exists $C$ depending only on $u_{0}$ and $V$ such that

[TABLE] 2. (ii)

There exists $C$ depending only on $u_{0}$, $V$, $\lambda_{m}$, $\lambda^{M}$ and $\theta^{\star}$ such that, for ${\rm size}(\mathcal{T})$ small enough, one has

[TABLE] 3. (iii)

There exists $C$ depending only on $\theta^{\star}$, $u_{0}$, $V$, $C_{\ref{CFL}}$, $\lambda_{m}$ and $\lambda^{M}$ such that

[TABLE]

Proof.

The first inequality in (48) is just a consequence of Jensen’s inequality because $H$ is a convex function. Moreover, since we assumed that $V\geq 0$ and proved the positivity of the discrete solution, we have:

[TABLE]

The last inequality in (48) is then a straightforward consequence of (30). The estimates in (49) are deduced from (31), (40) and Lemma 3.2.

It remains to prove (50). From Lemma 3.3, we have that

[TABLE]

We infer from the first inequality of (49) that $\widehat{\mathbb{I}}_{\mathcal{T}}^{n}\leq{C(1+T)}/{{\Delta t}}$ and the assumption (27) implies that

[TABLE]

This concludes the proof of Lemma 4.1.

∎

In order to get the compactness of a sequence of approximate solutions, it is also crucial to establish a discrete counterpart of a $L^{1}(0,T;W^{-1,1}(\Omega))$ estimate on the discrete time derivative. In what follows, we denote by

[TABLE]

Lemma 4.2**.**

There exists $C$ depending only on $V$, $u_{0}$, $T$, $\kappa$, $\theta^{\star}$, $C_{\ref{CFL}}$, $\lambda_{m}$ and $\lambda^{M}$ such that

[TABLE]

Proof.

We proceed as in [18, Lemma 3.4]. It follows from (20a) that for all $\psi_{\mathcal{T}}\in\mathbb{R}^{\mathcal{T}}$, one has

[TABLE]

The application $(g_{\mathcal{T}}^{n+1},\psi_{\mathcal{T}})\mapsto T_{\mathfrak{D}}(u_{\mathcal{T}}^{n+1},g_{\mathcal{T}}^{n+1},\psi_{\mathcal{T}})+\kappa\llbracket\mathcal{P}(g_{\mathcal{T}}^{n+1}),\psi_{\mathcal{T}}\rrbracket_{\mathcal{T}}$ is a scalar product on $\mathbb{R}^{\mathcal{T}}$, then it follows from Cauchy-Schwarz inequality that

[TABLE]

We can estimate the term $T_{\mathfrak{D}}(u_{\mathcal{T}}^{n+1},\psi_{\mathcal{T}},\psi_{\mathcal{T}})$ by

[TABLE]

Thanks to (50) and since $u_{h,\mathfrak{D}}^{n+1}\geq 0$, one gets that

[TABLE]

therefore

[TABLE]

Since $\psi_{\mathcal{T}}$ can be chosen arbitrarily, this implies that

[TABLE]

Thanks to Cauchy-Schwarz inequality once again, we obtain that

[TABLE]

One concludes the proof by using (31). ∎

Let $\left(\mathcal{T}_{m}\right)_{m\geq 1}$ be a sequence of meshes as in Section 2.1 such that ${\rm size}(\mathcal{T}_{m})$ tends to [math] as $m$ tends to $\infty$, and such that the regularity of the discretization $\mathcal{T}_{m}$ is uniformly bounded w.r.t. $m$, i.e.,

[TABLE]

Assume moreover that simultaneously the time step ${\Delta t}_{m}$ tends to [math] as $m$ tends to $\infty$ while satisfying (27). The estimates stated in this section are uniform w.r.t. $m$. We deduce from Lemma 4.1 that the sequences $\left(u_{h_{m},{\Delta t}_{m}}\right)_{m\geq 1}$, $\left(u_{h_{m},{\Delta t}_{m},\mathfrak{M}}\right)_{m\geq 1}$, and $\left(u_{h_{m},{\Delta t}_{m},\mathfrak{M}^{*}}\right)_{m\geq 1}$ are equi-integrable in $L^{1}(Q_{T})$ while uniformly bounded in $L^{\infty}(0,T;L^{1}(\Omega))$, hence there exists $u,u^{(1)},u^{(2)}\in L^{\infty}(0,T;L^{1}(\Omega))$ such that, for all $p\in[1,\infty)$, the following convergence holds up to the extraction of an unlabeled subsequence:

[TABLE]

Moreover, it follows from (28) and (49) that

[TABLE]

Thus, similarly to [2, Proposition 5.3] (which is strongly related to [3, Lemma 3.8]), we get the existence of $\xi\in L^{2}(0,T;H^{1}(\Omega))$ such that

[TABLE]

Up to now, we only have weak convergence results, which are not sufficient to pass to the limit in the scheme because of the nonlinearities. The purpose of the following proposition is to recover some strong convergence results. Our approach is based on the time-compactness toolbox presented in [4] (see also [33] for a closely related approach).

Proposition 4.3**.**

Up to the extraction of an unlabeled subsequence,

[TABLE]

for all $p\in[1,\infty)$. Moreover, let $\xi\in L^{2}(0,T,H^{1}(\Omega))$ be as in (53), then $\xi=\sqrt{u}$.

Proof.

The proof is divided into three steps. We remove the subscripts $m$ for the ease of reading.

Step 1. The goal of this part is to make use on both $\left(u_{h,{\Delta t},\mathfrak{M}}\right)$ and $\left(u_{h,{\Delta t},\mathfrak{M}^{*}}\right)$ of the time-compactness criterion of [4, Theorem 3.9].

It follows directly from their definitions that

[TABLE]

Thanks to discrete Poincaré-Sobolev Inequality [8], it follows from (52) that

[TABLE]

for some $C_{p}$ depending only on $p$, on $\Omega$ and on the regularity $\theta^{\star}$ of the mesh and therefore

[TABLE]

On the other hand, the mass conservation (28) and the positivity of the solutions yield

[TABLE]

It results from (55), (56) and from Riesz-Thorin interpolation theorem that

[TABLE]

thus in particular for $p=3/2$. The weak limits $u^{(1)}$, $u^{(2)}$ of $u_{h,{\Delta t},\mathfrak{M}}$ and $u_{h,{\Delta t},\mathfrak{M}^{*}}$ thus belong to $L^{3/2}(Q_{T})$, and the weak limits of $\xi_{h,{\Delta t},\mathfrak{M}}$ and $\xi_{h,{\Delta t},\mathfrak{M}^{*}}$ (up to the extraction of an unlabeled subsequence), denoted by $\xi^{(1)}$, $\xi^{(2)}$, belong to $L^{3}(Q_{T})$ of . The functions $u_{h,{\Delta t},\mathfrak{M}}$ and $\xi_{h,{\Delta t},\mathfrak{M}}$, as well as $u_{h,{\Delta t},\mathfrak{M}^{*}}$ and $\xi_{h,{\Delta t},\mathfrak{M}^{*}}$, are thus in duality.

We now want to apply [4, Theorem 3.9] in order to show that $\xi^{(1)}=\sqrt{u^{(1)}}$ and $\xi^{(2)}=\sqrt{u^{(2)}}$. Therefore, we have to verify that the three assumptions $\bf(A_{x}1)$, $\bf(A_{x}2)$ and $\bf(A_{x}3)$ of [4] are satisfied.

1. (i)

As a direct consequence of the arguments developed in the proofs of [3, Lemma 3.8] or [18, Proposition 4.3] and of Poincaré Sobolev embedding [8], any sequence $\left(v_{\mathcal{T}_{m}}\right)_{m}$ such that $\left\|v_{h}\right\|_{1,2,\mathcal{T}}\leq C$ is such that $v_{h,\mathfrak{M}}$ and $v_{h,\mathfrak{M}^{*}}$ converges strongly in $L^{3}(\Omega)$ (up to the extraction of an unlabeled subsequence). Therefore Assumption $\bf(A_{x}1)$ of [4] is fulfilled. 2. (ii)

The reconstructions $v_{h,\mathfrak{M}}$ and $v_{h,\mathfrak{M}^{*}}$ from $v_{\mathcal{T}}$ are piecewise constant, the functions being equal to nodal values of $v_{\mathcal{T}}$ a.e. in $\Omega$, then Assumption $\bf(A_{x}2)$ of [4] is fulfilled by these reconstructions. Note here that this is not the case of the reconstruction $v_{h}$. 3. (iii)

Let $\varphi\in C^{\infty}(\overline{\Omega})$ and let us define $\varphi_{\mathcal{T}}$ by

[TABLE]

Following the proof of [18, Proposition 4.2], there exists $C$ depending only on the regularity of the mesh $\theta^{\star}$ such that

[TABLE]

On the other hand, one can show that

[TABLE]

for some $C$ depending only on the regularity $\theta^{\star}$ of the mesh. Therefore,

[TABLE]

Assumption $\bf(A_{x}3)$ of [4] is thus fulfilled.

We can then make use of Lemma 4.2 and apply [4, Theorem 3.9] to claim that $\xi^{(i)}=\sqrt{u^{(i)}}$ and that

[TABLE]

Setting $u=\left(u^{(1)}+u^{(2)}\right)/2$, we get that

[TABLE]

Moreover, as the sequences $\left(u_{h,{\Delta t}}\right)_{m}$, $\left(u_{h,{\Delta t},\mathfrak{M}}\right)_{m}$, and $\left(u_{h,{\Delta t},\mathfrak{M}^{*}}\right)_{m}$ are uniformly equi-integrable in $L^{p}(0,T,L^{1}(\Omega))$ and converge point-wise, thanks to (62)–(63), we can apply Vitali’s convergence theorem to claim that the sequences converge strongly in $L^{p}(0,T,L^{1}(\Omega))$.

Step 2. Here, the goal is to show that the sequences $\left(u_{h,{\Delta t}}\right)_{m\geq 1}$, $\left(u_{h,{\Delta t},\mathfrak{M}}\right)_{m\geq 1}$ and $\left(u_{h,{\Delta t},\mathfrak{M}^{*}}\right)_{m\geq 1}$ share the same limit $u$, i.e., $u^{(1)}=u^{(2)}=u$. As a consequence of (31), one has

[TABLE]

hence

[TABLE]

The regularity of the exterior potential $V$ implies that

[TABLE]

so that one gets that

[TABLE]

Up to a subsequence, this ensures that $\log(u_{h,{\Delta t},\mathfrak{M}})-\log(u_{h,{\Delta t},\mathfrak{M}^{*}})$ tends to [math] a.e. in $Q_{T}$. But owing to (62), $\log(u_{h,{\Delta t},\mathfrak{M}})$ tends to $\log(u^{(1)})$ while $\log(u_{h,{\Delta t},\mathfrak{M}^{*}})$ tends to $\log(u^{(2)})$. Then $\log(u^{(1)})=\log(u^{(2)})$ a.e. in $Q_{T}$, thus $u^{(1)}=u^{(2)}=u$.

Step 3. In order to conclude the proof of Proposition 4.3, it remains to check that (54d) holds. To this end, we compute $u_{h,{\Delta t},\mathfrak{D}}-u_{h,{\Delta t}}$ on each quarter diamond to get

[TABLE]

Using the identity $|a-b|=\left(\sqrt{a}+\sqrt{b}\right)\left|\sqrt{a}-\sqrt{b}\right|$ for $a,b\geq 0$ , one gets

[TABLE]

Owing to Cauchy-Schwarz inequality, there holds

[TABLE]

It is easy to verify that $B_{\sigma}^{\mathcal{D}}\geq\frac{1}{\theta^{\star}}$ and $B_{\sigma^{*}}^{\mathcal{D}}\geq\frac{1}{\theta^{\star}}$ (see (32) for their definition). Hence, using (50), it provides

[TABLE]

Using (36) together with (49), we obtain that

[TABLE]

thus $u_{h,{\Delta t},\mathfrak{D}}$ also converges towards $u$ in $L^{1}(Q_{T})$. To get the convergence in $L^{p}(0,T,L^{1}(\Omega))$ for any $p\in[1,\infty)$, it only remains to write

[TABLE]

and to use (28) and (50). ∎

The purpose of the following statement is the uniform in time weak-$L^{1}$ in space convergence of the approximate solution towards $u$.

Proposition 4.4**.**

Up to the extraction of an additional subsequence,

[TABLE]

Moreover, the limit function $u$ satisfies

[TABLE]

for some $C$ depending only $u_{0}$ and $V$.

Proof.

Let $R>0$ be arbitrary. As a consequence of the de La Vallée Poussin theorem [22] the space

[TABLE]

is equi-integrable in $L^{1}(\Omega)$, thus it follows from Dunford-Pettis theorem that $E_{R}$ is relatively compact for the weak-$L^{1}(\Omega)$ topology. Since the function $f\mapsto\int_{\Omega}H(f){\rm d}{\boldsymbol{x}}$ is lower semi-continuous, any limit value for a sequence of $E_{R}$ also belongs to $E_{R}$, hence $E_{R}$ is closed, thus compact.

Since $\Omega$ is bounded and because $E_{R}$ is equi-integrable, the $L^{1}(\Omega)$-weak topology coincides with the topology corresponding to the narrow convergence of measures restricted to $E_{R}$. It can thus be endowed with the bounded-Lipschitz metric:

[TABLE]

In the above formula, $f$ and $f_{n}$ ($n\geq 1$) belong to $E_{R}$, and

[TABLE]

We refer for instance to [41, Theorem 5.9] for the equivalence of the topology induced by the bounded-Lipschitz distance with the one of narrow convergence of positive measures. The fact that this latter topology coincides with the weak-$L^{1}(\Omega)$ topology on $E_{R}$ results from its equi-integrability.

As a consequence of Lemma 4.1, there exists $R$ such that $u_{h_{m},{\Delta t}_{m}}(\cdot,t)$ belongs to $E_{R}$ for all $t\in[0,T]$ and all $m\geq 1$. Let $\tau\in(0,T)$, and let $t\in(0,T-\tau)$ and let $\varphi:\overline{\Omega}\to\mathbb{R}$ be Lipschitz continuous with $\|\nabla\varphi\|_{\infty}\leq 1$. Define $\varphi_{\mathcal{T}_{m}}$ as in (58), then (we remove the subscript $m$ for legibility)

[TABLE]

where $N^{(1)}$ and $N^{(2)}$ are the positive integers such that

[TABLE]

Using the scheme (20a), we obtain that

[TABLE]

Cauchy-Schwarz inequality yields

[TABLE]

The first term in the right-hand side is uniformly bounded by ${\mathbb{E}}(0)$ thanks to (31). On the other hand, Estimate (59) together with $\|\nabla\varphi\|_{\infty}\leq 1$ provide that

[TABLE]

Owing to (50), we get that

[TABLE]

Hereby, we deduce that

[TABLE]

and since $\varphi$ was chosen arbitrarily in $\{\varphi:\Omega\to\mathbb{R}\,|\,\|\nabla\varphi\|_{\infty}\leq 1\}$, we obtain that

[TABLE]

We can apply the refined version of Arzelà-Ascoli theorem [1, Proposition 3.3.1] (see also [26, Theorem 4.26]) and claim that $u_{h,{\Delta t}}$ converges uniformly towards $u\in C([0,T);E_{R})$, $E_{R}$ being endowed with the $L^{1}(\Omega)$-weak topology. ∎

4.2. Identification of the limit

The goal of this section is to show that the limit function $u$ exhibited in Proposition 4.3 is a weak solution in the sense of Definition 1.1.

The second and last point to check to complete the proof of Theorem 2.2 is the fact that $u$ is a solution to (1) in the distributional sense, i.e., that the weak formulation (8) is fulfilled. This is the purpose of the following statement.

Proposition 4.5**.**

Let $u$ be as in Proposition 4.3, then $u$ satisfies the weak formulation (8).

Proof.

Here again, we remove the subscript $m$ when it appears us to be detrimental for the readability. Let $\varphi\in C^{\infty}_{c}(\overline{\Omega}\times[0,T))$, then denote by

[TABLE]

Choosing the corresponding $\psi_{\mathcal{T}}=\varphi_{\mathcal{T}}^{n}$ in (20a), multiplying by ${\Delta t}$ and summing over $n\in\{0,\dots,N_{T-1}\}$ leads to

[TABLE]

where

[TABLE]

For the terms $A_{m}$ and $C_{m}$, we can proceed as is [18] to get

[TABLE]

Let us now detail the treatment of the term $B_{m}$ and start by splitting it into

[TABLE]

where

[TABLE]

and

[TABLE]

Since $V$ and $\varphi$ are smooth functions, one has

[TABLE]

whereas ${\boldsymbol{\Lambda}}_{h,\mathfrak{D}}$ converges a.e. towards ${\boldsymbol{\Lambda}}$. Then it follows from (54d) that

[TABLE]

The last term $B_{2,m}$ is treated following the method proposed in [14] that consists in writing

[TABLE]

with

[TABLE]

where we have set

[TABLE]

We know that, up to a subsequence, $\sqrt{u_{h,{\Delta t},\mathfrak{D}}}$ converges strongly in $L^{2}(Q_{T})$ towards $\sqrt{u}$, whereas $\nabla^{h}\xi_{h,{\Delta t}}$ converges weakly in $L^{2}(Q_{T})^{2}$ towards $\nabla\sqrt{u}$, and $\nabla^{h}\varphi_{h,{\Delta t}}(\cdot,\cdot-{\Delta t})$ converges uniformly towards $\nabla\varphi$. Thus we can pass to the limit in $B_{2,m}^{(1)}$ and obtain that

[TABLE]

In order to show that $B_{2,m}^{(2)}$ tends to [math], we need a few preliminaries. Owing to the definition (69) of $\xi_{KL}^{n}$ and $\xi_{K^{*}L^{*}}^{n}$, one always has

[TABLE]

so that, denoting by $\widetilde{\xi}_{h,{\Delta t},\mathfrak{D}}$ and $\widetilde{\xi}^{*}_{h,{\Delta t},\mathfrak{D}}$ the functions defined almost everywhere by

[TABLE]

one obtains that

[TABLE]

for some $C$ depending only on $\theta^{\star}$ and ${\boldsymbol{\Lambda}}$. Similarly, one has

[TABLE]

Bearing in mind that $\nabla^{h}\xi_{h,{\Delta t}}$ is uniformly bounded in $L^{2}(Q_{T})^{2}$ w.r.t. $m$, this ensures in particular that

[TABLE]

As a consequence of (54d), the function $\sqrt{u_{h,{\Delta t},\mathfrak{D}}}$ also converges towards $\xi$ in $L^{2}(Q_{T})$.

We now have at hand all the necessary material to study $B_{2,m}^{(2)}$. It results from Cauchy-Schwarz inequality that

[TABLE]

Thanks to the regularity of the mesh and of $\varphi$, one has

[TABLE]

for some $C>0$ depending only on $\theta^{\star}$ and on $\|\nabla\varphi\|_{\infty}$, whereas Lemma 3.4 and (31) ensure that

[TABLE]

Hence, we get

[TABLE]

Using (71) together with the fact that $\sqrt{u_{h,{\Delta t},\mathfrak{D}}}$ also converges towards $\xi$ in $L^{2}(Q_{T})$, we get that

[TABLE]

We conclude the proof by putting together the statements (66), (67), (68), (70) and (72) in (65). ∎

5. Numerical experiments

5.1. About the practical implementation

The nonlinear system (20) is solved thanks to Newton’s method. In order to avoid the singularity of the $\log$ near [math], the sequence ${(u_{\mathcal{T}}^{n+1,i})}_{i\geq 0}$ to compute $u_{\mathcal{T}}^{n+1}$ from the previous state ${(u_{\mathcal{T}}^{n})}_{i\geq 0}$ is initialized by $u_{\mathcal{T}}^{n+1,0}=\max(u_{\mathcal{T}}^{n},10^{-12})$. In practice, we observe that the threshold criterion is not used. As a stopping criterion, we require the $\ell^{1}$-norm of the residual to be smaller than $10^{-10}$.

5.2. Convergence w.r.t. to the discretization parameters

We test our method on a test case inspired from the one in [14]. We set $\Omega=(0,1)^{2}$, and $V(x_{1},x_{2})=-x_{2}$. The exact solution $u_{\rm ex}$ is then defined by

[TABLE]

with $\alpha=\pi^{2}+\frac{1}{4}$. We choose $u_{0}=u_{\rm ex}(\cdot,0)$. Note that $u_{0}$ vanishes on $\{x_{2}=1\}$.

In order to illustrate the convergence and the robustness of our method, we study its convergence on two sequences of meshes. The first sequence of primal meshes is made of successively refined Kershaw meshes. The second sequence of primal meshes is the so-called quadrangle meshes mesh_quad_i of the FVCA8 benchmark on incompressible flows. One mesh of each sequence is depicted in Figure 2. In the refinement procedure, the time step is divided by $4$ when the mesh size is divided by $2$.

We have introduced a penalization operator in order to prove that reconstruction on the primal mesh $u_{h,\mathfrak{M},{\Delta t}}$ and the reconstruction on the dual mesh $u_{h,\mathfrak{M}^{*},{\Delta t}}$ converge to the same limit. In Table 1, we compute normU the $L^{2}(\Omega\times(0,T))^{2}$ norm of the difference between the two different reconstructions and ordU the corresponding convergence order for different values of $\kappa$ the penalization parameter. We numerically observe the same result : the two reconstructions converge to the same limit even if $\kappa$ is zero.

In the following of this section, the penalization parameter $\kappa$ is set to zero. In Tables 2 and 3, the quantities erru and errgu respectively denote the $L^{\infty}((0,T);L^{2}(\Omega))$ error on the solution and the $L^{2}(\Omega\times(0,T))^{2}$ error on the gradient, whereas ordu and ordgu are the corresponding convergence orders. It appears that the method is slightly more than second order accurate w.r.t. space.

The maximal (resp. mean) number of Newton iterations by time step is denoted by $N_{\text{max}}$ (resp. $N_{\text{mean}}$). We observe that the needed number of Newton iterations starts from a reasonably small value and falls down to $1$ after a small number of time steps. Therefore, our method does not imply an important extra computational cost when compared to linear methods. Eventually, we can check that the minimal value min $u_{\mathcal{T}}^{n}$ remains strictly greater than 0, as proved in Lemma 3.5.

5.3. Long time behavior

In this section, the penalisation parameter $\kappa$ is set to zero. The discrete stationary solution $u_{\mathcal{T}}^{\infty}$ is defined by $u_{K}^{\infty}=\rho e^{-V(x_{K})}$ and $u_{K^{\ast}}^{\infty}=\rho^{\ast}e^{-V(x_{K^{\ast}})}$ for $K \in\overline{\mathfrak{M}}$ and $K^{\ast}\in\overline{\mathfrak{M}}^{\ast}$, the quantities $\rho$ and $\rho^{\ast}$ being fixed so that $\sum_{K\in\mathfrak{M}}u_{K}^{\infty}m_{K}=\sum_{K\in\overline{\mathfrak{M}}^{\ast}}u_{K^{\ast}}^{\infty}m_{K^{\ast}}=\int_{\Omega}u_{0}(x)dx.$ In order to give an evidence of the good large-time behavior of our scheme, we plot in Figure 3 the evolution of the relative energy

[TABLE]

computed on the Kershaw meshes and on the quadrangle meshes. We observe the exponential decay of the relative energy, recovering on general grids the behavior of the Scharfetter-Gummel scheme [19].

Appendix A A trace inequality

First, to a given vector $u_{\mathcal{T}}=\left(\left(u_{{{K}}}\right)_{{\scriptscriptstyle\mathcal{K}}\in{\overline{\mathfrak{M}}}},\left(u_{{{K}^{*}}}\right)_{{\scriptscriptstyle\mathcal{K}^{*}}\in{\overline{\mathfrak{M}}^{*}}}\right)\in\mathbb{R}^{\mathcal{T}}$ defined on a DDFV mesh $\mathcal{T}$, we associate its primal trace $\gamma_{{\partial}\mathfrak{M}}u_{\mathcal{T}}$ on ${\partial}\Omega$ defined by

[TABLE]

Theorem A.1** (Trace inequality).**

Let $\Omega$ be a convex polygonal domain of $\mathbb{R}^{2}$ and $\mathcal{T}$ a DDFV mesh of this domain. There exist $C>0$, depending only on $\Omega$ and $\theta^{\star}$, such that $\forall\ u_{\mathcal{T}}\in\mathbb{R}^{\mathcal{T}}$:

[TABLE]

Proof.

The calculations are similar to those followed in [28, Lemma 10.5] and in [18, Theorem 7.1] for $L^{1}$-norm. The difference comes from the fact that here we define $u_{\partial\mathfrak{M}}$ using the boundary primal mesh instead of the interior primal mesh. Adapting the proof of [18, Theorem 7.1] in the $L^{2}$-norm, we get for ${{K}}\in\mathfrak{M}$ such that $\overline{{{K}}}\cap\partial\Omega\neq 0$ the inequality

[TABLE]

It implies

[TABLE]

Using the fact that $u_{{{L}}}-u_{{{K}}}={\rm m}_{\sigma^{*}}(\nabla^{{\mathcal{D}}}u_{\mathcal{T}})\cdot\boldsymbol{{\boldsymbol{\tau}}}_{{K\!,L}}$, we conclude

[TABLE]

∎