# Développement et analyse de schémas d'ordre élevé pour des modèles de convection-diffusion : étude du comportement en temps long

Julien Moatti

► To cite this version:

Julien Moatti. Développement et analyse de schémas d'ordre élevé pour des modèles de convection-diffusion : étude du comportement en temps long. Analyse numérique [math.NA]. Université de Lille, 2023. Français. NNT : 2023ULILB021 . tel-04415432

HAL Id: tel-04415432

https://theses.hal.science/tel-04415432

Submitted on 24 Jan 2024

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

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

École doctorale MADIS

Unité de recherche Centre Inria de l'Université de Lille

Thèse présentée par **Julien MOATTI**

Soutenue le 26 septembre 2023

En vue de l'obtention du grade de docteur de l'Université de Lille

Discipline **Mathématiques**

# Développement et analyse de schémas d'ordre élevé pour des modèles de convection-diffusion, étude du comportement en temps long

Thèse dirigée par Claire CHAINAIS-HILLAIRET
Maxime HERDA
Simon LEMAIRE

## Composition du jury

Rapporteurs
Lourenço BEIRÃO DA VEIGA professeur à l'Université degli Studi di Milano-Bicocca
Jérôme DRONIOU professeur à Monash University, Melbourne

Examinateurs
Alexandre ERN professeur à l'ENPC - Université Paris Est
Raphaèle HERBIN professeure à l'Université d'Aix-Marseille
Stella KRELL MCF HDR à l'Université Côte d'Azur
Ilaria PERUGIA professeure à l'Universität Wien

président du jury

Invités
Maxime HERDA chargé de recherche chez Inria
Simon LEMAIRE chargé de recherche chez Inria

Directrice de thèse Claire CHAINAIS-HILLAIRET professeure à l'Université de Lille
<figure><img src="image_1.png" /></figure>
UNIVERSITÉ DE LILLE

École doctorale MADIS

Unité de recherche Centre Inria de l'Université de Lille

Thèse présentée par **Julien MOATTI**

Soutenue le 26 septembre 2023

En vue de l'obtention du grade de docteur de l'Université de Lille

Discipline **Mathématiques**

# Développement et analyse de schémas d'ordre élevé pour des modèles de convection-diffusion, étude du comportement en temps long

Thèse dirigée par Claire CHAINAIS-HILLAIRET
Maxime HERDA
Simon LEMAIRE

## Composition du jury

Rapporteurs
Lourenço BEIRÃO DA VEIGA professeur à l'Université degli Studi di Milano-Bicocca
Jérôme DRONIOU professeur à Monash University, Melbourne

Examinateurs
Alexandre ERN professeur à l'ENPC - Université Paris Est
Raphaèle HERBIN professeure à l'Université d'Aix-Marseille
Stella KRELL MCF HDR à l'Université Côte d'Azur
Ilaria PERUGIA professeure à l'Universität Wien

président du jury

Invités
Maxime HERDA chargé de recherche chez Inria
Simon LEMAIRE chargé de recherche chez Inria

Directrice de thèse Claire CHAINAIS-HILLAIRET professeure à l'Université de Lille
<figure><img src="image_2.png" /></figure>
UNIVERSITÉ DE LILLE

Doctoral School MADIS

University Department Centre Inria de l'Université de Lille

Thesis defended by **Julien MOATTI**

Defended on 26th September, 2023

In order to become Doctor from Université de Lille

Academic Field Mathematics

Development and numerical analysis
of high-order schemes for
convection-diffusion models, study of
their long-time behaviour

**Thesis supervised by** Claire CHAINAIS-HILLAIRET
Maxime HERDA
Simon LEMAIRE

**Committee members**

<table>
  <tbody>
    <tr>
      <td>Referees</td>
      <td>Lourenço BEIRÃO DA VEIGA<br>Jérôme DRONIOU</td>
      <td>Professor at Università degli Studi di Milano-Bicocca<br>Professor at Monash University, Melbourne</td>
    </tr>
    <tr>
      <td>Examiners</td>
      <td>Alexandre ERN<br>Raphaèle HERBIN<br>Stella KRELL<br>Ilaria PERUGIA</td>
      <td>Professor at ENPC - Université Paris Est<br>Professor at Université d'Aix-Marseille<br>hdr Associate Professor at Université Côte d'Azur<br>Professor at Universität Wien</td>
    </tr>
    <tr>
      <td>Guests</td>
      <td>Maxime HERDA<br>Simon LEMAIRE</td>
      <td>Junior Researcher at Inria<br>Junior Researcher at Inria</td>
    </tr>
    <tr>
      <td>Supervisor</td>
      <td>Claire CHAINAIS-HILLAIRET</td>
      <td>Professor at Université de Lille</td>
    </tr>
  </tbody>
</table>
<figure><img src="image_3.png" /></figure>
**Mots clés:** analyse numérique, méthode d'entropie, méthodes numériques sur maillages généraux, méthodes numériques d'ordre élevé, semi-conducteurs, edp paraboliques

**Keywords:** numerical analysis, entropy method, numerical methods on general meshes, high-order numerical methods, semiconductors, parabolic pdes
<figure><img src="image_4.png" /></figure>
Cette thèse a été préparée dans les laboratoires suivants.

**Centre Inria de l'Université de Lille**

Parc scientifique de la Haute-Borne
59650 Villeneuve d'Ascq
France

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

**Laboratoire Paul Painlevé**

Cité Scientifique
59655 Villeneuve-d'Ascq
France

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

<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /></figure>
# DÉVELOPPEMENT ET ANALYSE DE SCHÉMAS D'ORDRE ÉLEVÉ POUR DES MODÈLES DE CONVECTION-DIFFUSION, ÉTUDE DU COMPORTEMENT EN TEMPS LONG

## Résumé

Dans cette thèse, nous nous intéressons à l'approximation numérique de problèmes de convection-diffusion, potentiellement anisotropes, par des schémas d'ordre élevé sur maillages généraux. Notre objectif est de proposer des méthodes fiables, précises et efficaces : les solutions numériques doivent préserver les propriétés physiques des solutions continues (conservation de la masse, positivité des densités, comportement en temps long) tout en autorisant une large gamme de paramètres de discrétisation (pas de temps grands, maillages spatiaux généraux) et en optimisant la précision de calcul à coût donné. Les problèmes considérés sont des équations d'advection-diffusion ainsi que des systèmes couplés de dérive-diffusion qui modélisent les composants semi-conducteurs.

On se concentre d'abord sur une équation d'advection-diffusion seule, pour laquelle nous proposons et analysons trois méthodes d'ordre bas de type volumes finis hybrides (HFV). Cette comparaison met en avant la nécessité d'utiliser un schéma non-linéaire afin de préserver la positivité des solutions, tant d'un point de vue théorique que numérique. On s'intéresse alors à l'approximation d'un système de dérive-diffusion, constitué de deux équations d'advection-diffusion couplées avec une équation de Poisson. Pour ce problème, on introduit un schéma non-linéaire basé sur la méthode précédente, qui préserve les bornes des densités calculées (et en particulier leur positivité) et le comportement en temps long de la solution. Ce schéma HFV pour la dérive-diffusion est ensuite comparé numériquement avec un schéma présentant des propriétés similaires basé sur la méthode volumes finis en dualité discrète (DDFV).

Nous nous intéressons alors à des schémas d'ordre élevé (en espace). Ces schémas sont basés sur des méthodes de type hybride d'ordre élevé (HHO) qui peuvent être interprétées comme des extensions à l'ordre arbitraire des méthodes HFV. Nous introduisons deux méthodes pour les équations d'advection-diffusion linéaires. La première est linéaire, tandis que la seconde est non-linéaire, et permet de préserver la positivité. Pour ces deux schémas, nous prouvons l'existence de solutions discrètes et établissons des résultats de comportement en temps long. Nous confirmons également ces résultats numériquement, et mettons en avant la nécessité d'utiliser une méthode non-linéaire pour préserver la positivité. Par ailleurs, on observe que la méthode non-linéaire converge à l'ordre attendu. De plus, la montée en ordre permet un gain d'efficacité (précision de l'approximation pour un coût de calcul donné) conséquent par rapport aux méthodes d'ordre bas des premiers chapitres.

Ces travaux sont complétés par l'étude de problèmes de convection-diffusion avec convection très irrégulière, effectuée en collaboration avec des physiciens durant la thèse. Les recherches menées visent à comprendre comment concevoir des diodes électroluminescentes efficaces émettant dans l'ultraviolet profond, et soulèvent divers enjeux relatifs à la modélisation, l'analyse numérique et la simulation de ces problèmes.

**Mots clés** : analyse numérique, méthode d'entropie, méthodes numériques sur maillages généraux, méthodes numériques d'ordre élevé, semi-conducteurs, edp paraboliques
# DEVELOPMENT AND NUMERICAL ANALYSIS OF HIGH-ORDER SCHEMES FOR CONVECTION-DIFFUSION MODELS, STUDY OF THEIR LONG-TIME BEHAVIOUR

## Abstract

In this thesis, we are interested in the numerical approximation of anisotropic convection-diffusion problems using high-order schemes on general meshes. Our objective is to develop reliable, accurate, and efficient methods: the numerical solutions must preserve the physical properties of the continuous solutions (mass conservation, positivity of densities, long-time behaviour) while allowing for a wide range of discretisation parameters (large time steps, general spatial meshes) and optimising the accuracy at a given computational cost. The problems under study are advection-diffusion equations as well as coupled drift-diffusion systems that model semiconductor devices.

We first focus on a single advection-diffusion equation for which we propose and analyse three different hybrid finite volume (HFV) methods to approximate the solution. Their comparison highlights the necessity of using a nonlinear scheme to ensure the preservation of positivity, both from a theoretical and a numerical level. We then consider the numerical approximation of a drift-diffusion system, consisting of two coupled advection-diffusion equations with a Poisson equation. For this problem, we introduce a nonlinear scheme based on the previous method, which preserves the bounds on the computed densities (including their positivity) and the long-time behaviour of the solution. This HFV scheme for drift-diffusion is then numerically compared with another scheme with similar properties based on the discrete duality finite volume (DDFV) method.

We then focus on high-order (in space) schemes. These schemes are based on hybrid high-order (HHO) methods which can be interpreted as an arbitrary-order extension of HFV methods. We introduce two schemes for linear advection-diffusion: the first method is linear, while the second one is nonlinear and preserves the positivity of the solution. For both of these schemes, we prove the existence of discrete solutions and establish results about their long-time behaviour. We also confirm these results numerically and emphasise the need to use a nonlinear method to preserve positivity. We observe that the nonlinear method converges at the expected order. Furthermore, increasing the order leads to a significant gain in efficiency (accuracy of the approximation for a given computational cost) compared to the low-order methods discussed in the first chapters.

These works are complemented with the study of convection-diffusion problems with highly irregular convection, carried out in collaboration with physicists during the thesis. These investigations aim to understand how to design efficient deep ultraviolet-emitting electroluminescent diodes, and raise various issues related to the modeling, numerical analysis and simulation of these problems.

**Keywords:** numerical analysis, entropy method, numerical methods on general meshes, high-order numerical methods, semiconductors, parabolic pdes
# Sommaire

<table><tr><td>Résumé</td><td>xi</td></tr><tr><td>Sommaire</td><td>xiii</td></tr><tr><td>Liste des tableaux</td><td>xvii</td></tr><tr><td>Table des figures</td><td>xix</td></tr><tr><td>Avant-propos</td><td>1</td></tr><tr><td>Introduction générale</td><td>3</td></tr><tr><td>Phénomènes de convection-diffusion : modélisation, comportement et analyse</td><td>4</td></tr><tr><td>Approximation numérique des problèmes de convection-diffusion</td><td>11</td></tr><tr><td>Overview of the works and results</td><td>18</td></tr><tr><td>1 Long-time behaviour of hybrid finite volume schemes for advection-diffusion equations: linear and nonlinear approaches</td><td>31</td></tr><tr><td>1.1 Introduction</td><td>32</td></tr><tr><td>1.2 Hybrid finite volume discretisation of a variable diffusion problem</td><td>35</td></tr><tr><td>1.3 Definition of the schemes and well-posedness</td><td>40</td></tr><tr><td>1.4 Long-time behaviour</td><td>54</td></tr><tr><td>1.5 Numerical results</td><td>60</td></tr><tr><td>1.6 Conclusion</td><td>67</td></tr><tr><td>1.A Functional inequalities</td><td>67</td></tr><tr><td>1.B Nonlinear scheme for mixed Dirichlet-Neumann boundary conditions</td><td>71</td></tr><tr><td>1.C Proofs of technical results</td><td>73</td></tr><tr><td>2 A structure preserving hybrid finite volume scheme for semiconductor models with magnetic field on general meshes</td><td>77</td></tr><tr><td>2.1 Introduction</td><td>78</td></tr><tr><td>2.2 Discrete setting and schemes</td><td>83</td></tr><tr><td>2.3 Analysis of the stationary scheme</td><td>91</td></tr><tr><td>2.4 Analysis of the transient scheme</td><td>94</td></tr><tr><td>2.5 Numerical results</td><td>101</td></tr><tr><td>2.6 Conclusion</td><td>112</td></tr><tr><td>2.A Discrete boundedness by entropy and dissipation</td><td>113</td></tr></table>
<table><tr><td>3 A comparison of structure-preserving schemes for drift-diffusion systems on general meshes: DDFV vs HFV</td><td>117</td></tr><tr><td>3.1 Motivation</td><td>118</td></tr><tr><td>3.2 Description of the schemes</td><td>119</td></tr><tr><td>3.3 Numerical experiments</td><td>123</td></tr><tr><td>3.4 Conclusion</td><td>126</td></tr></table>

<table><tr><td>4 High-order polytopal schemes for advection-diffusion equations: linear and non-linear approaches</td><td>127</td></tr><tr><td>4.1 Motivations and context</td><td>128</td></tr><tr><td>4.2 Discrete setting and schemes</td><td>131</td></tr><tr><td>4.3 Main features of the schemes</td><td>137</td></tr><tr><td>4.4 Numerical results</td><td>149</td></tr><tr><td>4.5 Conclusion</td><td>160</td></tr></table>

<table><tr><td>Conclusion and perspectives</td><td>161</td></tr><tr><td>Further analysis of the schemes</td><td>161</td></tr><tr><td>More efficient and robust implementations</td><td>162</td></tr><tr><td>Applications to other complex problems</td><td>163</td></tr><tr><td>Numerical comparisons with existing methods</td><td>164</td></tr><tr><td>High-order discretisations in time</td><td>165</td></tr></table>

<table><tr><td>A Semiconductor models with varying band edge energies: an overview</td><td>167</td></tr><tr><td>A.1 Motivations and context</td><td>167</td></tr><tr><td>A.2 A comparison between numerical simulations and experimental results (Appendix B)</td><td>169</td></tr><tr><td>A.3 The question of thermodynamic consistency (Appendix C)</td><td>169</td></tr><tr><td>A.4 Application to the design of efficient LEDs (Appendix D)</td><td>170</td></tr><tr><td>A.5 Some mathematical issues raised by these works</td><td>171</td></tr></table>

<table><tr><td>B Impact of random alloy fluctuations on the carrier distribution in multi-color (In,Ga)N/GaN quantum well systems</td><td>173</td></tr><tr><td>B.1 Introduction</td><td>174</td></tr><tr><td>B.2 Model MQW structures and literature experimental findings</td><td>175</td></tr><tr><td>B.3 Theoretical framework</td><td>178</td></tr><tr><td>B.4 Results</td><td>182</td></tr><tr><td>B.5 Conclusions</td><td>186</td></tr></table>

<table><tr><td>C Importance of satisfying thermodynamic consistency in optoelectronic device simulations for high carrier densities</td><td>189</td></tr><tr><td>C.1 Introduction</td><td>190</td></tr><tr><td>C.2 Drift-diffusion equations and diffusion enhancement</td><td>191</td></tr><tr><td>C.3 Finite volume space discretization</td><td>192</td></tr><tr><td>C.4 Simulations</td><td>194</td></tr><tr><td>C.5 Conclusion</td><td>197</td></tr></table>

<table><tr><td>D Theoretical study of the impact of alloy disorder on carrier transport and recombination processes in deep UV (Al,Ga)N light emitters</td><td>199</td></tr><tr><td>Bibliography</td><td>207</td></tr></table>
<table><tr><td>Sommaire</td><td>xv</td></tr><tr><td>Table des matières</td><td>223</td></tr></table>
Liste des tableaux

<table><tr><td>1.1</td><td>Positivity of discrete solutions. Numerical results for T<sub>f</sub> = 5.10<sup>-4</sup> and Δt = 10<sup>-5</sup> on a tilted hexagonal-dominant mesh. At each time step, there are 4192 cell unknowns and 12512 edge unknowns.</td><td>64</td></tr><tr><td>2.1</td><td>Test-case 2. Comparison of the extremal values and costs for different magnetic fields</td><td>108</td></tr><tr><td>4.1</td><td>Positivity of discrete solutions.</td><td>154</td></tr><tr><td>4.2</td><td>Number of negative cell averages</td><td>155</td></tr><tr><td>B.1</td><td>Material parameters used in the different regions of the simulation supercell. Parameters denoted with † are taken from [194]; parameters denoted with ‡ are derived from [210] as described in the main text.</td><td>182</td></tr></table>
Table des figures

<table>
  <tbody>
    <tr>
      <td>1</td>
      <td>Géométrie d'une diode et maillage adapté associé.</td>
      <td>12</td>
    </tr>
    <tr>
      <td>2</td>
      <td>Long-time behaviour of discrete solutions to HFV schemes on a Kershaw mesh.</td>
      <td>20</td>
    </tr>
    <tr>
      <td>3</td>
      <td>Influence of the magnetic field intensity on the density of electrons of a non-equilibrium steady state.</td>
      <td>22</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Efficiency of non-linear HHO methods.</td>
      <td>25</td>
    </tr>
    <tr>
      <td>1.1</td>
      <td>Two-dimensional discretisation and corresponding notation.</td>
      <td>36</td>
    </tr>
    <tr>
      <td>1.2</td>
      <td><b>Implementation.</b> Coarsest meshes of each family used in the numerical tests.</td>
      <td>60</td>
    </tr>
    <tr>
      <td>1.3</td>
      <td><b>Long-time behaviour of discrete solutions.</b> Comparison of the long-time behaviour on Kershaw meshes for T<sub>f</sub> = 350 and Δt = 0.1.</td>
      <td>63</td>
    </tr>
    <tr>
      <td>1.4</td>
      <td><b>Accuracy of stationary solutions.</b> Relative errors in discrete L<sup>2</sup>-norm and H<sup>1</sup>-seminorm for the first test-case on triangular meshes.</td>
      <td>66</td>
    </tr>
    <tr>
      <td>1.5</td>
      <td><b>Accuracy of stationary solutions.</b> Relative errors in discrete L<sup>2</sup>-norm and H<sup>1</sup>-seminorm for the second test-case on Cartesian meshes.</td>
      <td>67</td>
    </tr>
    <tr>
      <td>2.1</td>
      <td>Two-dimensional discretisation and corresponding notations.</td>
      <td>84</td>
    </tr>
    <tr>
      <td>2.2</td>
      <td>The geometry of the PN diode.</td>
      <td>101</td>
    </tr>
    <tr>
      <td>2.3</td>
      <td><b>Implementation.</b> Coarsest meshes of each family used in the numerical tests.</td>
      <td>102</td>
    </tr>
    <tr>
      <td>2.4</td>
      <td><b>Test-case 1.</b> Evolution of the discrete density of holes.</td>
      <td>105</td>
    </tr>
    <tr>
      <td>2.5</td>
      <td><b>Test-case 2</b> (b = 1). Evolution of the discrete density of holes (note that the scale varies from a figure to the other).</td>
      <td>106</td>
    </tr>
    <tr>
      <td>2.6</td>
      <td><b>Test-case 2</b> (b = 1). Evolution of the discrete extremal values, time step and cost</td>
      <td>107</td>
    </tr>
    <tr>
      <td>2.7</td>
      <td><b>Test-case 3.</b> Evolution of the discrete density of holes, on a tilted hexagonal mesh.</td>
      <td>109</td>
    </tr>
    <tr>
      <td>2.8</td>
      <td><b>Test-case 3</b> (b = 1, tilted hexagonal mesh). Evolution of the discrete extremal values, time step and cost.</td>
      <td>110</td>
    </tr>
    <tr>
      <td>2.9</td>
      <td><b>Test-case 4.</b> Influence of the meshsize: entropy and L<sup>2</sup> distance to equilibrium</td>
      <td>111</td>
    </tr>
    <tr>
      <td>2.10</td>
      <td><b>Test-case 5.</b> Influence of the mesh geometry on the long-time behaviour: entropy and L<sup>2</sup> distance to equilibrium.</td>
      <td>111</td>
    </tr>
    <tr>
      <td>2.11</td>
      <td><b>Test-case 6.</b> Long-time behaviour for the Blakemore statistics and influence of the magnetic field</td>
      <td>112</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Definition of the diamonds D<sub>σ,σ*</sub> and related notations.</td>
      <td>120</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>PN diode geometry.</td>
      <td>124</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td><b>Positivity.</b> Evolution of the discrete minimal values, time step and cost</td>
      <td>125</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td><b>Long-time behaviour.</b> Evolution of the discrete relative entropies.</td>
      <td>126</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td><b>Accuracy of transient solutions.</b> Relative errors on triangular meshes.</td>
      <td>156</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td><b>Accuracy vs. computational cost.</b> Relative errors on triangular meshes.</td>
      <td>157</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>4.3</td>
      <td>Accuracy: nlhfv vs. nlhho. Relative errors on triangular meshes.</td>
      <td>158</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>Long-time behaviour of discrete solutions. Comparison of the long-time behaviour of the solutions to the nonlinear HHO schemes on Kershaw meshes for T<sub>f</sub> = 350 and Δt = 0.1.</td>
      <td>159</td>
    </tr>
    <tr>
      <td>4.5</td>
      <td>Long-time behaviour of discrete solutions. Comparison of the long-time behaviour of the solutions to the nonlinear HHO schemes on Kershaw meshes for T<sub>f</sub> = 350 and Δt = 0.1.</td>
      <td>159</td>
    </tr>
    <tr>
      <td>A.1</td>
      <td>Accuracy of schemes. Boltzmann statistics and heterojunctions.</td>
      <td>169</td>
    </tr>
    <tr>
      <td>B.1</td>
      <td>Conduction and valence band edges (black) along with the quasi-Fermi energies for electrons and holes (grey) in an (In,Ga)N/GaN multi-quantum well system described in virtual crystal approximation. The band edge profile and the quasi Fermi levels are shown at a current density of 50 A/cm<sup>2</sup>. The leftmost (In,Ga)N quantum well contains 12.5% indium while the other two (In,Ga)N wells (centre and right) contain 10% indium.</td>
      <td>176</td>
    </tr>
    <tr>
      <td>B.2</td>
      <td>Schematic illustration of multi-quantum well system. The n-doped region is shown in cyan, the p-doped is in red and undoped regions are in grey. The quantum wells are numbered starting from the n-side.</td>
      <td>177</td>
    </tr>
    <tr>
      <td>B.3</td>
      <td>Profile of (a) valence band edge energy, (b) conduction band edge energy, and (c) radiative recombination rate in the growth plane (c-plane) of an In<sub>0.1</sub>Ga<sub>0.9</sub>N quantum well; the current density is 50 A/cm<sup>2</sup> in all depicted figures. The slice displayed is the through the center well. The data are shown in all cases on a linear scale.</td>
      <td>177</td>
    </tr>
    <tr>
      <td>B.4</td>
      <td>Ratio of radiative recombination ρ, Eq. (B.1), from the shallow wells (In<sub>0.1</sub>Ga<sub>0.9</sub>N) to recombination from the deep well (In<sub>0.1</sub>Ga<sub>0.9</sub>N) calculated as a function of the position of the deep well in the multi-quantum well stack. Here ρ is evaluated using (a) nextnano excluding (purple) and including (green) quantum corrections via a self-consistent Schrödinger-Poisson-drift diffusion solver; results are shown when excluding (solid, filled circles) and including (dotted, open circles) an Al<sub>0.15</sub>Ga<sub>0.85</sub>N blocking layer, and (b) ddfermi excluding (purple), including (green) quantum corrections via localization landscape theory (LLT) using a virtual crystal approximation (VCA) and a random alloy calculation including LLT-based quantum corrections (blue); these calculations neglect the AlGaN blocking layer.</td>
      <td>183</td>
    </tr>
    <tr>
      <td>B.5</td>
      <td>Contribution of each quantum well (n-side; centre; p-side) in the (In,Ga)N multi-quantum well system to the total radiative recombination R<sub>Ω<sub>i</sub></sub><sup>RAD</sup> for i ∈ {n-side, center, p-side} as a percentage of the total radiative recombination from all 3 quantum wells for (a) virtual crystal approximation (VCA), (b) virtual crystal approximation with quantum corrections included via localization landscape theory (VCA + LLT) and (c) a random alloy calculation including localization landscape theory based quantum corrections (Random alloy + LLT). That data are shown as a function of the position of the deep quantum well (x-axis). Each bar contains the percentage recombination from the n-side quantum well (purple), the center quantum well (green) and the p-side quantum well (blue). Labelling is consistent with that introduced in Fig. B.2.</td>
      <td>184</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>B.6</td>
      <td>Hole density (black, solid), electron density (black, dashed), and radiative recombination rate (red, solid) averaged over each atomic plane along the transport direction. Results from calculations building on (i) a virtual crystal approximation (top), (ii) a virtual crystal including quantum corrections via localization landscape theory (LLT) (center) and a (iii) random alloy description including LLT-based quantum corrections (bottom); the deep well is located at (a) the n-side (left), (b) the center (middle) and (c) the p-side (right). The data are shown on a log scale.</td>
      <td>188</td>
    </tr>
    <tr>
      <td>C.1</td>
      <td>Conduction band edge (black) and quasi Fermi energy (red) at a bias of 3.3V (a) when using Boltzmann statistics, (b) when incorrectly using the Scharfetter-Gummel (SG) scheme with Fermi-Dirac (FD) statistics, and (c) when correctly using the SEDAN scheme with Fermi-Dirac statistics.</td>
      <td>194</td>
    </tr>
    <tr>
      <td>C.2</td>
      <td>Numerical electron flux averaged over each atomic plane at a bias of 3.3V shown for Boltzmann statistics (black), Fermi-Dirac statistics incorrectly using the Scharfetter-Gummel (SG) flux discretization (red) and Fermi-Dirac statistics correctly using the SEDAN flux discretization (blue).</td>
      <td>195</td>
    </tr>
    <tr>
      <td>C.3</td>
      <td>Electron density at a bias of 3.3V (a) when using Boltzmann statistics, (b) when incorrectly using the Scharfetter-Gummel (SG) scheme with Fermi-Dirac statistics, and (c) when correctly using the SEDAN scheme with Fermi-Dirac statistics.</td>
      <td>196</td>
    </tr>
    <tr>
      <td>C.4</td>
      <td>(a) Difference in magnitude of the Shockley-Read-Hall (SRH, black), radiative (red) and Auger (blue) recombination rates between Fermi-Dirac and Boltzmann statistics at a bias of 3.3V, calculated as described in the main text. (b) Current density-voltage (IV) curves using Boltzmann statistics (black, dashed), Fermi-Dirac statistics using the Scharfetter-Gummel scheme (SG, red) and Fermi-Dirac statistics using the SEDAN scheme (black, solid), shown on a log scale.</td>
      <td>197</td>
    </tr>
    <tr>
      <td>D.1</td>
      <td>Schematic illustration of the (Al,Ga)N-based p-i-n system underlying the simulations. The intrinsic barriers plus quantum well region is denoted as the "atomistic" mesh in the main text, while the n- and p-regions are described by a sparser device mesh.</td>
      <td>201</td>
    </tr>
    <tr>
      <td>D.2</td>
      <td>(a) Current-voltage curves, on a log-scale, for an Al<sub>0.7</sub>Ga<sub>0.3</sub>N/Al<sub>0.8</sub>Ga<sub>0.2</sub>N quantum well embedded in a p-n junction. The data are shown for simulations that (i) account for random alloy (RA) fluctuations in the well and barriers (blue), (ii) account for RA fluctuations in the well but treat the barrier in a virtual crystal approximation (VCA) (red), and (iii) use a VCA in the entire p-i-n structure (pink). The inset shows the I-V curves on a linear scale.</td>
      <td>202</td>
    </tr>
    <tr>
      <td>D.3</td>
      <td>3D isosurface plot of the current density in the intrinsic (Al,Ga)N quantum well and barrier regions at a bias of 5.5 V, using the fully atomistic description. The current density for the holes (blue) and electrons (red) are plotted at 700 A cm<sup>-2</sup>, reflecting locally high current densities due to alloy fluctuations; dashed black lines indicate well boundaries.</td>
      <td>203</td>
    </tr>
  </tbody>
</table>
# Avant-propos

Ce manuscrit présente les résultats obtenus durant la thèse réalisée au sein du centre Inria de l'Université de Lille et encadrée par Claire Chainais-Hillairet, Maxime Herda et Simon Lemaire. Ce document contient une introduction générale, qui décrit les modèles et méthodes numériques étudiés et dresse un panorama des contributions (en anglais). Cette introduction est suivie de quatre chapitres, correspondant chacun à des travaux publiés, et d'une conclusion dans laquelle nous exposons quelques perspectives. Le manuscrit est complété par quatre chapitres annexes, qui résultent de travaux réalisés en collaborations avec des mathématiciens et physiciens du Weierstrass Institute (Berlin, Allemagne) et du Tyndall National Institue (Cork, Irelande). L'Annexe A introduit les travaux et leur contextes physique et mathématique, tandis que les Annexes B-C-D rassemblent des travaux publiés ou soumis pour publication.

## Liste des publications :

Chapitre 1 : *Long-time behaviour of hybrid finite volume schemes for advection-diffusion equations : linear and nonlinear approaches* ; en collaboration avec C. Chainais-Hillairet, M. Herda et Simon Lemaire ; paru [70] dans *Numerische Mathematik*, 151, 963-1016.

Chapitre 2 : *A structure preserving hybrid finite volume scheme for semiconductor models with magnetic field on general meshes* ; paru [213] dans ESAIM : Mathematical Modelling and Numerical Analysis, 57, 2557-2593.

Chapitre 3 : *Structure-preserving schemes for drift-diffusion systems on general meshes : DDFV vs HFV* ; en collaboration avec S. Krell ; accepté pour publication [187] comme proceeding de la conférence *Finite Volume for Complex Applications 10*.

Chapitre 4 : version étendue du proceeding *A skeletal high-order structure-preserving scheme for advection-diffusion equations* ; accepté pour publication [212] à la conférence *Finite Volume for Complex Applications 10*.

Annexe B : *Impact of random alloy fluctuations on the carrier distribution in multi-color (In,Ga)N/GaN quantum well systems* ; en collaboration avec M. O'Donovan, P. Farrell, T. Streckenbach, T. Koprucki et S. Schulz ; soumis pour publication [219].

Annexe C : *Importance of satisfying thermodynamic consistency in optoelectronic device simulations for high carrier densities* ; en collaboration avec P. Farrell, M. O'Donovan, S. Schulz et T. Koprucki ; paru [126] dans *Optical and Quantum Electronics*, 55, 978.

Annexe D : *Theoretical study of the impact of alloy disorder on carrier transport and recombination processes in deep UV (Al,Ga)N light emitters* ; en collaboration avec R. Finn, M. O'Donovan, P. Farrell, T. Streckenbach, T. Koprucki et S. Schulz ; paru [133] dans *Applied Physics Letters* 122 (24): 241104.
# Introduction générale

## Sommaire du présent chapitre

<table><tr><td><strong>Phénomènes de convection-diffusion : modélisation, comportement et analyse</strong></td><td style="text-align: right;">4</td></tr><tr><td>Modèles anisotropes d'advection-diffusion et de dérive-diffusion</td><td style="text-align: right;">4</td></tr><tr><td>Positivité des solutions et conservation de la masse</td><td style="text-align: right;">7</td></tr><tr><td>Comportement en temps long et méthode d'entropie</td><td style="text-align: right;">8</td></tr><tr><td><strong>Approximation numérique des problèmes de convection-diffusion</strong></td><td style="text-align: right;">11</td></tr><tr><td>Enjeux de fiabilité de l'approximation</td><td style="text-align: right;">11</td></tr><tr><td>Volumes finis à deux points : forces et faiblesses</td><td style="text-align: right;">13</td></tr><tr><td>Schémas pour problèmes anisotropes sur maillages généraux</td><td style="text-align: right;">14</td></tr><tr><td>Méthodes d'ordre élevé</td><td style="text-align: right;">15</td></tr><tr><td>Méthodes étudiées dans cette thèse : HFV et HHO</td><td style="text-align: right;">16</td></tr><tr><td><strong>Overview of the works and results</strong></td><td style="text-align: right;">18</td></tr><tr><td>Long-time behaviour of Hybrid Finite Volume schemes for advection-diffusion (Chapter 1)</td><td style="text-align: right;">19</td></tr><tr><td>Structure-preserving schemes for semiconductor models (Chapters 2 and 3)</td><td style="text-align: right;">20</td></tr><tr><td>High-order schemes for advection-diffusion (Chapter 4)</td><td style="text-align: right;">23</td></tr><tr><td>Semiconductors models with irregular convection fields (Appendices A-B-C-D)</td><td style="text-align: right;">27</td></tr></table>

Cette thèse est dédiée à l'approximation numérique de problèmes de convection-diffusion, éventuellement anisotropes, par des schémas d'ordre élevé sur maillage généraux. L'originalité principale de ces travaux est de traiter simultanément plusieurs enjeux : préservation de la structure, maillages généraux et ordre élevé. Bien que de nombreuses avancées aient été réalisées pour aborder chacun de ces enjeux individuellement, concevoir des méthodes qui présentent toutes ces caractéristiques simultanément est un travail difficile. Notre objectif est donc de proposer des méthodes fiables et efficaces : les solutions numériques doivent préserver les propriétés physiques des solutions continues tout en autorisant une large gamme de paramètres de discrétisation et en maintenant des coûts de calcul raisonnables.

Dans cette introduction générale, nous présentons dans un premier temps les modèles considérés dans ce manuscrit, en partant de différents contextes physiques, et en décrivant les propriétés (positivité et comportement en temps long) de leurs solutions. Cette présentation est l'occasion d'introduire la méthode d'entropie, outil fondamental des travaux de cette thèse. Nous nous intéressons ensuite à l'approximation numérique de ces modèles, en introduisant les grands concepts et enjeux de cette procédure, puis en proposant un état de l'art des différentes méthodes d'approximation pour les modèles de convection-diffusion. Enfin, nous dressons un
panorama (en anglais) des travaux théoriques (analyse numérique de schémas) et numériques (implémentation de schémas, investigation numérique de leurs propriétés et comparaison) réalisés au cours de cette thèse.

## Phénomènes de convection-diffusion : modélisation, comportement et analyse

Les problèmes de convection-diffusion étudiés dans cette thèse constituent une brique de base pour décrire des phénomènes physiques et biologiques variés. Parmi les différents modèles complexes s’inscrivant dans le cadre discuté, on peut citer les modélisations d’écoulements en milieux poreux [17], l’étude de dispositif électronique comme les semi-conducteurs [258, 207, 208], les modélisations en biologie [135, 183, 179, 246], l’étude de la corrosion et de la dégradation de certains matériaux [3, 16, 53] ou encore des modèles de mécanique des fluides [206]. Tous ces phénomènes partagent une structure dissipative semblable. Ils échangent quelque chose avec leur environnement au cours du temps. Sous l’effet de ces échanges, le système va évoluer en vue d’atteindre un état plus stable. D’un point de vue mathématique, les modèles cités présentent tous un caractère diffusif [138, 130] et appartiennent à la grande classe des Équations aux Dérivées Partielles (EDP) paraboliques.

## Modèles anisotropes d’advection-diffusion et de dérive-diffusion

Le premier modèle qui nous intéresse est une équation linéaire d’advection-diffusion, étudiée dans les Chapitres 1 et 4. Il s’agit d’un modèle dissipatif jouet, qui contient la structure de base des modèles dissipatifs. Il décrit l’évolution au cours du temps d’une quantité de matière dans un domaine borné ($\Omega \subset \mathbb{R}^d, d \in \{1, 2, 3\}$) qui est soumise à deux phénomènes : la diffusion d’une part (transfert naturel des zones de forte occupation vers celles faiblement occupées) et l’advection par un champ extérieur $V : \Omega \to \mathbb{R}^d$ (transport de la matière à cause d’une force externe). Le système mathématique obtenu pour l’inconnue notée $u : \mathbb{R}_+ \times \Omega \to \mathbb{R}$, qui représente la quantité de matière (ou sa concentration/densité), est le suivant :

$$ egin{cases} 
\partial_t u - \text{div}(J(u)) = f & \text{dans } \mathbb{R}_+ \times \Omega, \\ 
J(u) = \Lambda(\nabla u + uV) & \text{dans } \mathbb{R}_+ \times \Omega, \\ 
u = u^D & \text{sur } \mathbb{R}_+ \times \Gamma^D, \\ 
J(u) \cdot n = g^N & \text{sur } \mathbb{R}_+ \times \Gamma^N, \\ 
u(0, \cdot) = u^{in} & \text{dans } \Omega. 
\end{cases} 
\quad (1) $$

Dans ce système, $u^{in}$ est la donnée initiale, $f$ est le terme source qui traduit l'apparition ou la disparition de matière dans le domaine. Les conditions de bord, de type Neumann sur $\Gamma^N$ et Dirichlet sur $\Gamma^D$ correspondent respectivement à une condition sur les flux au bord et une imposition des quantités de matière présentes. Le modèle est caractérisé par le flux d'advection-diffusion $J(u) = \Lambda(\nabla u + uV)$ qui décrit les échanges locaux de matière suivant la loi de Nernst-Planck [218, 227]. L'une des particularité de cette thèse est de s'intéresser à des modèles anisotropes : le flux dépend d'un tenseur $\Lambda : \Omega \to \mathbb{R}^{d \times d}$, qui traduit l'existence de directions privilégiées pour la diffusion, généralement à cause de la structure du milieu (comme c'est le cas dans les applications aux milieux poreux [17]) ou d'une force extérieure (comme dans les modèles de semi-conducteurs [143] présentés ci-dessous). Dans de nombreuses applications,
le champ d'advection dérive d'un potentiel $\phi : \Omega \to \mathbb{R}$ : on peut écrire $V = \nabla \phi$, et

$$J(u) = \Lambda(\nabla u + u \nabla \phi).$$

Dans cette thèse, on s'intéresse également à des modèles de transport de charge dans des matériaux semi-conducteurs. Ces matériaux sont les constituants de base de nombreux dispositifs électroniques, notamment sous forme de diodes et de transistors qui constituent les circuits intégrés (et donc, par exemple, les processeurs et mémoires des ordinateurs ou des téléphones). Ils sont aussi utilisés dans des dispositifs opto-électroniques, comme les Diodes Électroluminescentes (DEL, light-emitting diode LED en anglais) ou les cellules photovoltaïques. Le modèle le plus simple pour les étudier est le système de dérive-diffusion introduit par Van Roosbroeck dans [258] en 1950. Il modélise l'évolution de charges électriques dans un matériau possédant une structure de réseau cristallin, qui constitue le semi-conducteur. L'idée fondamentale régissant ces structures est la suivante : les électrons (charges négatives) peuvent facilement se déplacer au sein du matériau, qui possède ainsi des propriétés intermédiaires entre un isolant électrique et un conducteur. Le transport des électrons laisse libres certains sites dans les bandes de valence du semi-conducteur : ce sont les trous d'électrons, aussi appelés simplement trous, qui correspondent à des charges positives. Au niveau mathématique, après adimensionnement, on obtient un système couplé d'inconnues N, P et $\phi$ vérifiant

$$ egin{cases} 
\partial_t N - \text{div}(\nabla N - N \nabla \phi) = 0 & \text{dans } \mathbb{R}_+ \times \Omega, \\ 
\partial_t P - \text{div}(\nabla P + P \nabla \phi) = 0 & \text{dans } \mathbb{R}_+ \times \Omega, \\ 
-\lambda^2 \text{div}(\nabla \phi) = C + P - N & \text{dans } \mathbb{R}_+ \times \Omega, \\ 
N = N^D, P = P^D \text{ et } \phi = \phi^D & \text{sur } \mathbb{R}_+ \times \Gamma^D, \\ 
(\nabla N - N \nabla \phi) \cdot n = (\nabla P + P \nabla \phi) \cdot n = \nabla \phi \cdot n = 0 & \text{sur } \mathbb{R}_+ \times \Gamma^N, \\ 
N(0, \cdot) = N^{in} \text{ et } P(0, \cdot) = P^{in} & \text{dans } \Omega. 
\end{cases} 
\quad (2) $$

Dans ce système, $N : \mathbb{R}_+ \times \Omega \to \mathbb{R}$ représente la densité de charges négatives (électrons), $P : \mathbb{R}_+ \times \Omega \to \mathbb{R}$ la densité de charges positives (trous) et $\phi : \mathbb{R}_+ \times \Omega \to \mathbb{R}$ est un potentiel électrostatique, induit par l'inhomogénéité spatiale de charge électrique dans le domaine. Les densités de charges suivent des lois de conservation décrites par des équations d'advection-convection, dont le champ $V_E = \pm \nabla \phi$ dérive du potentiel électrostatique et dépend de la charge de l'espèce considérée. Ce potentiel est lui-même solution d'une équation de Poisson, qui prend en compte la densité de charges totale, constituée de plusieurs contributions : les trous, les électrons et les impuretés (charges intrinsèques ajoutées volontairement dans le réseau cristallin). Cet ajout de charge dans le matériau est nommé dopage et correspond à la fonction $C : \Omega \to \mathbb{R}$. Le dopage est un processus d'une importance capitale, car il permet de créer un dispositif aux propriétés voulues. Ainsi, le dopage C caractérise en partie le type de dispositif utilisé (diode, transistor, etc). Une illustration de géométrie typique de diode et son dopage associé est donnée en Figure 1, page 12.

Dans l'équation de Poisson qui régit le potentiel $\phi$, la constante $\lambda$ est nommée longueur de Debye (réduite). Elle est issue de l'adimensionnement du système, et dépend de certaines constantes physiques fondamentales (la conductivité diélectrique du milieu, la charge et la mobilité d'une particule élémentaire) et des caractéristiques du semi-conducteur étudié, dont la norme $L^\infty$ du dopage et la taille du dispositif. Les valeurs typiques de $\lambda$ varient de $10^{-3}$ à $10^{-5}$.

Pour les modèles de semi-conducteurs, les conditions de bord sont des conditions mixtes Dirichlet-Neumann : les conditions de Dirichlet modélisent des contacts ohmiques (i.e. un contact avec un métal, conducteur, qui va imposer la quantité de charges au bord, ou avec
d'autres semi-conducteurs) et celles de Neumann correspondent à des parties isolées (sans contact avec l'extérieur, les charges ne peuvent ni rentrer ni sortir).

Le modèle de Van Roosbroeck est un modèle simple, auquel ont peut ajouter des effets plus complexes. On peut notamment s'intéresser à l'effet d'un champ magnétique extérieur [143], à la prise en compte de diffusions non-linéaires [144, 257] ou encore prendre en compte les interactions entre les charges via des termes de réaction non-linéaires. La motivation principale de cette thèse est l'approximation du système ci-dessous, qui prend en compte ces différents effets :

$$
\left\{
\begin{array}{ll}
\partial_t N - \operatorname{div}(N \Lambda_N \nabla(h(N) - \phi)) = -R(N, P) & \text{dans } \mathbb{R}_+ \times \Omega \\
\partial_t P - \operatorname{div}(P \Lambda_P \nabla(h(P) + \phi)) = -R(N, P) & \text{dans } \mathbb{R}_+ \times \Omega \\
\qquad - \operatorname{div}(\Lambda_\phi \nabla \phi) = C + P - N & \text{dans } \mathbb{R}_+ \times \Omega \\
N = N^D, P = P^D \text{ et } \phi = \phi^D & \text{sur } \mathbb{R}_+ \times \Gamma^D \\
N \Lambda_N \nabla(h(N) - \phi) \cdot n = P \Lambda_P \nabla(h(P) + \phi) \cdot n = \Lambda_\phi \nabla \phi \cdot n = 0 & \text{sur } \mathbb{R}_+ \times \Gamma^N \\
N(0, \cdot) = N^{in} \text{ et } P(0, \cdot) = P^{in} & \text{dans } \Omega,
\end{array}
\right.
\quad (3)
$$

La nonlinéarité *h* dans les équations sur les charges est obtenue à partir de modèles ciné-
tiques [172]. On peut alors décrire les flux sur les charges comme étant le produit de densités
par des potentiels dits de quasi-Fermi *w*<sub>N</sub> et *w*<sub>P</sub> :

$$
J_N = -N \nabla w_N \quad \text{et} \quad J_P = -P \nabla w_P.
$$

Ces potentiels jouent un rôle fondamental dans le modèle, tant d'un point de vue physique que
mathématique. Ils sont liés aux densités des espèces à travers les équations d'état

$$
N = g(\phi - w_N) \quad \text{et} \quad P = g(-\phi + w_P), \tag{4}
$$

où la fonction $g : \mathbb{R} \rightarrow \mathbb{R}_+$ est la fonction de distribution associée à la statistique. C'est une fonction strictement monotone, de bijection réciproque $h$. On peut donc exprimer les potentiels de quasi-Fermi en fonction des densités

$$
w_N = h(N) - \phi \quad \text{et} \quad w_P = h(P) + \phi,
$$

et obtenir les équations de (3). Dans le cas du modèle de Van Roosbroeck, la statistique considérée est celle de Boltzmann, qui correspond à $g = \exp$ et $h = \log$; les termes de diffusion qui en découlent sont linéaires. Le type de statistique utilisée dans la modélisation dépend de la nature du semi-conducteur (organique ou non) et du cadre d'application [144, 257, 125]. Les expressions de $g$ sont généralement complexes, mais il existe des approximations analytiques couramment utilisées [31, 223].

La présence des tenseurs d'anisotropie Λ_N et Λ_P est liée à l'existence d'un champ magnétique extérieur, par exemple dans le cadre d'étude de capteurs [234]. Ce champ va faire tourner les charges et ainsi modifier le comportement global du semi-conducteur. Un modèle pour de tels dispositifs a été introduit (dans le cas de la statistique de Boltzmann) par Gajewski et Gärtner dans [143], puis obtenu comme limite rigoureuse d'un modèle cinétique [23, 155]. Le champ magnétique extérieur est décrit par un champ de vecteur B : Ω → ℝ<sup>d</sup>, qui modifie les flux d'électrons et de trous en rajoutant un terme dû à la contribution magnétique de la force de Lorentz. Des manipulations algébriques permettent ensuite d'exprimer directement les flux avec champ magnétique comme des produits de tenseurs Λ_N et Λ_P par les flux de charges sans champ magnétique. Pour un problème bidimensionnel, l'effet de rotation ne peut avoir lieu que dans le plan du semi-conducteur, et on peut se restreindre au cas B = be<sub>z</sub>, où e<sub>z</sub> est un vecteur
orthogonal au semi-conducteur et $b$ est l'intensité du champ magnétique. Les tenseurs sont alors
de la forme

$$
\Lambda_N = \frac{1}{1+b^2} \begin{pmatrix} 1 & b \\ -b & 1 \end{pmatrix} \quad \text{et} \quad \Lambda_P = \frac{1}{1+b^2} \begin{pmatrix} 1 & -b \\ b & 1 \end{pmatrix}. \tag{5}
$$

Pour un problème en trois dimensions de l'espace, en notant $B = (b_1, b_2, b_3)$, on a

$$
\Lambda_N = \frac{1}{1 + |B|^2} \begin{pmatrix} 1 + b_1^2 & b_1 b_2 - b_3 & b_1 b_3 + b_2 \\ b_1 b_2 + b_3 & 1 + b_2^2 & b_2 b_3 - b_1 \\ b_1 b_3 - b_2 & b_2 b_3 + b_1 & 1 + b_3^2 \end{pmatrix}
$$

et

$$
\Lambda_P = \frac{1}{1 + |B|^2} \begin{pmatrix} 1 + b_1^2 & b_1 b_2 + b_3 & b_1 b_3 - b_2 \\ b_1 b_2 - b_3 & 1 + b_2^2 & b_2 b_3 + b_1 \\ b_1 b_3 + b_2 & b_2 b_3 - b_1 & 1 + b_3^2 \end{pmatrix}.
$$

L'une des particularités de ces tenseurs est leur caractère non-symétrique.

Dans les équations de transport de charges, les termes $R(N,P)$ sont des termes de réaction, qui traduisent des interactions entre les charges et avec le milieu. En particulier, à cause de l'agitation thermique interne ou d'une stimulation externe (par exemple optique ou électrique) il existe un phénomène d'apparition de charges dans le dispositif, nommé génération : les électrons de valence sont arrachés pour devenir des électrons de conduction (libres). D'autre part, les électrons et les trous peuvent interagir entre eux via un phénomène nommé recombinaison : un électron de conduction peut devenir un électron de valence et neutraliser un trou. Ces termes ont une structure particulière, de la forme

$$
R(N,P) = r(N,P) (e^{h(N)+h(P)} - 1), \quad (6)
$$

où $r$ est une fonction positive. Ces effets de génération/recombinaison sont fondamentaux
lorsque l'on s'intéresse à des dispositifs opto-électroniques, et la compétition entre eux a une
influence sur le rendement des dispositifs. Ces questions font l'objet de l'Annexe D.

Positivité des solutions et conservation de la masse

L'équation linéaire (1) a été longuement étudiée et son analyse est maintenant classique [119, Section 7.1]. On sait en particulier montrer l'existence et l'unicité d'une solution globale, et évaluer sa régularité en fonctions des données.

Parmi les propriétés remarquables des solutions figure la positivité :

$$
\text{si les données } u^{in}, f, g^N \text{ et } u^D \text{ sont positives, alors } u \ge 0 \text{ dans } \mathbb{R}_+ \times \Omega.
$$

Il suffit d'imposer en plus que $\int_{\Omega} u^{in} > 0$ pour assurer la stricte positivité de la solution. Ces
propriétés sont essentielles pour assurer une cohérence avec la physique du problème.

Une autre propriété remarquable concerne la préservation de la masse : si l'on impose des
conditions de Neumann sur tout le bord du domaine ($\Gamma^N = \partial\Omega$), alors sous une hypothèse de
compatibilité entre les données $\int_{\Omega} f + \int_{\partial\Omega} g^N = 0$, la masse de la solution est conservée au cours
du temps :

$$
\forall t > 0, \quad \int_{\Omega} u(t) = \int_{\Omega} u^{in}.
$$

Concernant les systèmes de dérive-diffusion, l'existence et l'unicité de solutions globales ont été prouvées dans divers contextes, sous des hypothèses raisonnables sur les données (dont la
positivité des densités initiales et des densités sur les contacts ohmiques). On peut notamment
citer les travaux fondateurs [214, 142, 208] pour le modèle de Van Roosbroeck (2), des études
pour des statistiques non-Boltzmann [144, 149] et une analyse pour un modèle anisotrope [143].
L'étude du système (3) peut être abordée en suivant les arguments présentés dans ces travaux.
Comme pour le modèle d'advection-diffusion, les solutions des systèmes de dérive-diffusion
satisfont une propriété de positivité : les densités de charges N et P sont strictement positives
sur le domaine en tout temps. Pour certaines statistiques non-Boltzmann (typiquement une
statistique de Blakemore qui correspond à une non-linéarité de la forme h(s) = log(s/(1-s))), on a
également une borne supérieure a priori sur les densités donnée par a = sup<sub>R</sub> h. Ainsi, en général

$$
\forall t > 0, 0 < N(t) < a \text{ et } 0 < P(t) < a,
$$

$$
\text{où } a \in ]0, +\infty[.
$$

Comportement en temps long et méthode d'entropie

Comme annoncé précédemment, les phénomènes considérés sont caractérisés par une évolu-
tion vers des états stationnaires. L'objectif de cette section est de présenter la méthode d'analyse
emblématique pour étudier ce phénomène.

Pour l'équation d'advection-diffusion (1), l'état stationnaire associé $u^{\infty}$ est défini par

$$
\left\{
\begin{array}{ll}
-\operatorname{div}(\Lambda(\nabla u^{\infty} + u^{\infty}V)) = f & \text{dans } \Omega, \\
u^{\infty} = u^{D} & \text{sur } \Gamma^{D}, \\
(\Lambda(\nabla u^{\infty} + u^{\infty}V)) \cdot n = g^{N} & \text{sur } \Gamma^{N}.
\end{array}
\right.
\qquad (7)
$$

Dans le cas où Γ<sup>D</sup> = ∅, pour assurer une cohérence avec le problème évolutif, la masse de u<sup>∞</sup> doit aussi être égale à la masse initiale :

$$
\int_{\Omega} u^{\infty} = \int_{\Omega} u^{in} = M.
$$

L'état stationnaire est bien défini et unique [106, 114]. De plus, il a les mêmes propriétés de
positivité que la solution transitoire.

Notons qu'il existe des états stationnaire particuliers, qui correspondent à des états physiques
très stables, pour lesquels les flux sont nuls. Ces états stationnaires sont nommés équilibres
thermiques. Lorsque le champ dérive d'un potentiel $V = \nabla\phi$, les équilibres sont associés à des
potentiels de quasi-Fermi constants sur $\Omega$. En pratique, pour l'équation (1) de tels état station-
naires ne peuvent exister que si les données sont compatibles avec l'équilibre. Typiquement, un
terme source nul $f = 0$ et des conditions de Neumann homogènes sur tout le bord $(\Gamma^N = \partial\Omega$
avec $g^N = 0$) donnent un équilibre thermique. Lorsque ces conditions ne sont pas satisfaites, on
obtient un état stationnaire dit hors-équilibre.

Pour le modèle d'advection-diffusion (1), la solution $u(t)$ converge exponentiellement vite, quand $t \to +\infty$, vers $u^\infty$, solution de (7). Un tel résultat se généralise à d'autre modèles dissipatifs, dont les systèmes de dérive-diffusion (2)-(3). Pour démontrer un tel comportement, on utilise la méthode d'entropie. Celle-ci a initialement été introduite dans le cadre de la théorie cinétique des gaz par Boltzmann dans [34], à travers le théorème H. Comme présenté dans [8], cette méthode a par la suite été utilisée dans des cadres plus variés, notamment pour l'étude d'équations d'advection-diffusion [62, 254, 63, 33], ou plus généralement d'équations paraboliques [61]. On renvoie à l'ouvrage de Jüngel [174] pour une explication détaillée, et décrivons ci-dessous la
philosophie générale de la méthode.

(i) **Identification d'une entropie relative** : on définit une quantité, motivée par la physique, qui dépend de la solution et de l'état stationnaire, nommée entropie relative $\mathbb{E}$. Cette quantité est positive, et doit mesurer dans un certain sens l'écart entre $u$ et $u^\infty$. En particulier, l'entropie relative s'annule quand la solution transitoire coïncide avec l'état stationnaire.

(ii) **Obtention d'une relation d'entropie** : on dérive l'entropie relative par rapport au temps afin d'obtenir une relation sur l'entropie, parfois nommée relation d'entropie/dissipation d'entropie, ou simplement relation de dissipation d'entropie, de la forme :

$$ \frac{d}{dt} \mathbb{E}(t) = -\mathbb{D}(t), \qquad (8) $$

où $\mathbb{D}$ est une quantité positive nommée dissipation d'entropie ou information de Fischer. Cette relation implique la décroissance de l'entropie au cours du temps : l'entropie relative s'interprète alors comme une fonctionnelle de Lyapunov du problème.

(iii) **Obtention d'un contrôle de l'entropie par sa dissipation** : on compare l'entropie relative avec sa dissipation, pour obtenir une estimation du type

$$ \forall t > 0, \nu \mathbb{E}(t) \le \mathbb{D}(t), \qquad (9) $$

où $\nu$ est une constante positive. Cette estimation s'obtient à l'aide d'inégalités fonctionnelles. On en déduit une décroissance exponentielle de l'entropie :

$$ \forall t > 0, \mathbb{E}(t) \le e^{-\nu t} \mathbb{E}(0). \qquad (10) $$

(iv) **Obtention d'un contrôle de la distance entre $u$ et $u^\infty$ par l'entropie** : pour conclure, on trouve une norme $\|\cdot\|$ sur un espace adapté tel que

$$ \forall t > 0, \|u(t) - u^\infty\|^\gamma \le C \mathbb{E}(t), \qquad (11) $$

où $\gamma$ et $C$ sont des constantes strictement positives. On en déduit alors de (10) une convergence exponentiellement rapide de $\|u(t) - u^\infty\|$ vers 0 quand $t \to \infty$.

Pour l'équation (1), on peut par exemple partir de l'énergie quadratique pour définir une entropie relative

$$ \mathbb{E}_2 = \frac{1}{2} \|u - u^\infty\|_{L^2(\Omega)}^2, \text{ de dissipation } \mathbb{D}_2 = \int_{\Omega} \Lambda \nabla (u - u^\infty) \cdot \nabla (u - u^\infty). $$

Une autre possibilité est d'utiliser l'entropie de Boltzmann (ou plus généralement des entropies de Tsallis [33])

$$ \mathbb{E}_B = \int_{\Omega} u^{\infty} \Phi_1 \left( \frac{u}{u^{\infty}} \right), \text{ de dissipation } \mathbb{D}_B = \int_{\Omega} u(t) \Lambda \nabla (w - w^{\infty}) \cdot \nabla (w - w^{\infty}) $$

où $\Phi_1(s) = s \log(s) - s + 1$, $w = \log(u) + \phi$ et $w^\infty = \log(u^\infty) + \phi$. Cette entropie correspond à une quantité physique nommée énergie interne. Notons que la relation d'entropie (8) donne également des bornes sur l'entropie et sa dissipation. Il est souvent possible d'en déduire des estimations sur la solution elle-même. Cette stratégie est l'un des outils essentiels de cette thèse, permettant d'obtenir des estimations a priori sur des solutions discrètes.

Pour les modèles de semi-conducteur, on peut également étudier le comportement en temps long via une méthode d'entropie. La notion d'état stationnaire étant plus délicate à définir
(notamment de part leur non-unicité [215, 216]), l'analyse du comportement en temps long est restreinte à un cadre simple où l'état stationnaire est un équilibre thermique ($N^e$, $P^e$, $\phi^e$). Comme pour l'équation d'advection-diffusion, celui-ci est caractérisé par des potentiels de quasi-Fermi constants :

$$h(N^e) - \phi^e = \alpha_N \text{ et } h(P^e) + \phi^e = \alpha_P \text{ dans } \Omega,$$

où $(\alpha_N, \alpha_P) \in \mathbb{R}^2$. Cette relation impose naturellement une condition sur les données au bord :

$$h(N^D) - \phi^D = \alpha_N \text{ et } h(P^D) + \phi^D = \alpha_P \text{ sur } \Gamma^D. \qquad (12)$$

En lien avec l'intuition physique, on impose également qu'il n'y ait plus de phénomène de recombinaison/génération lorsque le semi-conducteur est à l'équilibre, ce qui signifie que $R(N^e, P^e) = 0$. En revenant à la définition (6), ceci impose une autre condition de compatibilité dans le cas où $r \neq 0$ :

$$\alpha_N + \alpha_P = 0. \qquad (13)$$

Sous ces conditions, on obtient une relation simple entre les densités de charge et le potentiel à l'équilibre :

$$N^e = g(\alpha_N + \phi^e) \quad \text{et} \quad P^e = g(\alpha_P - \phi^e).$$

On en déduit que le potentiel à l'équilibre $\phi^e$ est une solution de l'équation de Poisson non-linéaire (nommée équation de Poisson-Boltzmann dans le cas de la statistique de Boltzmann $g = \exp$) suivante :

$$\left\{ \begin{array}{ll} -\mathrm{div}(\Lambda_{\phi} \nabla \phi^e) = C + g(\alpha_P - \phi^e) - g(\alpha_N + \phi^e) & \text{dans } \Omega, \\ \phi^e = \phi^D & \text{sur } \Gamma^D, \\ \Lambda_{\phi} \nabla \phi^e \cdot n = 0 & \text{sur } \Gamma^N. \end{array} \right. \qquad (14)$$

Les propriétés de $g$ permettent alors de montrer que (14) admet une unique solution $\phi^e$. Ceci implique qu'il existe un unique équilibre thermique ($N^e$, $P^e$, $\phi^e$). Sous ces hypothèses, l'analyse du comportement en temps long est alors menée pour différentes situations grâce à une méthode d'entropie [144, 143]. L'entropie relative utilisée pour le système (3) est définie par

$$\mathbb{E}(t) = \int_{\Omega} H(N) - H(N^e) - h(N^e)(N - N^e) + \int_{\Omega} H(P) - H(P^e) - h(P^e)(P - P^e) \\ + \frac{1}{2} \int_{\Omega} \nabla(\phi - \phi^e) \cdot \Lambda_{\phi} \nabla(\phi - \phi^e),$$

où $H: x \to \int_1^x h(s)ds$. Les deux premiers termes sont des énergies internes (entropies de type Boltzmann sur les densités), et le dernier terme est une énergie électrique. On a alors la relation de dissipation d'entropie

$$\frac{d\mathbb{E}}{dt}(t) = -\mathbb{D}(t),$$

où le terme de dissipation d'entropie s'écrit

$$\mathbb{D}(t) = \int_{\Omega} N \nabla(h(N) - \phi) \cdot \Lambda_N \nabla(h(N) - \phi) + \int_{\Omega} P \nabla(h(P) + \phi) \cdot \Lambda_P \nabla(h(P) + \phi) \\ + \int_{\Omega} R(N, P) (h(N) + h(P)).$$
Le dernier terme est lié aux processus de recombinaison/génération. Il est positif en vertu de la structure (6) des termes de réaction.

## Approximation numérique des problèmes de convection-diffusion

Afin de comprendre les enjeux de l'approximation numérique, il est bon de rappeler que pour les modèles jouets (1), (2) et (3), les résultats théoriques ne fournissent pas d'expression explicite de la solution. C'est d'autant plus vrai pour des modèles plus réalistes utilisés dans l'industrie ou la recherche (tels que les modèles utilisés actuellement par des équipes de recherche en opto-électronique, c.f. Annexes B-C-D). Pour pouvoir exploiter pleinement les modèles mathématiques proposés, il faut donc être capable d'en tirer des informations quantitatives, en calculant une solution approchée qui pourra être exploitée au niveau industriel. A titre d'exemple, la conception de dispositifs modernes comme les LEDs émettant en UV profond (étudiées en Annexe D) nécessite beaucoup de tests par essai-erreur avant d'obtenir des résultats probants (par exemple, pour obtenir des LEDs efficaces en termes de consommation d'énergie). De tels tests ne peuvent se faire expérimentalement, à cause du nombre de configurations et du coût (économique et temporel) de chaque dispositif. La simulation numérique est alors la solution privilégiée pour réaliser ces tests à moindre coût, mais il faut être capable de garantir que la solution calculée est correcte.

Le travail réalisé dans cette thèse s'inscrit dans cette logique : nous proposons des méthodes numériques pour approcher les solutions des modèles (1), (2) et (3). Les méthodes introduites sont alors analysées théoriquement et étudiées au niveau numérique.

Dans cette section, nous présentons les enjeux majeurs de l'analyse numérique, et dressons un panorama des méthodes existantes pour approcher les problèmes étudiés dans cette thèse. L'accent est mis sur la question de la préservation des propriétés des solutions continues au niveau discret.

## Enjeux de fiabilité de l'approximation

Pour les méthodes numériques qui nous intéressent, les solutions approchées sont obtenues via la discrétisation du domaine spatial en un maillage (partition en sous-domaines polyédriques nommés mailles ou volumes de contrôle). Dans l'idéal, on souhaite avoir une grande liberté lors de la conception du maillage, pour pouvoir par exemple avoir des mailles de formes différentes, ou des nœuds orphelins. Les questions relatives aux maillages sont crucales dans certaines applications, comme en géosciences où le maillage est souvent imposé par la connaissance que l'on a du sol : il s'agit d'une donnée du problème. Un autre intérêt des maillages généraux est la facilité avec laquelle ils permettent de réaliser des raffinements locaux. Même pour des simulations de semi-conducteurs (associés à des dispositifs aux géométries assez simples, souvent cartésiennes), il peut être très intéressant d'avoir un maillage adapté localement (voir Figure 1 avec le cas d'une diode où les zones d'intérêt sont raffinées) afin d'optimiser les coûts de calculs. Pour finir, les capacités de prise en compte d'un maillage général ou de la diffusion anisotrope sont fortement liées.

La qualité d'un schéma numérique est souvent évaluée par la précision du schéma (sur un maillage fixé). Cependant, un tel critère ne saurait être suffisant pour garantir une approximation fiable des phénomènes simulés. En effet, l'analyse théorique des modèles continus fournit des propriétés qualitatives importantes, qu'il convient de préserver au niveau discret. Les schémas qui suivent ce principe sont nommés de "schémas préservant la structure" (structure-preserving schemes en anglais). L'enjeu majeur de cette thèse est le développement de schémas préservant la structure d'entropie des problèmes de convection-diffusion, et ce sans condition sur le maillage
<figure><img src="image_9.png" /><figcaption>FIGURE 1 – Géométrie d’une diode et maillage adapté associé.</figcaption></figure>

spatial ou le pas de temps utilisés. L'importance de la relation de dissipation d'entropie au niveau
continu permet en particulier d'obtenir de nombreuses autres préservations, que nous détaillons
ci-dessous.

**Positivé des solutions discrètes.** On peut considérer qu'obtenir une densité approchée négative est une aberration vis-à-vis de la physique. D'un point de vue pratique, l'exploitation d'une densité non positive est difficile, voire impossible. D'un point de vue théorique, la question de la préservation de la positivité peut aussi être très importante pour l'analyse. Une telle préservation s'obtient au moyen d'une entropie qui est affectée par le signe de la solution, comme l'entropie de Boltzmann.

**Préservation des équilibres thermiques.** La question des équilibres thermiques est elle aussi capitale, au vu du rôle joué par ces états stationnaires particuliers dans les modèles et leur analyse. Définir précisément une notion de schéma préservant les équilibres thermiques n'est pas aisé, et plusieurs définitions seront proposées et utilisées dans cette thèse, en fonction des méthodes considérées. Pour une équation d'advection-diffusion (1), une telle propriété correspond essentiellement au fait que l'état stationnaire discret pour des données compatibles avec l'équilibre est, en un certain sens, une interpolation de l'équilibre thermique continu. Pour les systèmes de dérive-diffusion, ceci correspond au fait que l'état stationnaire discret peut être interprété comme une approximation d'une solution de l'équation de Poisson-Boltzmann (14). En règle générale, l'annulation du terme de dissipation d'entropie discrète permet de caractériser les équilibres discrets.

**Comportement en temps long des solutions discrètes.** Pour les modèles étudiés ici, le comportement asymptotique $t \to \infty$ est très particulier : les solutions convergent exponentiellement vite vers des états stationnaires. Concevoir des schémas dont les solutions vérifient une propriété analogue est donc un gage conséquent de fiabilité. Pour certaines applications qui s'intéressent à des phénomènes évoluant sur des échelles de temps très longues, comme les phénomènes de corrosion [16, 53], la question du temps long est cruciale car elle assure la fiabilité des simulations. Un schéma avec un bon comportement en temps long permet également de calculer des états stationnaires à moindre coût. Si l'on considère l'exemple des semi-conducteurs, on s'intéresse généralement aux courbes courant-tension, qui concernent des états stationnaires. Pour simuler ces dispositifs, on a recours à des schémas (non-linéaires) qui nécessitent d'être initialisés par des données suffisamment proches de la solution. On a alors recours à des procédés de continuation. Une des possibilités consiste à utiliser cette initialisation comme une donnée initiale, et laisser évoluer
(en temps) le problème. Si la solution discrète converge exponentiellement vite vers un état stationnaire discret, on obtient une approximation de l'état stationnaire en très peu d'itérations en temps.

Enfin, pour les modèles où la masse est préservée (par exemple, advection-diffusion avec conditions de bord de Neumann homogènes), il est également important d'assurer cette caractéristique au niveau discret. Comme pour la positivité, cette conservation peut être utilisée à des fins d'analyse théorique des schémas.

## Volumes finis à deux points : forces et faiblesses

Nous présentons maintenant une classe de schémas qui répond à presque tous les critères discutés à la section précédente, massivement utilisée pour simuler des modèles de semi-conducteurs : les méthodes volumes finis à deux points. Ces méthodes sont basées sur une discrétisation du flux de matière. En définissant de manière adaptée les flux numériques, on obtient des schémas très robustes, puisqu'ils s'appuient sur des quantités qui portent un sens physique. Initialement introduites pour des problèmes hyperboliques [188, 150], l'utilisation de ces méthodes s'est ensuite répandue à des problèmes elliptiques et paraboliques [121].

L'idée des schémas volumes finis deux points (Two Points Flux Approximation, TPFA en anglais) pour les problèmes elliptiques est d'introduire, pour une discrétisation $\mathcal{D}$ du domaine donné, des inconnues discrètes

$$ \underline{u}_{\mathcal{D}} = (u_K)_{K \in \mathcal{M}}, $$

où $\mathcal{M}$ est le maillage et $u_K$ l'inconnue de maille. Le problème continu $\mathrm{div}(J(u)) = f$ est alors discrétisé en intégrant l'équation sur chaque maille $K \in \mathcal{M}$, et en utilisant une formule de Green :

$$ \forall K \in \mathcal{M}, \quad \int_{\partial K} J(u) \cdot n = \int_K f. $$

Muni des inconnues de mailles, on approche le flux normal $J(u) \cdot n$ sur la face $\sigma$ par un flux numérique $F_{K,\sigma}(\underline{u}_{\mathcal{D}})$. On obtient finalement le schéma

$$ \forall K \in \mathcal{M}, \quad \sum_{\sigma \in \mathcal{E}_K} F_{K,\sigma}(\underline{u}_{\mathcal{D}}) = \int_K f, $$

où $\mathcal{E}_K$ est l'ensemble des faces de la maille $K$. Dans le cas où le flux est diffusif ($J(u) = \Lambda \nabla u$), on utilise une discrétisation différences finies des gradients. Cette méthode est alors très simple à implémenter, et préserve les propriétés de la solution continue. Elle souffre néanmoins d'une limitation intrinsèque : pour être consistant avec les problèmes continus, il faut se placer sur un maillage dit admissible vérifiant une condition d'orthogonalité (pour le produit scalaire induit par le tenseur de diffusion $\Lambda$) [121, Définition 3.1]. En pratique, même pour des tenseurs anisotropes relativement simples, il est très compliqué de construire de tels maillages admissibles. Lorsque l'on considère des problèmes plus compliqués (diffusion anisotrope et hétérogène, plusieurs tenseurs de diffusion pour un système couplé), de tels maillages peuvent tout simplement ne pas exister. Ainsi, l'usage de ces schémas est essentiellement restreint à des cas isotropes, sur des maillages avec des géométries simples (et sans nœuds orphelins).

Les schémas TPFA sont l'objet d'une abondante littérature. Pour approcher l'équation d'advection-diffusion linéaire (1) (avec $\Lambda = I_d$), le flux de Scharfetter-Gummel [239, 160, 73] est une solution privilégiée : il permet de gérer des advections très fortes (ce flux peut être interprété comme une interpolation d'un flux parabolique et d'un flux hyperbolique, obtenue grâce à la fonction de Bernoulli $B(s) = \frac{s}{e^s-1}$), préserve la positivité et les équilibres thermiques au niveau
discret. Il présente également un bon comportement en temps long, tout comme certaines de ses généralisations [131, 69].

Concernant les modèles de dérive-diffusion isotropes, on peut également citer les travaux fondateurs de [239]. Suivant cette approche, de nombreuses recherches autour de schémas TPFA positifs préservant le comportement en temps longs ont été menées [68, 27, 28]. Notons également que la question de généraliser le flux de Scharfetter-Gummel à des statistiques non-Boltzmann a été un sujet actif [265, 120, 26, 140, 50]. Toutes ces méthodes préservent la positivité et l'équilibre thermique. De nombreux travaux comparatifs (numériques et théoriques) [125, 127, 50] ont également permis d'identifier les atouts et faiblesses de chacune de ces méthodes.

Malheureusement, tous les résultats énoncés plus haut ne sont valables que pour des problèmes isotropes sur maillages orthogonaux. Ainsi, il n'existe pas de schéma satisfaisant pour le modèle anisotrope (3). La seule tentative connue est le schéma TPFA de [143] pour la statistique de Boltzmann ($h = \log$), qui repose sur une perturbation des flux de Scharfetter-Gummel pour gérer l'anisotropie. Les auteurs parviennent à montrer que la positivité des solutions est préservée sous certaines conditions liant la géométrie du maillage avec l'intensité du champ magnétique. En pratique, ces conditions sont extrêmement restrictives et difficiles à vérifier : les résultats numériques présentés indiquent clairement une perte de la positivité au niveau discret.

## Schémas pour problèmes anisotropes sur maillages généraux

Afin de gérer l'anisotropie et régler le défaut majeur des schémas TPFA, de nombreuses méthodes de type volumes finis ont été développées depuis les vingt dernières années [157, 105, 40, 39, 110, 123, 124] pour discrétiser des problèmes de diffusion anisotrope. Par rapport aux schémas TPFA, ces méthodes reposent sur l'ajout d'inconnues auxiliaires qui permettent de gérer l'anisotropie de la diffusion et la généricité du maillage, via la reconstruction d'un opérateur gradient discret. Elles peuvent cependant s'exprimer grâce à des flux numériques, qui ne sont plus basés sur une formulation différences finies entre deux points. La plupart de ces méthodes peuvent s'inscrire dans le formalisme des "schémas gradients" [111]. Ces méthodes permettent d'approcher les problèmes anisotropes avec une grande robustesse, mais celle-ci s'obtient au prix d'une perte de la positivité des solutions numériques, comme obervé lors du benchmark de la conférence FVCA 5 [154].

L'utilisation de ces méthodes a ensuite été généralisée à des problèmes de convection-diffusion qui nous intéressent dans ce manuscrit. Pour le problème d'advection-diffusion linéaire (1), un schéma Hybride Mimétique Mixte (HMM) a été proposé et analysé [19]. Ce schéma, étudié dans le Chapitre 1, est conçu comme une généralisation des B-schémas TPFA [66]. Cependant, il ne préserve ni la positivité ni les équilibres thermiques. Concernant les modèles de semi-conducteurs isotropes, des schémas Volumes Finis en Dualité Discrète (Discrete Duality Finite Volume, DDFV en anglais) on été proposés [64] pour supporter les maillages généraux sur des domaines bi-dimensionnels. On peut aussi citer le schéma HMM [81] pour un modèle de plasma (qui est essentiellement un système de dérive diffusion en statistique de Boltzmann avec une anisotropie sur l'une des deux espèces). Ces schémas ne préservent pas la positivité des densités, ce qui restreint leur analyse.

Afin de palier le défaut de positivité des schémas précédemment cités, différentes stratégies ont été proposées. Elles reposent sur l'introduction d'une nonlinéarité dans le schéma. On peut notamment citer les méthodes [245, 113, 32, 48] qui reposent sur l'ajout d'une nonlinéarité dans le flux numérique. Ces méthodes s'avèrent robustes vis-à-vis de l'anisotropie et des maillages pour des problèmes de diffusion linéaire, et l'existence et la consistance des flux sont prouvées [32]. Plus récemment, une méthode pour des systèmes de Poisson-Nernst-Planck (un
modèle pour décrire l'évolution d'ions dans une solution, semblable au système (2) mais avec des conditions de bord de type Neumann sur tout le bord du domaine) sur maillage généraux a été introduit et partiellement analysé [250] en suivant cette méthodologie.

Une autre solution a été proposée par Cancès et Guichard [55, 56] pour une classe d'équations non-linéaires contenant l'équation d'advection-diffusion linéaire : discrétiser le flux d'advection-diffusion

$$J(u) = \Lambda (\nabla u + u \nabla \phi)$$

sous la forme non-linéaire suivante

$$J(u) = u \Lambda \nabla (\log(u) + \phi).$$

Par rapport aux méthodes présentées plus haut, l'avantage de cette stratégie est qu'elle repose sur une non-linéarité qui a un sens au niveau continu. Celle-ci donne notamment au schéma une structure d'entropie (de Boltzmann), qui permet de proposer une analyse (existence de solutions positives et convergence). De plus, il est assez simple de généraliser ce schéma à des diffusions non-linéaires (correspondant à des modèles avec des statistiques non-Boltzmann). Suivant ce principe, un schéma non-linéaire Volumes Finis en Dualité Discrète a été dévéloppé et analysé [52, 51]. Grâce à des inégalités fonctionnelles de type Sobolev logarithmique, une analyse du comportement en temps long reposant sur l'entropie de Boltzmann est menée.

## Méthodes d'ordre élevé

Nous nous intéressons maintenant aux méthodes d'ordre élevé. L'idée générale de ces méthodes est d'utiliser des inconnues discrètes polynomiales pour obtenir une précision accrue, pour peu que les solutions continues soient assez régulières. Ce gain de précision est généralement accompagné d'un gain d'efficacité : pour un coût de calcul donné, on pourra atteindre une meilleure précision avec une méthode d'ordre élevée qu'avec une méthode d'ordre bas.

Bien qu'il existe de nombreuses méthodes d'ordre élevé pour les problèmes diffusifs, dont les éléments finis $\mathbb{P}^k$, nous nous concentrons ici sur des méthodes conservatives, qui peuvent s'interpréter comme des généralisations des méthodes d'ordre bas présentées précédemment pour des problèmes de diffusion linéaires (pouvant être anisotropes et posés sur maillages généraux). Historiquement, les premières méthodes introduites sont les méthodes Galerkin Discontinues [13, 9, 10] (Discontinuous Galerkin, DG en anglais). On peut interpréter ces schémas comme une généralisation à l'ordre élevé des méthodes volumes finis TPFA : le schéma repose sur un équivalent discret du flux entre chaque maille. Cette caractéristique permet d'assurer le caractère conservatif de la méthode. Pour en améliorer l'efficacité, des versions hybrides ont été introduites [85, 87]. Ces méthodes de Galerkin Discontinue Hybrides (Hybridisable Discontinuous Galerkin, HDG en anglais) reposent sur l'ajout d'inconnues sur les faces afin de découpler les inconnues de mailles. Ceci permet de réaliser une élimination des inconnues de mailles (condensation statique) et ainsi obtenir un gain d'efficacité conséquent pour des ordres et dimensions élevés. Une approche similaire, en utilisant des stabilisations différentes, donne lieu aux méthodes Hybrides d'Ordre Elevé [100, 99, 97, 83] (Hybrid High Order, HHO en anglais). Elles permettent notamment d'obtenir un ordre de convergence optimal ($k+2$ en norme $L^2$ pour des inconnues de degré $k$), contrairement aux méthodes HDG originelles. Enfin, on peut citer les méthodes des Elements Virtuels [20, 21, 22, 7] (Virtual Element Method, VEM en anglais), dont l'analyse repose sur l'introduction d'espaces de fonctions non-polynomiales.

Les trois méthodes précédemment citées (HDG, VEM et HHO) sont fortement liées les unes aux autres, et partagent une caractéristique commune : des inconnues auxiliaires, placées sur le squelette du maillages (réunion des faces) sont utilisées pour réaliser une condensation statique,
qui permet d'éliminer les inconnues de maille. On parle parfois de méthodes squelettales. Les liens entre les méthodes HDG et HHO ont été étudiés et explicités dans [84], tandis que les similarités entre les méthodes VEM et HHO sont discutées dans [192].

L'utilisation de ces méthodes squelettales d'ordre élevé pour des problèmes de convection-diffusion est l'objet de travaux depuis plus d'une dizaine d'années. Concernant les équations d'advection-diffusion linéaires, de nombreuses méthodes HDG ont été développées [86, 78, 79, 139, 229]. On peut également citer la méthode HHO [98], qui est une extension à l'ordre élevé du schéma HMM [19]. Pour l'approximation de systèmes de dérive-diffusion, divers schémas DG et HDG ont aussi été proposés [200, 77]. Toutes ces méthodes, qui traitent la diffusion de manière linéaire, présentent un défaut de positivité.

A notre connaissance, la conception de schémas d'ordre élevé préservant la positivité est un travail encore assez peu exploré, qui est essentiellement restreint à des méthodes DG.

La première stratégie connue pour préserver la positivité repose sur l'utilisation de schémas DG dont la solution est positive en moyenne sur les mailles [198, 196]. Les inconnues peuvent prendre des valeurs négatives dans la maille, mais on s'assure que leurs moyennes sont positives en choisissant habilement des paramètres de stabilisation. Ce choix de stabilisation est délicat, et dépend de l'ordre de la méthode ainsi que de la dimension du problème. De plus, le fait que la solution puisse être ponctuellement négative limite pour le moment l'analyse des schémas. Cette méthode a aussi été appliquée à des systèmes de Poisson-Nernst-Planck [197]. Une autre approche est proposée dans [36], où un schéma DG pour une équation de Fisher-KPP (équation de la forme $∂_t u - Δu = u(1-u)$, avec des conditions de Neumann homogènes) est introduit et analysé. Ce schéma suit les idées des schémas non-linéaires de [55, 56], et est basé sur la structure d'entropie du problème, utilisée avec un changement d'inconnue non-linéaire qui permet d'assurer la positivité de la solution. La positivité et la structure d'entropie discrète permettent alors de fournir une analyse étendue du schéma, incluant l'existence de solutions, le comportement en temps long et la convergence des solutions vers un problème semi-discrétisé en temps. L'analyse du schéma est réalisée pour un ordre et une dimension arbitraire. Cette analyse repose sur un terme de stabilisation bien choisi, qui fait intervenir les normes $L^∞$ des inconnues discrètes sur les faces du maillage. Les résultats numériques présentés pour un problème unidimensionnel indiquent bien une préservation de la positivité et un bon comportement en temps long, mais la question de l'ordre de convergence du schéma n'est pas discutée.

## Méthodes étudiées dans cette thèse : HFV et HHO

Le positionnement initial de cette thèse est d'utiliser les stratégies introduites dans [55, 56, 52, 51] pour développer et analyser des méthodes pour le modèle anisotrope (3) permettant d'assurer la positivité des densités calculées. Parmi les méthodes volumes finis développées pour gérer l'anisotropie de problèmes linéaires, une classe de schémas, les méthodes HMM, se distingue des autres de par l'existence de généralisations d'ordre élevé. Puisque l'utilisation de méthodes d'ordre bas est un succès, la question de savoir s'il est possible de proposer et d'analyser des méthodes non-linéaires d'ordre élevé suivant les principes de [55, 56, 52, 51] se pose naturellement. Avant de pouvoir développer des méthodes d'ordre élevé, un travail préliminaire consiste à adapter les stratégies non-linéaires au cas des versions d'ordre bas des schémas. Les méthodes qui nous intéressent dans cette thèse sont dites hybrides, basées sur des inconnues dans les mailles et sur les faces.

Afin de concevoir des schémas préservant la structure, nous avons besoin d'utiliser des méthodes implicites (en temps) et généralement non-linéaires. Une des questions crucales concernant leur analyse est l'existence de solutions. Comme pour l'étude des EDP il n'existe
pas de méthodologie générale pour ce genre de résultats. La technique utilisée dans cette thèse repose sur une méthode dite de continuité : l'idée est de comparer le schéma étudié avec un schéma plus simple, dont on sait qu'il admet des solutions (typiquement, un schéma linéaire). Cette comparaison se fait en reliant de manière continue les deux schémas. Il faut alors obtenir des estimations a priori sur les solutions, qui permettent d'en déduire l'existence de solution par des arguments de type degré topologique (voir par exemple [107] pour une introduction au degré ainsi que des applications à des résultats d'existence pour des problèmes continus et discrets). Ainsi, l'enjeu majeur de l'étude de schémas non-linéaires est l'obtention d'estimations a priori. Il est intéressant de noter que, contrairement au cas de schémas linéaires, la question de l'unicité des solutions ne découle en général pas des résultats d'existence. Cependant, cette problématique ne se limite pas à des enjeux théoriques : la qualité d'une méthode doit aussi être questionnée en pratique, à travers son implémentation dans un code de simulation. La question de l'existence est aussi délicate en pratique qu'en théorie : puisque calculer une solution revient à résoudre un système non-linéaire, comment peut-on la déterminer en pratique ? Dans notre cas, l'implémentation repose sur des méthodes de Newton, pour lesquelles nous n'avons pas de garantie théorique de convergence globale. Par conséquent, un enjeu réel est de vérifier que l'on arrive bien à calculer une solution. L'implémentation d'un schéma non-linéaire qui permet effectivement d'obtenir des solutions approchées est donc un résultat conséquent. Afin de s'assurer que le schéma (et son implémentation) est robuste, il est également nécessaire de réaliser des tests dans différentes situations (maillages, données physiques). Cette étape de test doit aussi permettre de confirmer que les propriétés théoriques du schéma sont présentes en pratique, notamment pour les schémas préservant la structure (positivité, comportement en temps long).

Notons finalement qu'il existe en général tout un panel de méthodes pour approcher un problème donné. Il est donc important de proposer une comparaison numérique, qui permet d'identifier les points forts des méthodes et d'en comprendre les cadres d'utilisation les plus appropriés. Cet aspect comparatif est assez présent dans cette thèse (Chapitres 1, 3 et 4). En particulier nous proposons lorsque c'est possible un travail de comparaison entre des méthodes linéaires et des méthodes non-linéaires. Il est important de noter que de tels travaux sont délicats, car les résultats peuvent dépendre de l'implémentation des schémas plus que des schémas eux même.

Les méthodes d'ordre bas étudiées dans cette thèse sont basées sur les Volumes Finis Hybrides [123] (Hybrid Finite Volume, HFV en anglais, qui sont un cas particulier du schéma SUSHI), qui coïncident avec les versions d'ordre bas des méthodes Hybrides d'Ordre Élevé [100] (Hybrid High Order, HHO en anglais). Notons que les méthodes HFV partagent des liens très forts avec les méthodes Différences Finies Mimétiques [40, 39] et les méthodes Volumes Finis Mixtes [110]. Dans [112], les auteurs unifient l'analyse et la description de ces trois méthodes, sous le nom de méthode Hybride Mimétiques Mixte (HMM).

Nous donnons ici une description succincte de la philosophie des schéma HFV. Pour une discrétisation $\mathcal{D}$ du domaine donné, on introduit des inconnues discrètes

$$ \underline{u}_D = ((u_K)_{K \in \mathcal{M}}, (u_\sigma)_{\sigma \in \mathcal{E}}) $$

constituées d'inconnues de maille $u_K \in \mathbb{R}$ d'inconnues de face $u_\sigma \in \mathbb{R}$. L'approche de la méthode est locale : pour chaque maille $K$, on va définir un gradient discret sur $K$ qui dépend des inconnues locales : l'inconnue de maille $u_K$ et les inconnues de face de la maille $(u_\sigma)_{\sigma \in \mathcal{E}_K}$. Cette approche permet d'obtenir un gradient discret $\nabla_D \underline{u}_D$ sur $\Omega$, qui est utilisé comme analogue discret du gradient continu. De manière équivalente, on peut définir des flux numériques $F_{K,\sigma}(\underline{u}_D)$ associés à ce gradient discret, et écrire la méthode comme une méthode volumes finis. Comme pour le gradient, les flux sont locaux, dans le sens où $F_{K,\sigma}(\underline{u}_D)$ dépend uniquement de
$u_K$ et des inconnues de face locales $(u_\sigma)_{\sigma \in \mathcal{E}_K}$ :

$$F_{K,\sigma}(\underline{u}_D) = \sum_{\sigma' \in \mathcal{E}_K} A_K^{\sigma\sigma'}(u_K - u_{\sigma'}),$$

où les $A_K^{\sigma\sigma'}$ dépendent de la géométrie du maillage et du tenseur. L'une des forces majeures de cette approche locale consiste en la possibilité de réaliser une condensation statique lors de la résolution numérique des systèmes linéaires correspondant : l'inconnue de maille $K$ peut être exprimée en fonction des inconnues de face associées. Il est alors possible d'éliminer localement les inconnues de maille, pour se ramener à un système linéaire sur les inconnues de face uniquement. Cette procédure offre un gain d'efficacité conséquent, en diminuant de fait le nombre d'inconnues de la méthode. Enfin, notons que les schémas HFV coïncident avec les schémas TPFA dans des cas isotropes sur des maillages fortement contraints, dits suradmissibles [123].

Comme pour les méthodes HFV, les schémas HHO sont basés sur des inconnues hybrides

$$\underline{u}_D = ((u_K)_{K \in \mathcal{M}}, (u_\sigma)_{\sigma \in \mathcal{E}}),$$

mais ici les inconnues de maille et les inconnues de face sont polynomiales. Comme pour la méthode HFV, l'idée principale est de reconstruire un gradient sur chaque maille, en utilisant toutes les inconnues locales à disposition (inconnue de maille et inconnues de face). La définition de ce gradient suit une formule d'intégration par parties discrète, qui est testée contre des fonctions polynomiales. Ce gradient discret est consistant avec le gradient continu, mais il présente un défaut de stabilité. Pour obtenir un schéma stable, il faut également prendre en compte les sauts entre les inconnues de mailles et les inconnues de face, via un terme de stabilisation. Le choix de cette stabilisation est crucial, car il a une influence directe sur l'ordre de convergence. La stabilisation originelle de [100] permet, pour des inconnues de maille et de face de degré $k$ d'obtenir un schéma d'ordre optimal $k+2$ en norme $L^2$. Avec ce choix, le schéma HHO pour $k=0$ coïncide avec la méthode HFV sous certaines hypothèses sur les données et le maillage [100]. Notons que la stabilisation initialement utilisée pour les méthodes HDG permet d'obtenir une convergence d'ordre $k+1$ en norme $L^2$ pour la même configuration.
La stabilisation utilisée au Chapitre 4 est différente de la stabilisation HHO originelle. Elle repose sur des inconnues de maille "enrichies" de degré $k+1$ et des inconnues de face de degré $k$. Avec de telles inconnues discrètes, la définition du gradient discret ne change pas, mais on peut utiliser une stabilisation dite de Lehrenfeld-Schöberl [190, 191], qui permet d'obtenir un schéma d'ordre $k+2$ en norme $L^2$. Dans le même temps, la condensation statique permet de garder un système global qui ne dépend que des inconnues de face, de degré $k$. De telles méthodes sont nommées méthodes HHO d'ordre mixte $(k,k+1)$ ou méthodes HDG+. L'intérêt de la stabilisation de Lehrenfeld-Schöberl est qu'elle est bien plus simple à exprimer, manipuler et implémenter que la stabilisation originelle de [100].

## Overview of the works and results

In this section, we describe the main theoretical and numerical results presented in the manuscript. For each chapter, we also give some details about the codes implemented.
Long-time behaviour of Hybrid Finite Volume schemes for advection-diffusion (Chapter 1)

In the first chapter of this thesis, we focus on the numerical approximation of the linear advection-diffusion toy model

$$ \partial_t u - \operatorname{div}(\Lambda(\nabla u + u\nabla\phi)) = f $$

with Hybrid Finite Volume (HFV) methods preserving entropy structures of the continuous problem. We focus on HFV methods because they are a first step towards HHO schemes.

We first introduce two linear methods (namely the one introduced in [19] and an exponential fitting scheme based on [43], relying on the change of unknown $\rho = u e^{\phi}$). Then, we develop a new nonlinear scheme based on the reformulation of the flux $u\Lambda\nabla(\log(u) + \phi)$ as in [56, 52] for Vertex Approximate Gradient (VAG) and Discrete Duality Finite Volume (DDFV) methods. The three HFV schemes have an entropy structure, related to quadratic entropies for the linear schemes and Boltzmann entropy for the nonlinear one. Thanks to these entropy dissipation estimates, we show several theoretical results.

The first theoretical result of this chapter is related to the existence of (positive) solutions to the nonlinear scheme, stated in Theorem 1. The associated proof differs from the existing ones of [56, 52, 51, 57], and is based on a discrete entropy method alongside with a Brouwer fixed-point argument. The method of proof is generic and can be easily adapted to any schemes with an entropy structure, see for instance [156, 1] and the other Chapters of the present manuscript. Other analysis results relate to the discrete long-time behaviour of the three schemes. These results are stated in Theorems 2, 3 and 4, and assert that the discrete solutions converge exponentially fast towards their associated steady-states as time tends to infinity. The proofs essentially rely on discrete functional inequalities. For both linear schemes, an analysis for data out of equilibrium is available.

As an intermediate step of the long-time behaviour analysis, we have established a Logarithmic-Sobolev inequality adapted to problems with Dirichlet boundary conditions in Appendices 1.A.2 and 1.A.3. Up to our knowledge, these results are new, both in the continuous and discrete frameworks.

We numerically validate our theoretical findings on a set of test-cases (see Figure 2) and, for completeness, we also compare the accuracy of the three schemes on stationary problems. These tests reveal that the two linear schemes can indeed compute non-positive solutions in practice. Moreover, all three schemes are accurate at order two in space for $L^2$ norm.

This work is an article published in Numerische Mathematik [70], in collaboration with Claire Chainais-Hillairet, Maxime Herda and Simon Lemaire.

About the codes

The HFV codes related to this first work are implemented in C++. The codes rely on an already existing mesh class for 2 dimensional meshes. The implementation is based on the "Finite Volume" flavour, using numerical fluxes. The main difficulty of these codes lies in the nonlinear schemes and the related Newton methods. Two codes were developed for the need of Chapter 1:

(i) **hfv_evol**: a code for linear anisotropic advection-diffusion equations with homogeneous Neumann boundary condition, zero source term and symmetric diffusion tensor. It includes
<figure><img src="image_10.png" /><figcaption>Figure 2 – Long-time behaviour of discrete solutions to HFV schemes on a Kershaw mesh.</figcaption></figure>

the linear HMM scheme (1.26) of [19], the linear HFV exponential fitting scheme (1.35) and the nonlinear HFV scheme (1.45). The HMM scheme is implemented to handle general (not necessarily irrotational) advection fields. This code was used for the numerical results about positivity and long-time behaviour presented in Sections 1.5.2 and 1.5.3.

(ii) **hfv_sta:** a code for stationary linear anisotropic advection-diffusion equations with mixed Dirichlet-Neumann boundary condition (with a Dirichlet boundary of positive measure) with general boundary data, general source term and symmetric diffusion tensor. As for *hfv_evol*, this code includes the three schemes presented in Chapter 1. For the nonlinear scheme, in order to ensure the convergence of the Newton method, a continuation procedure from the thermal equilibrium towards the real problem was implemented. This code was used for the accuracy study of Section 1.5.4.

**Structure-preserving schemes for semiconductor models (Chapters 2 and 3)**

In Chapters 2 and 3, we assert the relevance of the previously selected method (the nonlinear HFV scheme) by using it on a complex coupled system of PDEs (3) modelling semiconductor devices. The aim of this work is threefold:

- use the nonlinear HFV scheme on a nonlinear and anisotropic problem which cannot be approximated thanks to TPFA schemes;

- give us a better understanding of the theoretical features needed to approximate this kind of problems (with the aim of developing high-order schemes);

- check that such schemes actually work in practice for more complex problems than linear advection-diffusion equations.

In Chapter 2, we are interested in the numerical approximation of a drift-diffusion system which models semiconductor devices. The model under consideration is a generalisation of the one introduced in [143]. Its main feature is the fact that it takes into account an external magnetic field, which induces a rotation of the charge carriers (see Figure 3) and therefore
anisotropy on the convection-diffusion equations. From a mathematical point of view, such a system is modelled by the following system, where the unknowns N and P are the density of negative (electrons) and positive (holes) charges, and $\phi$ is the electrostatic potential induced by the spatial inhomogeneity of the electric charges:

$$ 
\left\{ 
\begin{array}{l}
\partial_t N - \mathrm{div}(N \Lambda_N \nabla(h(N) - \phi)) = -R(N, P) \\
\partial_t P - \mathrm{div}(P \Lambda_p \nabla(h(P) + \phi)) = -R(N, P) \\
-\mathrm{div}(\Lambda_\phi \nabla \phi) = C + P - N.
\end{array}
\right. 
$$

In this system, $\Lambda_N$ and $\Lambda_p$ are (potentially) nonsymmetric anisotropic tensors related to the magnetic field, $h$ is a nonlinearity in the diffusion (related to the so-called statistics used) and $R(N, P)$ is a reaction term. Because of the anisotropy induced by the magnetic field, it is not possible to design TPFA-admissible mesh for this problem. In [143] Gajewski and Gärtner propose a scheme based on a perturbation of the TPFA Scharfetter-Gummel scheme in order to handle this anisotropy for Boltzmann statistics ($h = \log$), but the scheme proved to compute negative densities on some situations. On the other hand, the positivity of the density is crucial both from a practical (since negative densities holds no physical meaning) and theoretical (since the analysis of the schemes relies on entropy structure, which do not have any meaning if the densities are negative) point of view. Thus, we use the nonlinear HFV method introduced in Chapter 1 in order to tackle both anisotropy and positivity issues, and introduce a structure-preserving scheme for this model.

With respect to the Chapter 1, the main novelty and difficulty of this work is the support of the coupling between unknowns. Another difference comes from the nonsymmetry of the tensors $\Lambda_N$ and $\Lambda_p$, but it proves not to be a significant issue for the HFV method used. Last, since the diffusion considered for the equation on the charge carriers are nonlinear, one has to handle general nonlinearities and potentially two bounds on the densities. The analysis highlights the importance of the quasi-Fermi potentials $w_N = h(N) - \phi$ and $w_P = h(N) + \phi$, which are the key variables to handle the coupling. Thanks to this approach, we show that the scheme admits solution with relevant physical bounds on the densities (which are in particular positive), and preserve the long-time behaviour of the solution.

Our first result regarding the approximation of semiconductor models is Theorem 7, in which we state an existence result for the scheme. The computed densities are shown to be bounded in appropriate intervals (and in particular are positive) thanks to the use of a discrete entropy structure. We use the strategy employed in Chapter 1 together with a key correspondence discussed in Section 2.3.2 between discrete densities and discrete quasi-Fermi potentials through the nonlinear stationary scheme (2.42).

A second result of this work is Theorem 9, in which we establish the exponential decay of the discrete solutions towards the associated thermal equilibrium on a fixed mesh. The proof of this result relies on the use of a discrete Poincaré inequality on the discrete quasi-Fermi potentials.

From a numerical point of view, we give some evidence of our theoretical results, and investigate the relation between the cost of the scheme (number of Newton iterations) and the intensity of the magnetic field. Moreover, we assert the numerical robustness of the scheme with respect to the physical data, including the intensity of the magnetic field.

Following this work, we focus in Chapter 3 on the numerical comparison of two different structure-preserving methods for drift-diffusion systems, namely the HFV scheme introduced in Chapter 2 and a DDFV scheme based on the scheme for advection-diffusion [52, 51]. Up to our knowledge, such a comparative work between positivity preserving schemes on general meshes for coupled problems is new, but is fundamental in order to to understand the practical advantages and limitations of each method. For the sake of simplicity and conciseness, we
<figure><img src="image_11.png" /><figcaption>Figure 3 – Influence of the magnetic field intensity on the density of electrons of a non-equilibrium steady state.</figcaption></figure>

consider a simplified model with Boltzmann statistics, without anisotropy neither reaction terms. We present both schemes in a similar way, in order to point out the numerous theoretical similarities between both frameworks. Therefore, the notations used in this Chapter (for the mesh and the HFV scheme) differ a little bit from those used in the rest of the manuscript. We then perform a numerical comparison about the long-time behaviour and the preservation of the positivity. Overall, both schemes present a similar behaviour, both from a quantitative and qualitative point of view.

These works are published in M2AN [213] (Chapter 2) and in the book of proceedings of the conference FVCA X [187] (Chapter 3, in collaboration with Stella Krell).

## About the codes

The HFV code related to this part, `hfv_SC`, is implemented in C++ and follows the structure of `hfv_evol`. It is a code for transient anisotropic drift-diffusion system with mixed Dirichlet-Neumann boundary conditions. This code is developed to handle both Boltzmann and Blakemore statistics, realistic recombination/generation term (either Auger recombination, Shokley-Read-Hall term of a linear combination of both), anisotropic and non-symmetric tensors. With respect to the previous code, the main difficulty of the implementation lies in the coupling between the three different unknowns (positive and negative charges carrier density, electrostatic potential). Considering nonlinear diffusion (Blakemore statistics) for the densities also implies new challenges, including the fact that one has to handle two bounds on the densities (a lower one, which is 0 and upper one) in the Newton method. Another main difference lies in the stiff behaviour induced by the coupling (for small Debye lengths and strong magnetic fields), which requires a more refined choice of parameters (thresholds of truncation and stopping criteria). This code was used for the numerical results of Chapters 2 and 3.

Regarding the DDFV scheme studied in Chapter 3, it was implemented by Stella Krell.
High-order schemes for advection-diffusion (Chapter 4)

In the last chapter, we are interested in the discretisation of advection-diffusion equations
with high-order methods in space. The main goal of this part is to develop and analyse a high-
order (which means, at least of order 3 in L² norm in space) scheme which keeps the structure
preserving features of the low-order schemes introduced in the previous chapters, and especially
preserves the positivity. There is moreover a real need to investigate the numerical behaviour
and the efficiency of such schemes for 2D problems. Thus, our objective is twofold:

- develop and analyse high-order structure-preserving schemes for diffusive problems;

- exhibit reasonable evidence that these high-order methods are appealing in terms of computational efficiency (compared to low-order structure preserving methods).

As in Chapter 1, this work focuses on a single linear advection-diffusion equation, but is of
course intended to be used on more complex problems, including the anisotropic drift-diffusion
system studied in Chapter 2.

We now provide a brief discussion about the locks and the key we had to identify in order
to achieve the main goal of the thesis. The difficulties one needs to face in the design of a
positivity-preserving high-order scheme inspired from the HFV scheme of Chapter 1 essentially
lie in the following questions.

(i) **How can we ensure (both from a theoretical and a practical point of view) that the polynomial unknowns stay positive?**

Indeed, it is not easy a priori to infer the sign of a polynomial function (of arbitrary degree) from the degrees of freedom (values of the function at some points, or coefficients of the function in a given polynomial basis) which are the quantities available in practice in the code. Such a difficulty seems to be further reinforced by general meshes.

(ii) **How can we define high-order counterparts of the "reconstruction operators" used for the nonlinear HFV schemes?**

Indeed, the design and the analysis of the nonlinear HFV schemes rely on the definition of a local operator $r_K : \underline{V}_K \to \mathbb{R}$, which maps the local unknowns to a constant within the cell. A natural high-order generalisation of this operator (in order to get a consistent scheme at the expected order) should be an operator $r_K^k : \underline{V}_K^k \to \mathbb{P}^k(K)$, where $\underline{V}_K^k$ is the set of local (polynomial) unknowns, and $\mathbb{P}^k(K)$ is the space of polynomials functions on the cell $K$. However, using a simple mean between cell and face unknowns, as it was done for the low-order schemes, is not a relevant possibility anymore since the face unknowns are not defined on the cell.

(iii) **How can we generalise the method to obtain the $L^\infty$ a priori estimates on the discrete solution?**

Indeed, the existence results for the low-order schemes rely on entropy methods: a control
on the discrete entropy and its dissipation ensures that the discrete solution is bounded in
$L^\infty$ (with meshsize-dependant bounds). Up to now (and in similar works with different
space discretisations [56, 52, 57]) the proofs of such results are based on a propagation of
the information: one has a bound on a given cell unknown thanks to mass conservation, can
therefore compare this unknown with the surrounding face unknowns and get estimates on
these face unknowns. Then, proceeding in the same way, we compare the face unknowns
with nearly cell unknowns, and get a bound for every unknown. For high-order schemes, a
natural comparison can be performed on faces, between the trace of the cell unknown and
the surrounding face unknowns. However, once we have a bound on a given face unknown,
inferring a bound on the neighbouring cell unknown seems difficult, since we can only get
information on the trace of this unknown.
Of course, these issues are somehow interconnected: the reconstruction operators should also be positive on the cells, and should be the key to get the algebraic structure allowing to get $L^\infty$ bounds on the unknowns. The work presented in Chapter 4 proposes a first solution in order to answer the three previous points. We detail below the key ideas.

(1) **Potentials as polynomials and densities as exponentials of polynomials.** With regards to the analysis performed in Chapter 2, the key idea to solve the positivity issue (i) is to use the quasi-Fermi potential $w = \log(u) + \phi$ as the main unknown. Therefore, we want to approximate $\ell = \log(u)$ as a polynomial function (on each cell and face), and reconstruct a discrete density by taking the exponential of the discrete $\ell$. The resulting discrete solution will thus be positive on the cells and faces, as the exponential of some real functions. Note that this process is easily generalisable to nonlinear diffusion since one just has to use the distribution function $g$ of this statistics instead of the exponential.

(2) **Two different reconstructions for the consistent and stabilisation parts.** To address the reconstruction issue (ii), the solution chosen is to split the reconstruction operator into two different parts: a “consistent” reconstruction defined on the cells, which only takes into account the cell unknowns, and a “stabilising” reconstruction defined on the faces, which takes into account the face unknowns as well as the traces of the cell unknown. In order to simplify the design of the stabilising reconstruction, we choose to use a variant of the classical HHO scheme: instead of using cell and face unknowns of the same (maximal) degree $k$, we increase the maximal degree of the cell unknowns to $k+1$. Therefore, the scheme that we use is based on a mixed-order HHO$(k, k+1)$ scheme (see [83, Section 3.2.1]), in $HHO(l, p) l$ is the degree of face unknowns and $p$ is the degree of cell unknowns), and on a stabilisation whose characteristics are derived from the so-called Lehrenfeld-Schöberl stabilisation [190, 191] in the context of HDG methods (which are sometimes called HDG+ schemes). This strategy allows for a scheme of order $k+2$ (in $L^2$ norm) while the systems to solve depend only on face unknowns, of degree $k$.

(3) **Over-dissipation thanks to over-stabilisation.** To tackle the issue (iii), we first have to understand how such results are achieved for low-order schemes. For all the nonlinear schemes presented and analysed in [56, 52, 57, 70], the proofs of $L^\infty$ a priori estimates rely on the definition of the reconstruction operators, which ensure an over-dissipation: at the discrete level, one adds something to increase the dissipation (for the HFV scheme, we add the faces unknowns to get this over-dissipation, for the conforming $P^1$ scheme of [57], the authors use the fact that the mean of a positive affine function on a simplex controls every local DOF of the unknown). In order to get analogous results in the high-order framework, we introduce a new term in the discretisation which induces such an over-dissipation. This term, which corresponds to a discretisation of $-\epsilon \operatorname{div}(\Lambda \nabla(\log(u) + \phi))$ can be seen as an over-stabilisation, where the $\epsilon$ depends on the mesh and tends to zero as the mesh is refined. Such a term provides sufficient control on the quasi-Fermi potential to obtain discrete $L^\infty$ bounds from the entropy and its dissipation. Moreover, we scale this stabilisation parameter $\epsilon$ with respect to the size of the mesh so that it does not change the order of accuracy of the scheme. On the other hand, this stabilisation is designed to not disturb the mass preservation neither the thermodynamical consistency: the scheme has an entropy structure, and preserves the thermal equilibrium whether with or without this stabilisation. Up to now, we do not know if this over-stabilisation is really needed in the analysis, but numerical results indicates that the use of this term in the scheme does not have significant impact on its behaviour neither on its accuracy.

In order to demonstrate the interest of using high-order schemes, we introduce two arbitrary order skeletal schemes : a linear one (4.14), based on the exponential fitting strategy already
used in Chapter 1, and a nonlinear one (4.20), devised to preserve the positivity according to the previous discussion. Both schemes preserve the thermal equilibrium and enjoy discrete entropy structures. The exponential fitting schemes is essentially introduced for comparison purposes with respect to the nonlinear scheme.

In terms of analysis, our result is the existence of solutions (with positive densities) to the nonlinear scheme (4.20), stated in Theorem 10. This result is based on a discrete a priori estimate stated in Lemma 9, alongside with the entropy structure and the mass conservation property of the scheme. As previously discussed, the use of an over-dissipation term is the key ingredient to get these estimates.

We also state long-time behaviour results for both schemes in Propositions 14 and 16. Regarding the exponential fitting scheme, we obtained the exponential convergence of the solution towards the thermal equilibrium by the use an adapted discrete Poincaré inequality. For the nonlinear scheme, an analysis relying on discrete logarithmic Sobolev inequalities is not yet available, and the convergence of solution in time is based on a compactness argument alongside with the fact that the dissipation tends to vanish.

At the numerical level, we investigate the behaviour of both linear and nonlinear schemes, using the same test-cases as in Chapter 1, and compare some of the results with the low-order schemes of Chapter 1. In particular, we show that our theoretical results about long-time behaviour are satisfied for both scheme, and that the non-linear scheme preserves the positivity of the solution, while exhibiting an optimal order of accuracy (convergence of order $k + 2$ for face unknowns of degree $k$). Moreover, this comparative work indicates two important outcomes. First, the use of the high-order version of the nonlinear scheme allows for significant gains in computational efficiency despite the fact that the time discretisation is merely of order one (see Figure 4). At the same time, when using high-order schemes, the use of nonlinear strategies seems to be needed to preserve the positivity of the discrete solution. Indeed, the solutions computed with the linear exponential fitting scheme demonstrates spurious oscillations (inducing some negative values) whose intensity increases with the order of approximation.

<figure><img src="image_12.png" /><figcaption>Figure 4 – Efficiency of non-linear HHO methods.</figcaption></figure>
At the end, these promising results pave the way for the development and analysis of similar high-order and efficient schemes for complex dissipative systems.

The Chapter 4 is an extended version of the proceeding [212] for the conference FVCA X.

## About the codes

The HHO codes related to Chapter 4 are implemented in C++. The codes rely on an existing code developed by Simon Lemaire for a linear Poisson equation, with isotropic diffusion and homogeneous Dirichlet boundary conditions on 2 dimensional meshes. This initial code relies on equal-order HHO method, and the consistent part used is the (continuous) gradient of the potential reconstruction. Three codes were developed for the need of the last chapter:

(i) **hfv.fe**: a test code which was a preliminary work in order to understand how to write the finite volume scheme **hfv_evol** in a finite element flavour, with an implementation allowing a natural high-order extension. The nonlinear schemes implemented aim to approximate solution to transient linear advection-diffusion equations with homogeneous Neumann boundary conditions. The main attempt was to write separately the consistent term and the stabilisation term, in such a way that polynomial unknowns could be use instead of constant ones, while keeping an accuracy of order 2. The unknown used was either the density $u$ or the potential $\ell = \log(u)$. The numerical results obtained with this code are not shown in this manuscript, but they indicated the importance of using an appropriate stabilisation term. These results also motivated the use of mixed-order HHO methods instead of equal-order one, in view of the complexity of both analysis and implementation of a relevant stabilisation term.

(ii) **hho_expfitt**: a code for linear anisotropic advection-diffusion equations with mixed Dirichlet-Neumann boundary conditions and general source term. The scheme implemented (4.14) is based on the exponential fitting strategy already used in Chapter 1, and is therefore a linear method. The main work regarding this code was to adapt the existing code to mixed-order method, implement general boundary conditions and discrete gradient operators. An other issue lies in the non-polynomial reconstruction of the discrete density (see Remark 19). There are two versions of this code: one for stationary problems (with does not handle pure Neumann boundary conditions) and one for transient problems. This code was used for the numerical results of Chapter 4.

(iii) **hho.nl**: a code for transient anisotropic advection-diffusion equation with homogeneous Neumann boundary conditions. The scheme implemented is the nonlinear HHO scheme (4.20), main goal of this PhD. The main difficulties lies in the use of non-polynomial reconstructed densities and the nonlinearity of the scheme. Moreover, the use of high-order polynomials induces many numerical instabilities which hinders the convergence the Newton method. Various strategies had to be implemented in order to get a reliable solver, including truncation, filtering of oscillations and loop break for divergent Newton methods. The different tricks used are discussed in Section 4.4.1. The structure of this code is different from the previous ones, and essentially relies on two classes: one for representing the discrete unknowns (both potentials and densities), and an other one to handle the Newton method. This code was used for the numerical results of Chapter 4.
Semiconductors models with irregular convection fields (Appendices A-B-C-D)

The Appendices are dedicated to the numerical approximation of realistic modern semi-
conductor devices. This work began following Patricio Farrell’s visit to Lille in October 2021,
during which he presented to the Rapsodi team various projects and open questions concerning
numerical approximations of real-life models of modern semiconductors. One of them was about
semiconductor models which takes into account what is sometimes called random fluctuations
in the alloys constituting the devices.

In a very simplified way, this consists in a drift-diffusion system with Boltzmann statistics
(as the one studied in Chapter 3) where the convection fields are perturbed by some irregular
quantities, called band-edge energies (BEEs). It reads

$$
\left\{
\begin{array}{l}
\partial_t N - \text{div}(\nabla N - N \nabla(\phi - E_c)) = 0 \\
\partial_t P - \text{div}(\nabla P + P \nabla(\phi + E_v)) = 0 \\
-\lambda^2 \text{div}(\nabla \phi) = C + P - N,
\end{array}
\right.
$$

where $E_c$ and $E_v$ are the conduction and valence energies band, which are a characteristic of the material used in the semiconductor device. For classical models and devices, these BEEs are often considered to be constant functions, and therefore do not really appear in the equation since they only intervene trough their gradients. The main originality of the “random alloys fluctuations model” considered here lies in the fact that we do not consider constant BEEs anymore: these functions are now irregular ones. This irregularity comes from intern fluctuation inside the materials, and could even be described by a function which exhibits some jumps between each pair of atoms. This problem can then be reduced in a first time to the numerical approximation of solution to the advection-equation

$$
\partial_t u - \mathrm{div}(\nabla u + u \nabla(\phi + E)) = 0,
$$

where $\phi$ is a regular potential and $E$ is an irregular band edge energy (which can be, to fix ideas, a piecewise constant function). The approach used by P. Farrell and his collaborators [221, 222] to get numerical simulations of such equations relied on a TPFA-Scharfetter–Gummel scheme, used as if the BEE $E$ was regular, either with the discontinuous $E$ or with a regularised version of $E$, obtained thanks to a Gaussian convolution. With this method, simulations provide numerical solutions in accordance with the physical intuition. On the other hand, this approach was lacking in a more robust analytical framework and some guarantees about the quality and the robustness of the approximation. Indeed, even at the continuous level, if $E$ is a discontinuous function, its is not straightforward to define a natural notion of solution to such an equation (for example in $1D$, with a piecewise constant BEE, the gradient of $E$ is a sum of Dirac masses and the corresponding advection field is extremely irregular). In the case one uses a regularisation on $E$, the question of knowing how to regularise in a good way to keep a similar behaviour to the non-regular model also needs to be answered.

Following the presentation of this problem, I suggested that one could use a strategy inspired
from the exponential fitting procedure used in Chapter 1 in order to overcome the irregularity
issue. Indeed, formally, using the change of variable $\rho = \frac{u}{e^{-E}}$, we get the following equation on $\rho$

$$
e^{-E} \partial_t \rho - \text{div} \left( e^{-E} (\nabla \rho + \rho \nabla \phi) \right) = 0,
$$

which does not take into account the undefined $\nabla E$. Therefore, this formal computation suggests
to use the new "exponentially fitted" formulation to define a notion of solution to this irregular
problem, and then define a solution of the irregular problem as $u = \rho e^E$. At the numerical level,
this strategy implies to solve problems with a heterogeneous and discontinuous tensor $e^{-E}$. I then
spent 10 weeks at the WIAS in Berlin during the summer of 2022 to work with Patricio Farrell
on these issues. This stay led to some discussions and collaborations with physics researchers
who worked on these models: Thomas Koprucki and Timo Streckenbach from the WIAS (Berlin,
Germany), and Stefan Schulz, Michael O'Donovan and Robert Finn form the National Tyndall
institute (Cork, Ireland). The content of the appendices presents the work resulting from this
stay and these collaborations.

From a numerical point of view, this work is quite different from the previous ones, since we
focus on TPFA finite volume methods. The main reason to this choice lies in the necessity to be
able to use the methods developed in the existing codes currently used by the collaborators for
the simulation of LEDs. However, we believe that the nonlinear hybrid methods studied during
this PhD could be used in this framework. In fact, in the light of certain difficulties encountered
to generate TPFA-admissible meshes for these problems in the situations studied in practice, the
use of methods handling general meshes could be a solution in order to gain some computational
efficiency, and therefore a first step towards the simulation of bigger devices. The same remark
should apply to high-order schemes, in regards of the encouraging results of Chapter 4 about
efficiency and computational costs. Thus, one can see this project as a preliminary work in order
to apply state of the art numerical methods to current problems in physics and engineering.

In Appendix A, we are first interested in a the simple model of advection-diffusion with
irregular advection field suggested by the physicists, and we propose a theoretical framework to
analyse the equation, as well as different TPFA schemes to approximate the solution. Then, we
discuss and summarise the works done with the physicists about the design of efficient LEDs.

For completeness, the full versions of the articles corresponding to these works are presented
in the last appendices. In Appendix B, we compare the numerical results obtained thanks to
our strategy with numerical results obtained by commercial packages and real experiments.
In Appendix C, we investigate some issues related to the approximation of models with non-
Boltzmann statistics. These issues are well-known in the mathematical community, but are
rather poorly documented in physical communities. Especially, the main interest of this work is
to show that unsuitable approximations lead to non-correct behaviours, which can even be seen
on I-V curves. Last, in Appendix D we investigate numerically the way of designing efficient
LEDs.

The Appendix B is an article submitted [219], in collaboration with Michael O'Donovan, Patricio Farrell, Timo Streckenbach, Thomas Koprucki and Stefan Schulz. Appendix C is an article published in Optical and Quantum Electronics [126], in collaboration with Patricio Farrell, Michael O'Donovan, Stefan Schulz and Thomas Koprucki. Last, the Appendix D is an article published in Applied Physics Letters [133], in collaboration with Robert Finn, Michael O'Donovan, Patricio Farrell, Timo Streckenbach, Thomas Koprucki and Stefan Schulz.

## About the codes

The work of implementation related to this part differs from the previous ones, since it
covers the modification of an existing code for realistic problems and the development of a
demonstration code for toy problems.

In order to investigate the different possible discretisations of advection-diffusion with
irregular advection, a C++ code was implemented from scratch: **tpfa_qw**. It is developed
to approximate the solutions of the model problem $-\text{div}(u\nabla(h(u) + \phi + E)) = f$ in one space
dimension. It can handle both Boltzmann and Blakemore statistics, generic electrostatic potential
$\phi$ and piecewise constant band edge energy $E$. The main specificity of the code lies in the fine
handling of the meshes and the discontinuities. Indeed, one of the strategies considered (the
one before this work for the simulations) is a scheme where the quantities are not continuous on
a given cell. With respect to the HFV codes this code marks the debut of the Object Oriented
Programming in the implementation, and it introduces the general structure used to solve the
nonlinear problem of Chapter 4. In particular, the code includes a class representing the discrete
unknowns, and a class which handles the Newton method. This last class also allows one to test
if the Jacobian is well computed.

On the other hand, an implementation work was performed on an existing code, **ddfermi** [104],
which is designed to simulate drift-diffusion processes in semiconductor devices. This code
handles 3D simulations, and various TPFA discretisations of the system, including the SEDAN
flux [265] and the Bessemoulin-Chatard flux [26]. My main contribution was to add some imple-
mentations of numerical fluxes for non-Boltzmann statistics with varying band-edge energies.
The simulations presented in Appendices B-C-D were obtained via **ddfermi**.
Chapter 1

# Long-time behaviour of hybrid finite volume schemes for advection-diffusion equations: linear and nonlinear approaches

## Outline of the current chapter

<table><tr><td>1.1</td><td>Introduction</td><td>32</td></tr><tr><td>1.2</td><td>Hybrid finite volume discretisation of a variable diffusion problem</td><td>35</td></tr><tr><td>1.2.1</td><td>Mesh and discrete unknowns</td><td>35</td></tr><tr><td>1.2.2</td><td>Discrete problem</td><td>38</td></tr><tr><td>1.2.3</td><td>Well-posedness</td><td>40</td></tr><tr><td>1.3</td><td>Definition of the schemes and well-posedness</td><td>40</td></tr><tr><td>1.3.1</td><td>Standard HFV scheme</td><td>41</td></tr><tr><td>1.3.2</td><td>Exponential fitting HFV scheme</td><td>44</td></tr><tr><td>1.3.3</td><td>Nonlinear HFV scheme</td><td>47</td></tr><tr><td>1.4</td><td>Long-time behaviour</td><td>54</td></tr><tr><td>1.4.1</td><td>Asymptotic behaviour of the standard HFV scheme</td><td>54</td></tr><tr><td>1.4.2</td><td>Asymptotic behaviour of the exponential fitting scheme</td><td>56</td></tr><tr><td>1.4.3</td><td>Asymptotic behaviour of the nonlinear scheme</td><td>57</td></tr><tr><td>1.5</td><td>Numerical results</td><td>60</td></tr><tr><td>1.5.1</td><td>Implementation</td><td>60</td></tr><tr><td>1.5.2</td><td>Long-time behaviour of discrete solutions</td><td>62</td></tr><tr><td>1.5.3</td><td>Positivity of discrete solutions</td><td>64</td></tr><tr><td>1.5.4</td><td>Accuracy of stationary solutions</td><td>65</td></tr><tr><td>1.6</td><td>Conclusion</td><td>67</td></tr><tr><td>1.A</td><td>Functional inequalities</td><td>67</td></tr><tr><td>1.A.1</td><td>Discrete Poincaré inequalities</td><td>67</td></tr></table>
<table><tbody><tr><td>1.A.2 Logarithmic Sobolev inequalities</td><td>68</td></tr><tr><td>1.A.3 Discrete logarithmic Sobolev inequalities</td><td>70</td></tr><tr><td><b>1.B Nonlinear scheme for mixed Dirichlet-Neumann boundary conditions</b></td><td><b>71</b></td></tr><tr><td>1.B.1 Scheme and well-posedness</td><td>71</td></tr><tr><td>1.B.2 Long-time behaviour</td><td>72</td></tr><tr><td><b>1.C Proofs of technical results</b></td><td><b>73</b></td></tr><tr><td>1.C.1 Discrete boundedness by mass and dissipation</td><td>73</td></tr><tr><td>1.C.2 A local comparison result</td><td>75</td></tr></tbody></table>

This chapter is a work in collaboration with Claire Chainais-Hillairet, Maxime Herda and Simon Lemaire published in Numerische Mathematik [70].

We are interested in the long-time behaviour of approximate solutions to anisotropic and heterogeneous linear advection-diffusion equations in the framework of hybrid finite volume (HFV) methods on general polygonal/polyhedral meshes. We consider two linear methods, as well as a new, nonlinear scheme, for which we prove the existence and the positivity of discrete solutions. We show that the discrete solutions to the three schemes converge exponentially fast in time towards the associated discrete steady-states. To illustrate our theoretical findings, we present some numerical simulations assessing long-time behaviour and positivity. We also compare the accuracy of the schemes on some numerical tests in the stationary case.

## 1.1 Introduction

We are interested in the numerical approximation of linear advection-diffusion equations on bounded domains. These equations constitute the main building block in the modelling of more complex problems stemming from physics (e.g., porous media flows [17], or corrosion models [16]), biology, or electronics (semi-conductor devices modelling [258]). Thus, designing reliable numerical schemes to approximate their solutions is a pre-requisite before discretising more complex models. Our aim here is the preservation of some key physical properties of these equations at the discrete level, on a large variety of meshes.

Let $\Omega$ be an open, bounded, connected polytopal subset of $\mathbb{R}^d$, $d \in \{2, 3\}$, with boundary $\partial\Omega$ divided into two disjoint open subsets $\Gamma^D$ and $\Gamma^N$, in such a way that $\partial\Omega = \overline{\Gamma^D} \cup \overline{\Gamma^N}$. We consider the following problem: Find $u : \mathbb{R}_+ \times \Omega \to \mathbb{R}$ solution to

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

where $n$ is the unit normal vector to $\partial\Omega$ pointing outward $\Omega$, and the data satisfy:

- $\Lambda \in L^\infty(\Omega; \mathbb{R}^{d \times d})$ is a symmetric and uniformly elliptic diffusion tensor: there exist $\lambda_b, \lambda_\# > 0$ such that, for a.e. $x$ in $\Omega$, $\xi \cdot \Lambda(x)\xi \ge \lambda_b|\xi|^2$ and $|\Lambda(x)\xi| \le \lambda_\#|\xi|$ for all $\xi \in \mathbb{R}^d$;
- $\phi \in C^1(\bar{\Omega})$ is a regular potential from which derives the advection field $V^\phi = -\Lambda \nabla \phi$, assumed to satisfy $V^\phi \in H(\text{div}; \Omega)$;

- $f \in L^2(\Omega)$ is a source term;

- $g^D \in H^{\frac{1}{2}}(\Gamma^D)$ is a Dirichlet datum, assumed to be the trace on $\Gamma^D$ of $u^D \in H^1(\Omega)$ satisfying $\|u^D\|_{H^1(\Omega)} \le C\|g^D\|_{H^{\frac{1}{2}}(\Gamma^D)}$ for a given $C > 0$;

- $g^N \in L^2(\Gamma^N)$ is a Neumann datum;

- $u^{in} \in L^2(\Omega)$ is an initial datum.

When $|\Gamma^D| = 0$, we assume that the compatibility condition $\int_{\Omega} f + \int_{\partial\Omega} g^N = 0$ holds true, and we denote by $M$ the initial mass such that $M = \int_{\Omega} u^{in}$, which is known to be preserved along time: $\int_{\Omega} u(t) = M$ for almost every $t > 0$. For further use, we also let in that case $u^M = \frac{M}{|\Omega|} \in \mathbb{R}$, and we refer to this quantity as the mass lifting. Advection-diffusion models of the form (1.1) enjoy certain structural properties. First, when the data $f$, $g^D$, $g^N$, and $u^{in}$ are positive, then the solution $u$ is also positive. Second, the asymptotics $t \to \infty$, the so-called long-time behaviour of the solutions, is well understood (see [33, 62, 63, 254] for related models). Indeed, the solution $u$ to (1.1) converges exponentially fast when $t \to \infty$ towards the steady-state $u^\infty$, solution to the stationary problem

$$ \left\{ \begin{array}{ll} -\text{div}(\Lambda(\nabla u^\infty + u^\infty \nabla \phi)) = f & \text{in } \Omega, \\ u^\infty = g^D & \text{on } \Gamma^D, \\ \Lambda(\nabla u^\infty + u^\infty \nabla \phi) \cdot n = g^N & \text{on } \Gamma^N, \end{array} \right. \qquad (1.2) $$

with additional constraint $\int_{\Omega} u^\infty = M$ when $|\Gamma^D| = 0$. The question of the long-time behaviour has been widely studied in the context of many-particle systems, for which the second law of thermodynamics ensures a relaxation of the transient phenomena towards an equilibrium. From a mathematical point of view, this evolution is strongly related to the dissipation of an entropy functional. Such a vision based on entropy dissipation has given birth to the so-called entropy method. As highlighted by Arnold et al. in [8], the successful use of the entropy method in kinetic theory paves the way to extended applications on various dissipative systems. We refer the reader to the book [174] of Jüngel for a presentation of some of these applications. In [33], Bodineau et al. proposed an entropy functional adapted to drift-diffusion equations with non-homogeneous Dirichlet boundary conditions. A direct adaptation of their method allows to conclude in the present case on the exponential convergence in time of the solution to Problem (1.1) towards the solution to Problem (1.2).

Under appropriate assumptions on the data (a sufficient condition, also valid for more general advection fields, is to assume that $\text{div} V^\phi \ge 0$ a.e. in $\Omega$ and $V^\phi \cdot n \le 0$ a.e. on $\Gamma^N$), the stationary Problem (1.2) is coercive and its well-posedness is straightforward. It turns out that, even if such assumptions on the data are not fulfilled, for an advection field of the form $V^\phi = -\Lambda \nabla \phi$, the problem is still coercive in the new unknown $\rho^\infty = u^\infty e^\phi$, so that one can conclude on well-posedness by solving the problem in the new unknown. Concerning the evolution Problem (1.1), the same arguments show the existence and uniqueness of a global weak solution. For general advection fields (not necessarily deriving from a potential), we refer the reader to the results of Droniou [106] (for mixed Dirichlet-Neumann boundary conditions), and Droniou and Vázquez [114] (for pure Neumann boundary conditions) for detailed statements about well-posedness and regularity of the solutions.

When it comes to numerical approximation, the accuracy of the method is not the only important feature. In some applications (e.g., in subsurface modelling, where the mesh often results from seismic analysis), the mesh must be taken as a datum of the problem, and the
numerical method needs to be adapted so as to handle potentially fairly general meshes. In some other applications (e.g., power plant simulation), the preservation of the positivity of the solutions (or better, of the monotonicity properties of the equation) is an important quality criterion. In yet some other applications (e.g., nuclear waste repository management), finally, the reliability of the simulations in very large time proves to be crucial for sustainability purposes. The positivity and long-time behaviour of discrete solutions have been closely studied in the context of standard two-point flux approximation (TPFA) finite volume schemes, for isotropic diffusion (i.e., $\Lambda = \lambda I_d$ with $\lambda : \Omega \to \mathbb{R}_+$) on orthogonal meshes. In [131], Filbet and Herda studied the long-time behaviour of a TPFA scheme for nonlinear boundary-driven Fokker-Planck equations, adapting to the discrete setting the arguments of [33]. In [69], Chainais-Hillairet and Herda proved on a variety of models that a whole family of TPFA schemes (the so-called B-schemes) preserves the exponential decay towards discrete steady-states. The results of [131] and [69] are valid for general advection fields, and a choice of data $|\Gamma^D| > 0$, $f = 0$, $g^D > 0$, $g^N = 0$, and $u^{in} \ge 0$. We also refer to [195, 73, 160, 46] for related schemes and similar issues. However, these TPFA schemes suffer from an intrinsic limitation: the mesh needs to be $\Lambda$-orthogonal, which, in practice, restricts their use to isotropic diffusion tensors and (standard) orthogonal meshes. In order to overcome this limitation, several linear finite volume methods using auxiliary unknowns have been designed (cf. [109] for a presentation of some of these schemes). As highlighted by Droniou in [109], these methods however suffer from a lack of monotonicity, and so do not preserve the positivity of discrete solutions. As a possible remedy, Cancès and Guichard introduced in [56] (see also the seminal paper [55]), for a class of models encompassing (1.1) for pure Neumann boundary conditions and a choice of data $f = 0$, $g^N = 0$, and $u^{in} \ge 0$ with $M > 0$, a nonlinear vertex approximate gradient (VAG) scheme, designed so as to preserve at the discrete level the positivity of the solutions and the entropy structure of the models, for arbitrary anisotropic diffusions and general meshes. Following the same ideas, Cancès et al. devised and analysed in [52] a (nonlinear) positivity-preserving discrete duality finite volume (DDFV) scheme, whose discrete entropy structure and long-time behaviour were fully studied in [51], based on the adaptation to the discrete setting of nonlinear functional inequalities. The DDFV scheme at hand is however limited to the two-dimensional case, and its adaptation to a three-dimensional framework seems difficult (cf. [109]). Let us also mention the work [241] (and the references therein), in which a general framework for the convergence analysis of positivity-preserving nonlinear cell-centred finite volume methods on general meshes is introduced. On another level, it is known that, given adequate assumptions hold on the data, the solutions to Problem (1.1) are regular in space (at least locally). This suggests that the use of high-order methods shall be an interesting track in order to increase the accuracy at fixed computational cost. Recently introduced by Di Pietro et al. in [100], hybrid high-order (HHO) methods can be seen as an arbitrary-order generalisation of hybrid finite volume (HFV) schemes, that were introduced by Eymard et al. in [123] as yet another way to overcome the limitations of TPFA schemes. HFV methods hinge on cell and face unknowns (whence the vocable hybrid), and as such benefit from a unified 2D/3D formulation. HFV methods have also been bridged to the larger family of hybrid mimetic mixed (HMM) methods in [112]. In view of the above elements, the study of HFV methods appears to be a natural first step in order to design structure-preserving high-order (HHO) schemes for Problem (1.1), that shall both increase the accuracy at fixed computational burden, and preserve the key properties (positivity and long-time behaviour) of the model at hand.

In this article, we study three different HFV schemes for Problem (1.1). The first one is the HFV variant of the HMM family of schemes introduced and analysed in the stationary setting by Beirão da Veiga et al. in [19, 108] (note that an arbitrary-order (HHO) generalisation of this scheme has been proposed in [98]). It is a linear scheme, based on a discretisation of the diffusive
and advective fluxes, that is well-posed under a coercivity condition. The second scheme is also a linear one. Its construction is based on exponential fitting, and takes inspiration from ideas in [43] (it also shares some features with the works [189, 205] and [131], which cover general advection fields). This scheme is unconditionally coercive. These two linear schemes are not expected to preserve positivity, which motivates the introduction of the third method. For pure Neumann boundary conditions, and a choice of data $f = 0$, $g^N = 0$, and $u^{in} \ge 0$ with $M > 0$ (see Appendix 1.B for the case of mixed Dirichlet-Neumann boundary conditions), we introduce a nonlinear HFV scheme, that is devised along the lines of the nonlinear VAG and DDFV schemes of [56] and [52, 51], so as to guarantee the positivity of discrete solutions. Our first result, stated in Theorem 1, is the existence of (positive) solutions to this nonlinear scheme. In a second time, we investigate the long-time behaviour of the three schemes at hand. We establish in Theorems 2, 3, and 4 the exponential decay in time of their discrete solutions towards the associated discrete steady-states. We numerically validate our theoretical findings on a set of test-cases and, for completeness, we also compare the accuracy of the three schemes on stationary problems.

The article is organised as follows. In Section 1.2, we introduce the HFV framework (mesh, discrete unknowns and discrete operators) on a steady variable diffusion problem. In Section 1.3, we introduce the three schemes for the transient advection-diffusion problem, and we discuss their well-posedness. In Section 1.4, we study the long-time behaviour of the three schemes, and prove exponential decay to equilibrium. In Section 1.5, we discuss the implementation of the schemes, and provide a numerical validation of our theoretical results, as well as a comparison of the stationary schemes in terms of accuracy. Appendices 1.A, 1.B, and 1.C finally collect some functional inequalities and the proofs of supplementary and auxiliary results.

## 1.2 Hybrid finite volume discretisation of a variable diffusion problem

The aim of this section is to recall the HFV framework on a steady variable diffusion problem, which corresponds to (1.2) without advection term ($V^\phi = 0$). For a detailed presentation of the method, we refer to [123].

### 1.2.1 Mesh and discrete unknowns

The definitions and notation we adopt for the discretisation are essentially the same as in [123]. A discretisation of the (open, bounded) polytopal set $\Omega \subset \mathbb{R}^d$, $d \in \{2, 3\}$, is defined as a triplet $\mathcal{D} = (\mathcal{M}, \mathcal{E}, \mathcal{P})$, where:

- $\mathcal{M}$ (the mesh) is a partition of $\Omega$, i.e., a finite family of nonempty disjoint (open, connected) polytopal subsets $K$ of $\Omega$ (the mesh cells) such that (i) for all $K \in \mathcal{M}$, $|K| > 0$, and (ii) $\bar{\Omega} = \bigcup_{K \in \mathcal{M}} \bar{K}$.

- $\mathcal{E}$ (the set of faces) is a partition of the mesh skeleton $\bigcup_{K \in \mathcal{M}} \partial K$, i.e., a finite family of nonempty disjoint (open, connected) subsets $\sigma$ of $\bar{\Omega}$ (the mesh faces, or mesh edges if $d=2$) such that (i) for all $\sigma \in \mathcal{E}$, $|\sigma| > 0$ and there exists $\mathcal{H}_\sigma$ affine hyperplane of $\mathbb{R}^d$ such that $\sigma \subset \mathcal{H}_\sigma$, and (ii) $\bigcup_{K \in \mathcal{M}} \partial K = \bigcup_{\sigma \in \mathcal{E}} \bar{\sigma}$. We assume that, for all $K \in \mathcal{M}$, there exists $\mathcal{E}_K \subset \mathcal{E}$ (the set of faces of $K$) such that $\partial K = \bigcup_{\sigma \in \mathcal{E}_K} \bar{\sigma}$. For $\sigma \in \mathcal{E}$, we let $\mathcal{M}_\sigma = \{K \in \mathcal{M} | \sigma \in \mathcal{E}_K\}$ be the set of cells whose $\sigma$ is a face. Then, for all $\sigma \in \mathcal{E}$, either $\mathcal{M}_\sigma = \{K\}$ for a cell $K \in \mathcal{M}$, in which case $\sigma$ is a boundary face ($\sigma \subset \partial \Omega$) and we note $\sigma \in \mathcal{E}_{ext}$, or $\mathcal{M}_\sigma = \{K, L\}$ for two cells $K, L \in \mathcal{M}$, in which case $\sigma$ is an interface and we note $\sigma = K|L \in \mathcal{E}_{int}$.
— $\mathcal{P}$ (the set of cell centres) is a finite family $\{x_K\}_{K \in \mathcal{M}}$ of points of $\Omega$ such that, for all $K \in \mathcal{M}$, (i) $x_K \in K$, and (ii) $K$ is star-shaped with respect to $x_K$. Moreover, we assume that the Euclidean (orthogonal) distance $d_{K,\sigma}$ between $x_K$ and the affine hyperplane $\mathcal{H}_\sigma$ containing $\sigma$ is positive (equivalently, the cell $K$ is strictly star-shaped with respect to $x_K$).

For a given discretisation $\mathcal{D}$, we denote by $h_D > 0$ the size of the discretisation (the meshsize), defined by $h_D = \sup_{K \in \mathcal{M}} h_K$ where, for all $K \in \mathcal{M}$, $h_K = \sup_{x,y \in \bar{K}} |x-y|$ is the diameter of the cell $K$. For all $\sigma \in \mathcal{E}$, we let $\bar{x}_\sigma \in \sigma$ be the barycentre of $\sigma$. Finally, for all $K \in \mathcal{M}$, and all $\sigma \in \mathcal{E}_K$, we let $n_{K,\sigma} \in \mathbb{R}^d$ be the unit normal vector to $\sigma$ pointing outward $K$, and $P_{K,\sigma}$ be the (open) pyramid of base $\sigma$ and apex $x_K$ (notice that, when $d=2$, $P_{K,\sigma}$ is always a triangle). Since $|\sigma|$ and $d_{K,\sigma}$ are positive, we have $|P_{K,\sigma}| = \frac{|\sigma|d_{K,\sigma}}{d} > 0$. We depict on Figure 1.1 an example of discretisation. Notice that the mesh cells are not assumed to be convex, neither $x_K$ is assumed to be the barycentre of $K \in \mathcal{M}$. Notice that hanging nodes are seamlessly handled with our assumptions, so that meshes with non-conforming cells are allowed (see the orange cross in Figure 1.1; the cell $K$ therein is treated as an hexagon). We consider the following measure of regularity for the discretisation (which is

<figure><img src="image_13.png" /><figcaption>Figure 1.1 – Two-dimensional discretisation and corresponding notation.</figcaption></figure>

slightly stronger than the ones advocated in [123, Eq. (4.1)] or in [111, Eq. (7.8)-(7.9)]):

$$ \theta_D = \max \left( \max_{K \in \mathcal{M}, \sigma \in \mathcal{E}_K} \frac{h_K}{d_{K,\sigma}}, \max_{\sigma \in \mathcal{E}, K \in \mathcal{M}_\sigma} \frac{h_K^{d-1}}{|\sigma|} \right). \quad (1.3) $$
Notice that $\theta_D \ge 1$, and that for all $K \in \mathcal{M}$,

$$h_K^d \ge |K| = \sum_{\sigma \in \mathcal{E}_K} |P_{K,\sigma}| = \sum_{\sigma \in \mathcal{E}_K} \frac{|\sigma| d_{K,\sigma}}{d} \ge \sum_{\sigma \in \mathcal{E}_K} \frac{h_K^d}{d \theta_D^2} = \frac{|\mathcal{E}_K|}{d \theta_D^2} h_K^d.$$

Thus, the number of faces of any mesh cell is uniformly bounded:

$$\forall K \in \mathcal{M}, \quad |\mathcal{E}_K| \le d\theta_D^2. \qquad (1.4)$$

Also, it is an easy matter to verify that $\max_{\sigma=K|L \in \mathcal{E}_{int}} \max(\frac{d_{K,\sigma}}{d_{L,\sigma}}, \frac{d_{L,\sigma}}{d_{K,\sigma}}) \le \theta_D^{d-1}$. Given $\mathcal{F}$ a family of discretisations, we say that $\mathcal{F}$ is uniformly regular if there exists $\theta \ge 1$ such that for all $D \in \mathcal{F}$, $\theta_D \le \theta$.

We now introduce the set of (hybrid, cell- and face-based) discrete unknowns:

$$\underline{v}_D = \{\underline{v}_D = (v_K)_{K \in \mathcal{M}}, (v_\sigma)_{\sigma \in \mathcal{E}} : v_K \in \mathbb{R}, \forall K \in \mathcal{M}, v_\sigma \in \mathbb{R}, \forall \sigma \in \mathcal{E}\}.$$

Given a mesh cell $K \in \mathcal{M}$, we let $\underline{v}_K = \mathbb{R} \times \mathbb{R}^{|\mathcal{E}_K|}$ be the restriction of $\underline{v}_D$ to $K$, and $\underline{v}_K = (v_K, (v_\sigma)_{\sigma \in \mathcal{E}_K}) \in \underline{v}_K$ be the restriction of a generic element $\underline{v}_D \in \underline{v}_D$ to $K$. Also, for $\underline{v}_D \in \underline{v}_D$, we let $v_M : \Omega \to \mathbb{R}$ and $v_\mathcal{E} : \bigcup_{K \in \mathcal{M}} \partial K \to \mathbb{R}$ be the piecewise constant functions such that

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

In what follows, for any set $X \subset \bar{\Omega}$, we denote by $(\cdot, \cdot)_X$ the inner product in $L^2(X; \mathbb{R}^l)$, for $l \in \{1; d\}$.

In particular, we have $(w_M, v_M)_\Omega = \sum_{K \in \mathcal{M}} |K| w_K v_K$ and $(w_\mathcal{E}, v_\mathcal{E})_{\partial \Omega} = \sum_{\sigma \in \mathcal{E}_{ext}} |\sigma| w_\sigma v_\sigma$. For further

use, we let $\underline{v}_D$ denote the element of $\underline{v}_D$ with all coordinates equal to 1. Also, given a function $f : \mathbb{R} \to \mathbb{R}$, and with a slight abuse in notation, we denote by $f(v_D)$ the element of $\underline{v}_D$ whose coordinates are the $(f(v_K))_{K \in \mathcal{M}}$ and the $(f(v_\sigma))_{\sigma \in \mathcal{E}}$. Finally, we let the product $\underline{w}_D \times \underline{v}_D$ denote the element of $\underline{v}_D$ whose $i$-th coordinate is the product of the $i$-th coordinates of $\underline{w}_D$ and $\underline{v}_D$.

When considering mixed Dirichlet-Neumann boundary conditions, we assume that the discretisation $\mathcal{D}$ is compliant with the partition $\partial \Omega = \bar{\Gamma}^D \cup \bar{\Gamma}^N$ of the boundary of the domain, in the sense that the set $\mathcal{E}_{ext}$ can be split into two (necessarily disjoint) subsets $\mathcal{E}_{ext}^D = \{\sigma \in \mathcal{E}_{ext} | \sigma \subset \Gamma^D\}$ and $\mathcal{E}_{ext}^N = \{\sigma \in \mathcal{E}_{ext} | \sigma \subset \Gamma^N\}$ such that $\mathcal{E}_{ext} = \mathcal{E}_{ext}^D \cup \mathcal{E}_{ext}^N$. Notice that as soon as $|\Gamma^D| > 0$, $|\mathcal{E}_{ext}^D| \ge 1$. We define the following subspace of $\underline{v}_D$, enforcing strongly a homogeneous Dirichlet boundary condition on $\Gamma^D$:

$$\underline{v}_{D,0}^D = \{\underline{v}_D \in \underline{v}_D : v_\sigma = 0 \ \forall \sigma \in \mathcal{E}_{ext}^D\}.$$

In view of the upcoming analysis, we define a discrete counterpart of the $H^1$ seminorm. Locally to any cell $K \in \mathcal{M}$, we let, for any $\underline{v}_K \in \underline{v}_K$, $|\underline{v}_K|^2_{1,K} = \sum_{\sigma \in \mathcal{E}_K} \frac{|\sigma|}{d_{K,\sigma}} (v_K - v_\sigma)^2$. At the global level, for any $\underline{v}_D \in \underline{v}_D$, we let

$$|\underline{v}_D|_{1,D} = \sqrt{\sum_{K \in \mathcal{M}} |\underline{v}_K|^2_{1,K}}.$$

Notice that $|\cdot|_{1,D}$ does not define a norm on $\underline{v}_D$, but if $|\underline{v}_D|_{1,D} = 0$, then there is $c \in \mathbb{R}$ such that $\underline{v}_D = c \underline{v}_D$ ($\underline{v}_D$ is constant). Thus, $|\cdot|_{1,D}$ defines a norm on the space $\underline{v}_{D,0}^D$ as soon as $|\Gamma^D| > 0$, as
well as on the space of zero-mass vectors

$$ \underline{V}_{D,0}^{N} = \{ \underline{v}_{D} \in \underline{V}_{D} : \int_{\Omega} v_{M} = 0 \}. $$

For further use, and to allow for a seamless treatment of pure Neumann boundary conditions, we introduce the notation $\underline{v}_{D,0}$, to denote either $\underline{v}_{D,0}^N$ whenever $|\Gamma^D| = 0$, or $\underline{v}_{D,0}^D$ otherwise.

### 1.2.2 Discrete problem

The HFV method hinges on the definition of a discrete gradient operator $\nabla_D$, that maps any element $\underline{v}_D \in \underline{V}_D$ to a piecewise constant $\mathbb{R}^d$-valued function on the pyramidal submesh of $\mathcal{M}$ formed by all the $P_{K,\sigma}$'s, for $K \in \mathcal{M}$ and $\sigma \in \mathcal{E}_K$. More precisely, for all $K \in \mathcal{M}$, and all $\sigma \in \mathcal{E}_K$,

$$ \nabla_D \underline{v}_{D|K} = \nabla_K \underline{v}_K \quad \text{with} \quad \nabla_K \underline{v}_{K|_{P_{K,\sigma}}} = \nabla_{K,\sigma} \underline{v}_K = G_K \underline{v}_K + S_{K,\sigma} \underline{v}_K \in \mathbb{R}^d, $$

where $G_K \underline{v}_K$ is the consistent part of the gradient given by

$$ G_K \underline{v}_K = \frac{1}{|K|} \sum_{\sigma' \in \mathcal{E}_K} |\sigma'| (v_{\sigma'} - v_K) n_{K,\sigma'} = \frac{1}{|K|} \sum_{\sigma' \in \mathcal{E}_K} |\sigma'| v_{\sigma'} n_{K,\sigma'}, $$

and $S_{K,\sigma} \underline{v}_K$ is a stabilisation given, for some free parameter $\eta > 0$, by

$$ S_{K,\sigma} \underline{v}_K = \frac{\eta}{d_{K,\sigma}} (v_{\sigma} - v_K - G_K \underline{v}_K \cdot (\bar{x}_{\sigma} - x_K)) n_{K,\sigma}. \quad (1.5) $$

**Remark 1** (Choice of $\eta$). There are two specific values of the stabilisation parameter $\eta > 0$ for which one recovers known numerical schemes from the literature:

(i) for $\eta = \sqrt{d}$, one recovers the original HFV scheme of [123], that coincides with the TPFA scheme on super-admissible meshes (see [123, Lemma 2.10]);

(ii) for $\eta = d$, one recovers the Discrete Geometric Approach (DGA) of [88], later bridged to the non-conforming finite element setting in [101].

The influence of the value of $\eta$ on the numerical results has been investigated in [35] for anisotropic diffusion problems. It is shown that the above two values are appropriate choices (neither under- nor over-penalised).

Let us consider the stationary problem (1.2), without advection term ($V^\phi = 0$). Our aim is to write an HFV discretisation of this steady variable diffusion problem. Locally to any cell $K \in \mathcal{M}$, we introduce the discrete bilinear form $a_K^\Lambda : \underline{V}_K \times \underline{V}_K \to \mathbb{R}$ such that, for all $\underline{u}_K, \underline{v}_K \in \underline{V}_K$,

$$ a_K^\Lambda(\underline{u}_K, \underline{v}_K) = \sum_{\sigma \in \mathcal{E}_K} |P_{K,\sigma}| \nabla_{K,\sigma} \underline{v}_K \cdot \Lambda_{K,\sigma} \nabla_{K,\sigma} \underline{u}_K = (\Lambda \nabla_K \underline{u}_K, \nabla_K \underline{v}_K)_K, \quad (1.6) $$

where we set $\Lambda_{K,\sigma} = \frac{1}{|P_{K,\sigma}|} \int_{P_{K,\sigma}} \Lambda$. At the global level, we let $a_D^\Lambda : \underline{V}_D \times \underline{V}_D \to \mathbb{R}$ be the discrete bilinear form such that, for all $\underline{u}_D, \underline{v}_D \in \underline{V}_D$,

$$ a_D^\Lambda(\underline{u}_D, \underline{v}_D) = \sum_{K \in \mathcal{M}} a_K^\Lambda(\underline{u}_K, \underline{v}_K) = (\Lambda \nabla_D \underline{u}_D, \nabla_D \underline{v}_D)_{\Omega}. $$
The discrete HFV problem then reads: Find $\underline{u}_D^z \in \underline{V}_{D,0}$ such that

$$a_D^\Lambda(\underline{u}_D^z, \underline{v}_D) = (f, v_M)_\Omega + (g^N, v_\epsilon)_{\Gamma^N} - a_D^\Lambda(\underline{u}_D^l, \underline{v}_D) \quad \forall \underline{v}_D \in \underline{V}_{D,0}, \quad (1.7)$$

where $\underline{u}_D^l \in \underline{V}_D$ is equal

(i) either, when $|\Gamma^D| > 0$, to the HFV interpolate $\underline{u}_D^D$ of the known lifting $u^D$ of the Dirichlet datum $g^D$ (satisfying $|\underline{u}_D^D|_{1,D} \le C_{l,\Gamma^D} \|g^D\|_{H^{1/2}(\Gamma^D)}$, where $C_{l,\Gamma^D} > 0$ only depends on the discretisation $\mathcal{D}$ through $\theta_D$),

(ii) or, when $|\Gamma^D| = 0$, to $\underline{u}_D^M = u^M \underline{1}_D$, where we recall that $u^M = \frac{M}{|\Omega|}$ is the mass lifting (remark that $a_D^\Lambda(\underline{u}_D^M, \underline{v}_D) = 0$ for all $\underline{v}_D \in \underline{V}_D$),

and the approximation of the solution to (1.2), denoted $\underline{u}_D^\infty \in \underline{V}_D$, is finally defined as

$$\underline{u}_D^\infty = \underline{u}_D^z + \underline{u}_D^l. \quad (1.8)$$

Let us note that the superscript z stands for “zero”, while l stands for “lifting”.

Problem (1.7) defines a finite volume method, in the sense that it can be equivalently rewritten under a conservative form, with local mass balance, flux equilibration at interfaces, and boundary conditions. For all $K \in \mathcal{M}$, and all $\sigma \in \mathcal{E}_K$, the normal diffusive flux $-\int_\sigma \Lambda \nabla u^\infty \cdot n_{K,\sigma}$ is approximated by the following numerical flux:

$$F_{K,\sigma}^{\Lambda}(\underline{u}_K) = \sum_{\sigma' \in \mathcal{E}_K} A_K^{\sigma\sigma'}(u_K - u_{\sigma'}), \quad (1.9)$$

where the $A_K^{\sigma\sigma'}$ are defined by

$$A_K^{\sigma\sigma'} = \sum_{\sigma'' \in \mathcal{E}_K} |P_{K,\sigma''}| y_K^{\sigma''\sigma} \cdot \Lambda_{K,\sigma''} y_K^{\sigma''\sigma'}, \quad (1.10)$$

and the $y_K^{\sigma\sigma'} \in \mathbb{R}^d$ only depend on the geometry of the discretisation $\mathcal{D}$ (see, for example, [123, Eq. (2.22)]) for an exact definition with $\eta = \sqrt{d}$. For all $K \in \mathcal{M}$, one can express the local discrete bilinear form $a_K^\Lambda$ in terms of the local fluxes $(F_{K,\sigma}^\Lambda)_{\sigma \in \mathcal{E}_K}$: for all $\underline{u}_K, \underline{v}_K \in \underline{V}_K$,

$$a_K^\Lambda(\underline{u}_K, \underline{v}_K) = \sum_{\sigma \in \mathcal{E}_K} F_{K,\sigma}^\Lambda(\underline{u}_K)(\underline{v}_K - \underline{v}_\sigma). \quad (1.11)$$

As for the VAG [56] and DDFV [52, 51] schemes, we can also express the local discrete bilinear form in a different way, which will be useful in the sequel:

$$a_K^\Lambda(\underline{u}_K, \underline{v}_K) = \delta_K \underline{v}_K \cdot A_K \delta_K \underline{u}_K, \quad (1.12)$$

where, for all $\underline{v}_K \in \underline{V}_K$, $\delta_K \underline{v}_K \in \mathbb{R}^{|\mathcal{E}_K|}$ is defined by

$$\delta_K \underline{v}_K = (\underline{v}_K - \underline{v}_\sigma)_{\sigma \in \mathcal{E}_K},$$

and $A_K \in \mathbb{R}^{|\mathcal{E}_K| \times |\mathcal{E}_K|}$ is the symmetric (because $\Lambda$ is) positive semi-definite matrix whose entries are the $A_K^{\sigma\sigma'}$, that can actually be proved to be nonsingular (cf. Lemma 4).
### 1.2.3 Well-posedness

As for the continuous case, the well-posedness of HFV methods for diffusion problems relies on a coercivity argument. Let $K \in \mathcal{M}$, and reason locally. By definition (1.6) of the local discrete bilinear form $a_K^\Lambda$, and from the bounds on the diffusion coefficient, we have $\lambda_b \|\nabla_K \underline{v}_K\|_{L^2(K;\mathbb{R}^d)}^2 \le a_K^\Lambda(\underline{v}_K, \underline{v}_K) \le \lambda_\# \|\nabla_K \underline{v}_K\|_{L^2(K;\mathbb{R}^d)}^2$ for all $\underline{v}_K \in \underline{V}_K$. Furthermore, the following comparison result holds (cf. [111, Lemma 13.11, p = 2] and its proof): there exist $\alpha_b, \alpha_\#$ with $0 < \alpha_b \le \alpha_\# < \infty$, only depending on $\Omega, d$, and $\theta_D$ such that $\alpha_b|\underline{v}_K|_{1,K}^2 \le \|\nabla_K \underline{v}_K\|_{L^2(K;\mathbb{R}^d)}^2 \le \alpha_\#|\underline{v}_K|_{1,K}^2$ for all $\underline{v}_K \in \underline{V}_K$. Combining both estimates, we infer a local coercivity and boundedness result:

$$ \forall \underline{v}_K \in \underline{V}_K, \quad \lambda_b \alpha_b |\underline{v}_K|_{1,K}^2 \le a_K^\Lambda(\underline{v}_K, \underline{v}_K) \le \lambda_\# \alpha_\# |\underline{v}_K|_{1,K}^2. \quad (1.13) $$

Summing over $K \in \mathcal{M}$, we get the following global estimates:

$$ \forall \underline{v}_D \in \underline{V}_D, \quad \lambda_b \alpha_b |\underline{v}_D|_{1,D}^2 \le a_D^\Lambda(\underline{v}_D, \underline{v}_D) \le \lambda_\# \alpha_\# |\underline{v}_D|_{1,D}^2. \quad (1.14) $$

The well-posedness of Problem (1.7)-(1.8) follows.

**Proposition 1** (Well-posedness). There exists a unique solution $\underline{u}_D^\infty \in \underline{V}_D$ to Problem (1.7)-(1.8), which satisfies $|\underline{u}_D^\infty|_{1,D} \le C(\|f\|_{L^2(\Omega)} + \|g^N\|_{L^2(\Gamma^N)} + \|g^D\|_{H^{1/2}(\Gamma^D)})$, for some $C > 0$ depending on the data, and on the discretisation $\mathcal{D}$ only through $\theta_D$.

*Proof.* The existence/uniqueness of $\underline{u}_D^z \in \underline{V}_{D,0}$ solution to (1.7) (and in turn of $\underline{u}_D^\infty = \underline{u}_D^z + \underline{u}_D^l$) is a direct consequence of the coercivity estimate (1.14), and of the fact that $|\cdot|_{1,D}$ defines a norm on $\underline{V}_{D,0}$ (recall that $\underline{V}_{D,0}$ denotes either $\underline{V}_{D,0}^N$ when $|\Gamma^D| = 0$, or $\underline{V}_{D,0}^D$ otherwise). To prove the bound on $|\underline{u}_D^\infty|_{1,D}$, we use the triangle inequality:

$$ |\underline{u}_D^\infty|_{1,D} \le |\underline{u}_D^z|_{1,D} + |\underline{u}_D^l|_{1,D}. $$

To estimate the first term, we test Problem (1.7) with $\underline{v}_D = \underline{u}_D^z \in \underline{V}_{D,0}$, we use (1.14), and we apply the Cauchy-Schwarz inequality. We get

$$ \lambda_b \alpha_b |\underline{u}_D^z|_{1,D}^2 \le \|f\|_{L^2(\Omega)} \|u_M^z\|_{L^2(\Omega)} + \|g^N\|_{L^2(\Gamma^N)} \|u_\epsilon^z\|_{L^2(\Gamma^N)} + \lambda_\# \alpha_\# |\underline{u}_D^l|_{1,D} |\underline{u}_D^z|_{1,D}. $$

Using a discrete Poincaré inequality, recalled in Proposition 7, and applied to $\underline{u}_D^z \in \underline{V}_{D,0}$, as well as the discrete trace inequality of [111, Eq. (B.58), p = 2] combined with a discrete Poincaré inequality, we obtain

$$ |\underline{u}_D^z|_{1,D} \le C_1 (\|f\|_{L^2(\Omega)} + \|g^N\|_{L^2(\Gamma^N)}) + C_2 |\underline{u}_D^l|_{1,D}. $$

It remains to estimate the norm of the lifting $|\underline{u}_D^l|_{1,D}$:

(i) if $|\Gamma^D| > 0$, one has $\underline{u}_D^l = \underline{u}_D^D$ (interpolate of the lifting $u^D$), therefore $|\underline{u}_D^l|_{1,D} \le C_{l,\Gamma^D} \|g^D\|_{H^{1/2}(\Gamma^D)}$;

(ii) if $|\Gamma^D| = 0$, since $\underline{u}_D^l = \underline{u}_D^M = u^M \underline{1}_D$, the lifting is constant and $|\underline{u}_D^l|_{1,D} = 0$. $\square$

## 1.3 Definition of the schemes and well-posedness

In this section, we introduce and study the well-posedness of three HFV schemes, two linear ones and a nonlinear scheme, for the time-dependent advection-diffusion problem (1.1). For
the first two (linear) schemes, we introduce and study in the first place their steady versions on Problem (1.2). For the nonlinear scheme, by anticipation of the asymptotic analysis of Section 1.4.3, we restrict our study to the case where the (positive) solution to Problem (1.1) converges in long time towards the so-called thermal equilibrium (see (1.25)). However, as it will be verified numerically in Section 1.5.4 in the stationary setting, our scheme is applicable to more general data. We consider a fixed spatial discretisation $\mathcal{D}$ of $\Omega$, which satisfies the conditions detailed in Section 1.2.1, and a fixed time step $\Delta t > 0$ for the time discretisation.

**Remark 2 (Linear schemes and nonhomogeneous data).** The linearity of Problem (1.1) implies that, (i) if $|\Gamma^D| > 0$, the shifted variable $u^z = u - u^D$ (recall that $u^D$ is a known lifting of the Dirichlet datum $g^D$) satisfies an advection-diffusion equation with zero Dirichlet boundary condition on $\Gamma^D$, and (ii) otherwise, the shifted variable $u^z = u - u^M$ (recall that $u^M = \frac{M}{|\Omega|}$ is the mass lifting) satisfies a (compatible) pure Neumann advection-diffusion equation with zero-mass constraint. Thus, without loss of generality, we can restrict our study to the homogeneous case $g^D = 0$ or $M = 0$. For an example (in the steady, purely diffusive case) of how to handle at the discrete level nonhomogeneous data $g^D$ or $M$, we refer the reader to Problem (1.7)-(1.8) and Proposition 1 above. Notice that such manipulations are possible for linear schemes only.

### 1.3.1 Standard HFV scheme

We present here the HFV variant of the HMM family of schemes introduced in [19, 108].

#### Stationary problem

We consider Problem (1.2) with $g^D = 0$ when $|\Gamma^D| > 0$, or $M = \int_{\Omega} u^\infty = 0$ otherwise (cf. Remark 2). Locally to any cell $K \in \mathcal{M}$, we introduce the discrete bilinear form $a_K: \underline{V}_K \times \underline{V}_K \to \mathbb{R}$ such that, for all $\underline{u}_K, \underline{v}_K \in \underline{V}_K$,

$$a_K(\underline{u}_K, \underline{v}_K) = a_K^\Lambda(\underline{u}_K, \underline{v}_K) + a_K^\phi(\underline{u}_K, \underline{v}_K), \quad (1.15)$$

where the diffusive part $a_K^\Lambda$ is defined by (1.6) (and rewrites as (1.11) in terms of the local diffusive fluxes $F_{K,\sigma}^\Lambda(\underline{u}_K)$ given by (1.9)), and the advective part $a_K^\phi$ is defined by

$$a_K^\phi(\underline{u}_K, \underline{v}_K) = \sum_{\sigma \in \mathcal{E}_K} F_{K,\sigma}^\phi(\underline{u}_K)(\underline{v}_K - \underline{v}_\sigma), \quad (1.16)$$

with $F_{K,\sigma}^\phi(\underline{u}_K)$ an approximation of the normal advective flux $\int_\sigma u^\infty V^\phi \cdot n_{K,\sigma}$. In order to define the numerical advective fluxes, we need to introduce some data. We set $V_{K,\sigma}^\phi = \frac{1}{|\sigma|} \int_\sigma V^\phi \cdot n_{K,\sigma}$, and $\mu_\sigma = \min(1, \min_{K \in \mathcal{M}_\sigma} \text{Sp}(\Lambda_K)) > 0$ where $\Lambda_K = \frac{1}{|K|} \int_K \Lambda$. We could as well use the finer local Péclet number introduced in [80], namely consider the value $\mu_\sigma = \min(1, \min_{K \in \mathcal{M}_\sigma} (n_{K,\sigma} \cdot \Lambda_K n_{K,\sigma}))$, but we choose here to stick to the formula advocated in [19, 108]. We also consider a Lipschitz continuous function $A: \mathbb{R} \to \mathbb{R}$, satisfying the following conditions:

$$\begin{aligned} A(0) &= 0, \\ \forall s \in \mathbb{R}, \quad &A(-s) - A(s) = s, \\ \forall s \in \mathbb{R}, \quad &A(-s) + A(s) \ge 0. \end{aligned} \quad (1.17)$$
Notice that $A = B - 1$, where $B$ is the classical function used for the $B$-schemes introduced in [66]. In the $B$-schemes framework, advection and diffusion are simultaneously treated in the definition of the numerical flux. Here, as in [19], only the advective part is considered, whence the fact that $A = B - 1$. Standard choices of $A$ functions include:

- the centred discretisation: $A : s \mapsto -\frac{s}{2}$;

- the upwind discretisation: $A : s \mapsto \max(-s, 0)$;

- the Scharfetter-Gummel discretisation: $A : s \mapsto \begin{cases} \frac{s}{e^s - 1} - 1 & \text{if } s \neq 0 \\ 0 & \text{if } s = 0 \end{cases}$.

We eventually define, for all $K \in \mathcal{M}$, and all $\sigma \in \mathcal{E}_K$, the numerical advective flux: for all $\underline{u}_K \in \underline{V}_K$,

$$F_{K,\sigma}^{\phi}(\underline{u}_K) = |\sigma| \frac{\mu_{\sigma}}{d_{K,\sigma}} \left( A \left( -\frac{d_{K,\sigma}}{\mu_{\sigma}} V_{K,\sigma}^{\phi} \right) u_K - A \left( \frac{d_{K,\sigma}}{\mu_{\sigma}} V_{K,\sigma}^{\phi} \right) u_{\sigma} \right). \quad (1.18)$$

Letting $a_D: \underline{V}_D \times \underline{V}_D \to \mathbb{R}$ be the (global) discrete bilinear form such that, for all $\underline{u}_D, \underline{v}_D \in \underline{V}_D$,

$$a_D(\underline{u}_D, \underline{v}_D) = \sum_{K \in \mathcal{M}} a_K(\underline{u}_K, \underline{v}_K), \quad (1.19)$$

and recalling that $\underline{V}_{D,0}$ denotes either $\underline{V}_{D,0}^N$ whenever $|\Gamma^D| = 0$, or $\underline{V}_{D,0}^D$ otherwise, the discrete problem reads: Find $\underline{u}_D^\infty \in \underline{V}_{D,0}$ such that

$$a_D(\underline{u}_D^\infty, \underline{v}_D) = (f, v_M)_{\Omega} + (g^N, v_\varepsilon)_{\Gamma^N} \quad \forall \underline{v}_D \in \underline{V}_{D,0}. \quad (1.20)$$

Remark that, for pure Neumann boundary conditions with $M \neq 0$, as opposed to the purely diffusive case of Problem (1.7), $a_D(\underline{u}_D^M, \underline{v}_D) \neq 0$ a priori for $\underline{v}_D \in \underline{V}_{D,0}^N$. The well-posedness of (1.20) is discussed in the following proposition.

**Proposition 2 (Well-posedness).** Let A be a Lipschitz continuous function satisfying (1.17). If the advection field $V^\phi$ satisfies the two following conditions:

$$(i) \text{ almost everywhere on } \Gamma^N, V^\phi \cdot n \le 0, \quad (1.21a)$$

$$(ii) \text{ there exists } \beta < \frac{2\lambda_b\alpha_b}{C_P^2} \text{ such that, almost everywhere in } \Omega, \text{ div } V^\phi \ge -\beta, \quad (1.21b)$$

where $\lambda_b\alpha_b$ is the coercivity constant of (1.14), and $C_P$ is either equal to $C_{PW}$ if $|\Gamma^D| = 0$ or to $C_{P,\Gamma^D}$ otherwise (where $C_{PW}, C_{P,\Gamma^D}$ are the Poincaré constants of Proposition 7), then there exists $\kappa > 0$, only depending on $\Lambda, \beta, \Omega, d, \Gamma^D$, and $\theta_D$ such that

$$\forall \underline{v}_D \in \underline{V}_{D,0}, \quad a_D(\underline{v}_D, \underline{v}_D) \ge \kappa |\underline{v}_D|^2_{1,D}. \quad (1.22)$$

Consequently, there exists a unique solution $\underline{u}_D^\infty \in \underline{V}_{D,0}$ to Problem (1.20). Moreover, one has

$$|\underline{u}_D^\infty|_{1,D} \le C (\|f\|_{L^2(\Omega)} + \|g^N\|_{L^2(\Gamma^N)}),$$

for some $C > 0$ depending on the data, and on the discretisation $D$ only through $\theta_D$.

*Proof.* Let $K \in \mathcal{M}$, and $\sigma \in \mathcal{E}_K$. Let $s_{K,\sigma} = \frac{d_{K,\sigma}}{\mu_{\sigma}} V_{K,\sigma}^{\phi}$ and $\zeta_{K,\sigma} = A(-s_{K,\sigma}) + A(s_{K,\sigma})$. According to (1.17), $\zeta_{K,\sigma} \ge 0$, and we have $A(-s_{K,\sigma}) = \frac{s_{K,\sigma} + \zeta_{K,\sigma}}{2}$ and $-A(s_{K,\sigma}) = \frac{s_{K,\sigma} - \zeta_{K,\sigma}}{2}$. Consequently, for
all $\underline{v}_K \in \underline{V}_K$,

$$ (A(-s_{K,\sigma})v_K - A(s_{K,\sigma})v_{\sigma})(v_K - v_{\sigma}) = \frac{1}{2}s_{K,\sigma}(v_K^2 - v_{\sigma}^2) + \frac{1}{2}\zeta_{K,\sigma}(v_K - v_{\sigma})^2 \geq \frac{1}{2}s_{K,\sigma}(v_K^2 - v_{\sigma}^2). $$

Recalling (1.16) and (1.18), we infer that, for all $\underline{v}_D \in \underline{V}_D$,

$$ a_D^{\phi}(\underline{v}_D, \underline{v}_D) = \sum_{K \in \mathcal{M}} a_K^{\phi}(\underline{v}_K, \underline{v}_K) \geq \frac{1}{2} \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{E}_K} |\sigma| V_{K, \sigma}^{\phi} (v_K^2 - v_{\sigma}^2). $$

Since, for all $K \in \mathcal{M}$, $\sum_{\sigma \in \mathcal{E}_K} |\sigma| V_{K, \sigma}^{\phi} = \int_K \mathrm{div} V^{\phi}$ and, for all $\sigma \in \mathcal{E}_{int}$, $\sum_{K \in \mathcal{M}_{\sigma}} |\sigma| V_{K, \sigma}^{\phi} = 0$, we have, for all $\underline{v}_D \in \underline{V}_{D,0}$,

$$ a_D^{\phi}(\underline{v}_D, \underline{v}_D) \geq \frac{1}{2}(\mathrm{div} V^{\phi}, v_M^2)_{\Omega} - \frac{1}{2}(V^{\phi} \cdot n, v_{\varepsilon}^2)_{\Gamma^N}. \quad (1.23) $$

Combining (1.15) with (1.19), (1.21) with (1.23), and the coercivity result (1.14), we deduce that, for all $\underline{v}_D \in \underline{V}_{D,0}$,

$$ a_D(\underline{v}_D, \underline{v}_D) \geq \lambda_b \alpha_b |\underline{v}_D|^2_{1,D} - \frac{\beta}{2} \|v_M\|_{L^2(\Omega)}^2. \quad (1.24) $$

Using a Poincaré inequality from Proposition 7, one has, for all $\underline{v}_D \in \underline{V}_{D,0}$,

$$ a_D(\underline{v}_D, \underline{v}_D) \geq \left( \lambda_b \alpha_b - \frac{\beta C_p^2}{2} \right) |\underline{v}_D|^2_{1,D}, $$

therefore the estimate (1.22) holds for $\kappa = \lambda_b \alpha_b - \frac{\beta C_p^2}{2}$, which is positive according to (1.21b). The existence/uniqueness of $\underline{u}_D^\infty \in \underline{V}_{D,0}$ solution to (1.20) is a direct consequence of the coercivity estimate (1.22), and of the fact that $|\cdot|_{1,D}$ defines a norm on $\underline{V}_{D,0}$. The continuous dependency of $\underline{u}_D^\infty$ with respect to the data can then be proved as in the proof of Proposition 1. $\square$

**Remark 3** (Assumptions on the advection field). The well-posedness result of Proposition 2 does not use the fact that $V^{\phi}$ is related to the gradient of a potential, and thus extends to general advection fields. Even better, under a smallness assumption on the meshsize $h_D$, it is actually possible to prove well-posedness for Problem (1.20) without assumptions (1.21) on the advection field. The starting point to prove so is a discrete Gårding inequality like (1.24), which can be easily obtained in full generality from (1.23) (in the case $|\Gamma^N| > 0$, it is obtained from the multiplicative discrete trace inequality of [111, Eq. (B.57), p = 2] and holds for $h_D$ sufficiently small). The proof then proceeds by contradiction, assuming that a discrete inf-sup condition does not hold in the limit $h_D \to 0$, and using a compactness argument (cf. [111, Lemmas B.27-B.33, p = 2]), together with the unconditional well-posedness of the continuous problem (1.2) (cf. [106]).

**Remark 4** (Choice of A). The choice of the function A is of great importance. In particular, the Scharfetter-Gummel approximation is rather classical in various contexts. First introduced in [239] in the framework of TPFA schemes, this approximation of the flux ensures the preservation of the so-called thermal (or Gibbs) equilibrium at the discrete level, which has the form:

$$ u_{th}^{\infty} = \bar{\rho} e^{-\phi} \quad (1.25) $$

where $\bar{\rho} \in \mathbb{R}_+$ is prescribed by the data. For instance, for pure Neumann boundary conditions, $f = 0$, and $g^N = 0$, we have $\bar{\rho} = M / \int_{\Omega} e^{-\phi}$. The discrete solution obtained with the TPFA scheme is then the interpolate of $u_{th}^\infty$. This property is no more true for the hybrid scheme and we observe numerically
that, for $\bar{\rho} \neq 0$, the discrete solution $u_D^\infty \in \underline{V}_D$ is in general not the HFV interpolate of $u_{th}^\infty$. However, as explained in [108, pp. 553-554], provided the parameters $(\mu_\sigma)_{\sigma \in \mathcal{E}}$ are well-chosen, the Scharfetter-Gummel flux ensures an automatic upwinding to the scheme in the advection-dominated regime (whereas it degenerates towards the centred scheme in the diffusion-dominated regime).

## Evolution problem

We consider Problem (1.1) with $g^D = 0$ when $|\Gamma^D| > 0$, or $M = \int_{\Omega} u^{in} = 0$ otherwise (see Remark 2). We use a backward Euler discretisation in time, and the HFV discretisation introduced in Section 1.3.1 in space. The discrete problem reads: Find $(u_D^n \in \underline{V}_{D,0})_{n \ge 1}$ such that

$$ \left\{ \begin{array}{ll} \frac{1}{\Delta t}(u_M^n - u_M^{n-1}, v_M)_\Omega + a_D(\underline{u}_D^n, \underline{v}_D) = (f, v_M)_\Omega + (g^N, v_\mathcal{E})_{\Gamma^N} & \forall \underline{v}_D \in \underline{V}_{D,0}, \\ u_K^0 = \frac{1}{|K|} \int_K u^{in} & \forall K \in \mathcal{M}, \end{array} \right. \quad (1.26a) $$

$$ u_K^0 = \frac{1}{|K|} \int_K u^{in} \quad \forall K \in \mathcal{M}, \quad (1.26b) $$

where $a_D$ is defined by (1.15) and (1.19). Since $a_D$ is coercive, the bilinear form in (1.26a) is also coercive, so the scheme (1.26) is well-posed under the assumptions (1.21) on the advection field, as a straightforward consequence of Proposition 2.

**Remark 5 (Pure Neumann case).** When considering pure Neumann boundary conditions, and contrary to the stationary case, one can actually seek at each time step for a solution to Problem (1.26) in $\underline{V}_D$, i.e., it is not necessary to seek for a solution in the constrained space $\underline{V}_{D,0}^N$. Indeed, testing (1.26a) by $\underline{1}_D \in \underline{V}_D$, and using that $\int_{\Omega} f + \int_{\partial\Omega} g^N = 0$ and $a_D(\underline{u}_D^n, \underline{1}_D) = 0$, one can automatically infer that $\int_{\Omega} u_M^n = \int_{\Omega} u_M^{n-1}$ for all $n \ge 1$, that is, $\int_{\Omega} u_M^n = \int_{\Omega} u_M^0 = \int_{\Omega} u^{in} = M = 0$ for all $n \ge 1$, i.e., $\underline{u}_D^n \in \underline{V}_{D,0}^N$ for all $n \in \mathbb{N}^*$.

## 1.3.2 Exponential fitting HFV scheme

Following ideas in [43] (cf. also [189, 205] and [131] for general advection fields) in the context of finite element methods, we aim to design an unconditionally (i.e., without the need for assumptions (1.21) on the advection field $V^\phi$) coercive scheme for the advection-diffusion problem in the HFV framework.

### Stationary problem

We consider Problem (1.2). The strategy advocated in [43] is based on the following observation: at the continuous level, if $u^\infty$ is a solution to (1.2), letting

$$ \omega = e^{-\phi}, \quad (1.27) $$

we can introduce the Slotboom change of variable $\rho^\infty = \frac{u^\infty}{\omega}$ (see [207, 208]). Then, noticing that $\nabla u^\infty + u^\infty \nabla \phi = \omega \nabla \rho^\infty - \rho^\infty \omega \nabla \phi + \rho^\infty \omega \nabla \phi = \omega \nabla \rho^\infty$, the Slotboom variable $\rho^\infty$ equivalently solves the following pure diffusion problem:

$$ \left\{ \begin{array}{ll} -\operatorname{div}(\omega \Lambda \nabla \rho^\infty) = f & \text{in } \Omega, \\ \rho^\infty = \omega^{-1} g^D & \text{on } \Gamma^D, \\ \omega \Lambda \nabla \rho^\infty \cdot n = g^N & \text{on } \Gamma^N, \end{array} \right. \quad (1.28) $$
with additional constraint $\int_{\Omega} \omega \rho^{\infty} = M$ when $|\Gamma^D| = 0$. Following Remark 2 (with $\rho$ instead of $u$, $\rho^D$ lifting of $\omega^{-1}g^D$ instead of $u^D$ lifting of $g^D$, and $\rho^M = \frac{M}{\int_{\Omega} \omega} \in \mathbb{R}$ instead of $u^M$), we consider Problem (1.28) with $\rho^{\infty} = 0$ on $\Gamma^D$ when $|\Gamma^D| > 0$, or $\int_{\Omega} \omega \rho^{\infty} = 0$ otherwise (which is equivalent to consider Problem (1.2) with $g^D = 0$ or $M = \int_{\Omega} u^{\infty} = 0$). Since $\phi$ is continuous on $\bar{\Omega}$, there exist $\omega_b, \omega_{\#}$ with $0 < \omega_b \le \omega_{\#} < \infty$, only depending on $\phi$ and $\Omega$, such that $\omega_b \le \omega(x) \le \omega_{\#}$ for all $x \in \bar{\Omega}$. We then denote by $L_{\omega}^2(\Omega)$ the $\omega$-weighted $L^2$ space on $\Omega$.

At the discrete level, instead of discretising (1.2), we approximate the solution to (1.28). For any $K \in \mathcal{M}$, we let $a_K^\omega: \underline{V}_K \times \underline{V}_K \to \mathbb{R}$ be the discrete bilinear form such that, for all $\underline{\rho}_K, \underline{v}_K \in \underline{V}_K$,

$$a_K^\omega(\underline{\rho}_K, \underline{v}_K) = (\omega \Lambda \nabla_K \underline{\rho}_K, \nabla_K \underline{v}_K)_K, \quad (1.29)$$

and, classically, we let $a_D^\omega: \underline{V}_D \times \underline{V}_D \to \mathbb{R}$ be the corresponding global discrete bilinear form obtained by sum of the local contributions. To account for the change of variable, we let $\underline{V}_{D,0}^\omega$ be the space $\underline{V}_{D,0}^D$ when $|\Gamma^D| > 0$, and the space $\{\underline{v}_D \in \underline{V}_D: \int_{\Omega} \omega v_M = 0\}$ otherwise. The discrete problem reads: Find $\underline{\rho}_D^\omega \in \underline{V}_{D,0}^\omega$ such that

$$a_D^\omega(\underline{\rho}_D^\omega, \underline{v}_D) = (f, v_M)_{\Omega} + (g^N, v_{\mathcal{E}})_{\Gamma^N} \quad \forall \underline{v}_D \in \underline{V}_{D,0}^\omega. \quad (1.30)$$

Remark that, for pure Neumann boundary conditions with $M \ne 0$, as for Problem (1.7), letting $\underline{\rho}_D^M = \rho^M \underline{1}_D$, $a_D^\omega(\underline{\rho}_D^M, \underline{v}_D) = 0$ for all $\underline{v}_D \in \underline{V}_D$. Letting $\underline{\omega}_D \in \underline{V}_D$ be the HFV interpolate of $\omega$, i.e.,

$$\omega_K = \frac{1}{|K|} \int_K \omega \quad \forall K \in \mathcal{M}, \quad \omega_\sigma = \frac{1}{|\sigma|} \int_\sigma \omega \quad \forall \sigma \in \mathcal{E}, \quad (1.31)$$

the approximation of the solution to Problem (1.2) is finally defined as the product $\underline{u}_D^\infty = \underline{\omega}_D \times \underline{\rho}_D^\infty$, that is

$$u_K^\infty = \omega_K \underline{\rho}_K^\infty \quad \forall K \in \mathcal{M}, \quad u_\sigma^\infty = \omega_\sigma \underline{\rho}_\sigma^\infty \quad \forall \sigma \in \mathcal{E}. \quad (1.32)$$

Remark that $\underline{u}_D^\infty \in \underline{V}_{D,0}$. Reasoning as in Section 1.2.3, one can easily prove that, for all $\underline{v}_D \in \underline{V}_D$,

$$a_D^\omega(\underline{v}_D, \underline{v}_D) \ge \omega_b \lambda_b \alpha_b |\underline{v}_D|^2_{1,D}. \quad (1.33)$$

This estimate is instrumental to infer well-posedness for Problem (1.30)-(1.32).

**Proposition 3** (Well-posedness). There exists a unique solution $\underline{u}_D^\infty \in \underline{V}_{D,0}$ to Problem (1.30)-(1.32), which satisfies $|\underline{u}_D^\infty|_{1,D} \le C (\|f\|_{L^2(\Omega)} + \|g^N\|_{L^2(\Gamma^N)})$, for some $C > 0$ depending on the data, and on the discretisation $\mathcal{D}$ only through $\theta_D$.

*Proof.* The existence/uniqueness of $\underline{\rho}_D^\infty \in \underline{V}_{D,0}^\omega$ solution to (1.30) (and in turn of $\underline{u}_D^\infty \in \underline{V}_{D,0}$) is a direct consequence of the coercivity estimate (1.33), and of the fact that $|\cdot|_{1,D}$ clearly defines a norm on $\underline{V}_{D,0}^\omega$. To prove the bound on $|\underline{u}_D^\infty|_{1,D}$, we use the fact that $\underline{u}_D^\infty = \underline{\omega}_D \times \underline{\rho}_D^\infty$. For all $K \in \mathcal{M}$,

$$|\underline{u}_K^\infty|_{1,K}^2 = \sum_{\sigma \in \mathcal{E}_K} \frac{|\sigma|}{d_{K,\sigma}} (\omega_K \underline{\rho}_K^\infty - \omega_\sigma \underline{\rho}_\sigma^\infty)^2 \le 2 \sum_{\sigma \in \mathcal{E}_K} \frac{|\sigma|}{d_{K,\sigma}} (\omega_\sigma^2 (\underline{\rho}_K^\infty - \underline{\rho}_\sigma^\infty)^2 + (\underline{\rho}_K^\infty)^2 (\omega_K - \omega_\sigma)^2).$$

By definition (1.31) of $\underline{\omega}_D$, and local stability of the HFV interpolant (cf. [111, Proposition B.7,
$p=2]$ and its proof), we infer

$$|\underline{u}_K^\infty|_{1,K}^2 \le 2 \omega_\#^2 |\underline{\rho}_K^\infty|_{1,K}^2 + 2 (\rho_K^\infty)^2 C_{\text{sta}}^2 \|\nabla \omega\|_{L^2(K;\mathbb{R}^d)}^2$$

with $C_{\text{sta}} > 0$ only depending on $d$ and $\theta_D$. Since $\nabla\omega = -\omega\nabla\phi$, $\phi \in C^1(\bar{\Omega})$, and $\omega > 0$, we have

$$\|\nabla \omega\|_{L^2(K;\mathbb{R}^d)}^2 \le \omega_\# \sup_{x \in \bar{\Omega}} |\nabla \phi(x)|^2 |K| \omega_K.$$

Summing over $K \in \mathcal{M}$ then yields

$$|\underline{u}_D^\infty|_{1,D}^2 \le 2 \omega_\#^2 |\underline{\rho}_D^\infty|_{1,D}^2 + 2 C_{\text{sta}}^2 \omega_\# \sup_{x \in \bar{\Omega}} |\nabla \phi(x)|^2 \|\rho_M^\infty\|_{L^2(\Omega)}^2.$$

When $|\Gamma^D| > 0$, $\|\rho_M^\infty\|_{L^2(\Omega)}^2 \le \omega_\# \|\rho_M^\infty\|_{L^2(\Omega)}^2 \le \omega_\# C_{P,\Gamma^D}^2 |\underline{\rho}_D^\infty|_{1,D}^2$, where we have applied the discrete Poincaré inequality (1.83) to $\underline{\rho}_D^\infty \in \underline{V}_{D,0}^D$. Otherwise, $\underline{\rho}_D^\infty \in \underline{V}_{D,0}^\omega$ satisfies $\int_\Omega \omega \rho_M^\infty = 0$, and one can use [51, Lemma 5.2] to infer that $\|\rho_M^\infty\|_{L^2(\Omega)} \le 2\|\rho_M^\infty - \frac{1}{|\Omega|} \int_\Omega \rho_M^\infty\|_{L^2(\Omega)}\|$, and finally get that $\|\rho_M^\infty\|_{L^2(\Omega)} \le 2 \sqrt{\omega_\#} C_{PW} |\underline{\rho}_D^\infty|_{1,D}$ applying the discrete Poincaré inequality (1.82) to $\underline{\rho}_D^\infty - \frac{1}{|\Omega|} \int_\Omega \rho_M^\infty \underline{1}_D \in \underline{V}_{D,0}^N$. In any case, we end up bounding $|\underline{u}_D^\infty|_{1,D}$ by $|\underline{\rho}_D^\infty|_{1,D}$, with multiplicative constant depending on the data, and on the discretisation $\mathcal{D}$ only through $\theta_D$. The rest of the proof consists in bounding $|\underline{\rho}_D^\infty|_{1,D}$, and proceeds as in the proof of Proposition 1, using that

$$\|\rho_M^\infty\|_{L^2(\Omega)} \le 2 \sqrt{\frac{\omega_\#}{\omega_\#}} C_{PW} |\underline{\rho}_D^\infty|_{1,D} \text{ when } |\Gamma^D| = 0. \quad \square$$

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

In the sequel, the HFV scheme (1.30)-(1.32) will be referred to as “exponential fitting”. Notice that no assumption on $V^\phi$ is needed to ensure its coercivity.

**Remark 6** (Preservation of the thermal equilibrium). As for the exponential fitting scheme of [43], the method (1.30)-(1.32) preserves the thermal equilibrium, in the sense that $\underline{u}_D^\infty = \bar{\rho}\underline{\omega}_D$ when $u^\infty = u_{th}^\infty$ (see (1.25)). This property is analogous to what holds true for the TPFA Scharfetter-Gummel scheme.

## Evolution problem

We consider Problem (1.1). Following the previous strategy, letting $\rho = \frac{u}{\omega}$, one can show that $\rho$ equivalently solves the following transient pure diffusion problem:

$$\left\{
\begin{array}{ll}
\omega \partial_t \rho - \operatorname{div}(\omega \Lambda \nabla \rho) = f & \text{in } \mathbb{R}_+ \times \Omega, \\
\rho = \omega^{-1} g^D & \text{on } \mathbb{R}_+ \times \Gamma^D, \\
\omega \Lambda \nabla \rho \cdot n = g^N & \text{on } \mathbb{R}_+ \times \Gamma^N, \\
\rho(0, \cdot) = \rho^{in} & \text{in } \Omega,
\end{array}
\right.
\quad (1.34)$$

where $\rho^{in} = \frac{u^{in}}{\omega}$. Following Remark 2, we consider Problem (1.34) with $\rho = 0$ on $\Gamma^D$ when $|\Gamma^D| > 0$, or $\int_\Omega \omega \rho^{in} = 0$ otherwise (which is equivalent to consider Problem (1.1) with $g^D = 0$ or $M = \int_\Omega u^{in} = 0$).

At the discrete level, instead of discretising (1.1), we approximate the solution to (1.34). We use a backward Euler discretisation in time, and the HFV discretisation introduced in
Section 1.3.2 in space. The discrete problem reads: Find $(\rho_D^n \in \underline{V}_{D,0}^\omega)_{n \ge 1}$ such that

$$ \left\{ \begin{array}{ll} \frac{1}{\Delta t}(\rho_M^n - \rho_M^{n-1}, \omega_M v_M)_\Omega + a_D^\omega(\rho_D^n, \underline{v}_D) = (f, v_M)_\Omega + (g^N, v_\varepsilon)_{\Gamma^N} & \forall \underline{v}_D \in \underline{V}_{D,0}^\omega, \quad (1.35a) \\ \rho_K^0 = \frac{1}{\omega_K |K|} \int_K u^{in} & \forall K \in \mathcal{M}, \quad (1.35b) \end{array} \right. $$

where $a_D^\omega$ is defined (locally) by (1.29), and the approximation $\underline{u}_D^n \in \underline{V}_{D,0}$ of the solution to Problem (1.1) is finally defined as $\underline{u}_D^n = \underline{\omega}_D \times \rho_D^n$, i.e., according to (1.32) (with superscript $n$ instead of $\infty$). Once again, because of the coercivity of $a_D^\omega$, the scheme (1.35) is unconditionally well-posed, as a straightforward consequence of Proposition 3.

**Remark 7 (Pure Neumann case).** When considering pure Neumann boundary conditions, as for the standard HFV scheme (see Remark 5), one can seek for a solution to Problem (1.35) in $\underline{V}_D$.

### 1.3.3 Nonlinear HFV scheme

We are interested in the evolution problem (1.1). We restrict our study to the pure Neumann case ($|\Gamma^D| = 0$), and to the choice of data $f = 0$, $g^N = 0$, and $u^{in} \ge 0$ with $M = \int_\Omega u^{in} > 0$ (see Appendix 1.B for the case of mixed Dirichlet-Neumann thermal equilibrium boundary conditions). Under these assumptions, it is known that the solution $u$ to Problem (1.1) is strictly positive on $\mathbb{R}_+^* \times \Omega$. Furthermore, in long time, $u(t)$ converges towards the thermal equilibrium. Indeed, we easily verify that the function $u^\infty = \rho^M e^{-\phi}$, where we recall that $\rho^M = \frac{M}{\int_\Omega e^{-\phi}} > 0$ (cf. Section 1.3.2), solves the steady problem (1.2) with same data. Since $u > 0$ on $\mathbb{R}_+^* \times \Omega$, we can rewrite the flux $J = -\Lambda(\nabla u + u\nabla\phi)$ under the nonlinear form

$$ J = -u\Lambda\nabla(\log(u) + \phi) = -u\Lambda\nabla\log\left(\frac{u}{u^\infty}\right). $$

At the continuous level, introducing this nonlinearity enables to highlight the following entropy/dissipation structure of the model at hand: testing the equation against $\log(\frac{u}{u^\infty})$, we get

$$ \frac{d}{dt} \mathbb{E}(t) + \mathbb{D}(t) = 0, \quad (1.36) $$

where

$$ \mathbb{E}(t) = \int_{\Omega} u^{\infty} \Phi_1 \left( \frac{u(t)}{u^{\infty}} \right) \quad \text{and} \quad \mathbb{D}(t) = \int_{\Omega} u(t) \Lambda \nabla \log \left( \frac{u(t)}{u^{\infty}} \right) \cdot \nabla \log \left( \frac{u(t)}{u^{\infty}} \right), \quad (1.37) $$

with $\Phi_1(s) = s \log(s) - s + 1$ for all $s > 0$. Since $\Phi_1 \ge 0$, the relative entropy $\mathbb{E}(t)$ is a non-negative quantity (as well as the relative dissipation $\mathbb{D}(t)$). The entropy/dissipation structure (1.36)-(1.37) is instrumental to prove the exponential convergence in time of the solution $u(t)$ to Problem (1.1) towards the equilibrium $u^\infty$. From the above nonlinear expression of the flux $J$, we build a nonlinear hybrid discretisation of the problem, leading to a scheme designed along the same principles as the nonlinear VAG and DDFV schemes of [56] and [52, 51]. This scheme is devised so as to ensure the positivity of discrete solutions, as well as to preserve at the discrete level the entropy/dissipation structure (and the long-time behaviour) of the model. The choice of designing a nonlinear HFV scheme is driven by the prospect of the design of hybrid high-order (HHO) schemes which could have similar features.
Definition of the scheme and key properties of discrete solutions

In the sequel, a vector of discrete unknowns $\underline{v}_D \in \underline{V}_D$ will be called positive if and only if, for all $K \in \mathcal{M}$ and all $\sigma \in \mathcal{E}$, $v_K > 0$ and $v_\sigma > 0$. Recall the definition (1.27) of $\omega = e^{-\phi}$, as well as the definition (1.31) of the HFV interpolate $\underline{\omega}_D \in \underline{V}_D$ of $\omega$. Remark that $\underline{\omega}_D$ is positive. If $\underline{u}_D \in \underline{V}_D$ is positive, one can then define $\underline{w}_D$ as the element of $\underline{V}_D$ such that

$$w_K = \log\left(\frac{u_K}{\omega_K}\right) \quad \forall K \in \mathcal{M}, \qquad w_\sigma = \log\left(\frac{u_\sigma}{\omega_\sigma}\right) \quad \forall \sigma \in \mathcal{E}. \qquad (1.38)$$

In what follows, to emphasise the dependency of $\underline{w}_D$ upon $\underline{u}_D$, we sometimes write $\underline{w}_D(\underline{u}_D)$. Locally to any cell $K \in \mathcal{M}$, we define an approximation of

$$(u, v) \mapsto - \int_K J \cdot \nabla v = \int_K u \Lambda \nabla \log \left( \frac{u}{u^\infty} \right) \cdot \nabla v$$

under the form

$$T_K(\underline{u}_K, \underline{w}_K, \underline{v}_K) = \int_K r_K(\underline{u}_K) \Lambda \nabla_K \underline{w}_K \cdot \nabla_K \underline{v}_K,$$

for all $\underline{u}_K \in \underline{V}_K$ positive and all $\underline{v}_K \in \underline{V}_K$, where $r_K : (\underline{V}_K)_+^\star \to \mathbb{R}_+^\star$ is a local reconstruction operator. Since $r_K(\underline{u}_K)$ is a (positive) constant on $K$, we have

$$T_K(\underline{u}_K, \underline{w}_K, \underline{v}_K) = r_K(\underline{u}_K) a_K^\Lambda(\underline{w}_K, \underline{v}_K), \qquad (1.39)$$

where $a_K^\Lambda$ is defined by (1.6). Following (1.12), one can equivalently reformulate (1.39) using the local matrix $A_K$ defined by (1.10):

$$T_K(\underline{u}_K, \underline{w}_K, \underline{v}_K) = r_K(\underline{u}_K) \delta_K \underline{v}_K \cdot A_K \delta_K \underline{w}_K. \qquad (1.40)$$

As already pointed out in the analysis of the nonlinear DDFV scheme of [52, 51], the definition of the local reconstruction operator is crucial to guarantee the existence of solutions and a good long-time behaviour to the scheme. The most natural choice in the HFV context would obviously be $r_K(\underline{u}_K) = u_K$, however it turns out that such a reconstruction embeds too few information on $\underline{u}_K$ to conclude, as already suggested in [49]. Therefore, we use a richer reconstruction, described below, which embeds information from both the local cell and face unknowns. For $\underline{u}_K \in \underline{V}_K$ positive, we let

$$r_K(\underline{u}_K) = f_{\varepsilon_K} \left| \left( (m(u_K, u_\sigma))_{\sigma \in \varepsilon_K} \right) \right|, \qquad (1.41)$$

with $m : (\mathbb{R}_+^\star)^2 \to \mathbb{R}_+^\star$ and, for $k \ge 1$ integer, $f_k : (\mathbb{R}_+^\star)^k \to \mathbb{R}_+^\star$, such that

$$m \text{ is non-decreasing with respect to both its variables,} \qquad (1.42a)$$

$$m(x, x) = x \text{ for all } x \in \mathbb{R}_+^\star \text{ and } m(y, x) = m(x, y) \text{ for all } (x, y) \in (\mathbb{R}_+^\star)^2, \qquad (1.42b)$$

$$m(\lambda x, \lambda y) = \lambda m(x, y) \text{ for all } \lambda > 0 \text{ and all } (x, y) \in (\mathbb{R}_+^\star)^2, \qquad (1.42c)$$

$$\frac{y-x}{\log(y) - \log(x)} \le m(x, y) \le \max(x, y) \text{ for all } (x, y) \in (\mathbb{R}_+^\star)^2, x \ne y, \qquad (1.42d)$$
and

$$f_k(x_1, \dots, x_k) = \frac{1}{k} \sum_{i=1}^{k} x_i \quad \text{or} \quad f_k(x_1, \dots, x_k) = \max(x_1, \dots, x_k). \qquad (1.43)$$

Note that, for all $(x, y) \in (\mathbb{R}_+^*)^2$, one has

$$\frac{y-x}{\log(y)-\log(x)} \le \left( \frac{\sqrt{x}+\sqrt{y}}{2} \right)^2 \le \frac{x+y}{2} \le \max(x,y),$$

and each expression of the previous sequence is a mean function $m$ satisfying the properties (1.42). Heuristically, $r_K(u_K)$ computes an average of the unknowns attached to the cell $K$, especially it contains information about all the local face unknowns. As far as the properties (1.42)-(1.43) are concerned, they will be instrumental to prove Lemma 2 and Proposition 6 below. As now standard, we finally let $T_D$ be such that, for all $\underline{u}_D \in \underline{V}_D$ positive, and all $\underline{v}_D \in \underline{V}_D$,

$$T_D(\underline{u}_D, \underline{w}_D, \underline{v}_D) = \sum_{K \in \mathcal{M}} T_K(\underline{u}_K, \underline{w}_K, \underline{v}_K), \qquad (1.44)$$

where the local contributions $T_K$ are defined by (1.39).

Using a backward Euler discretisation in time, and the HFV discretisation we have just introduced in space, our discrete problem reads: Find $(\underline{u}_D^n \in \underline{V}_D)_{n \ge 1}$ such that

$$\left\{ \frac{1}{\Delta t} (u_M^n - u_M^{n-1}, v_M)_{\Omega} + T_D(\underline{u}_D^n, \underline{w}_D(\underline{u}_D^n), \underline{v}_D) = 0 \quad \forall \underline{v}_D \in \underline{V}_D, \right. \qquad (1.45a)$$

$$u_K^0 = \frac{1}{|K|} \int_K u^{in} \quad \forall K \in \mathcal{M}. \qquad (1.45b)$$

Notice that if $(\underline{u}_D^n)_{n \ge 1}$ solves Problem (1.45), then, necessarily, $\underline{u}_D^n$ is positive for all $n \ge 1$. Therefore, in the sequel, we will speak about the positive solutions to (1.45). Notice also that $u_M^0$ may vanish in some cells of the mesh, since we only impose that $u^{in} \ge 0$ (but $u_M^0$ cannot be identically zero in $\Omega$ since $M > 0$). Notice finally that $u_E^0$ needs not be defined, as the scheme only uses $u_M^0$.

Testing (1.45a) with $\underline{v}_D = \underline{1}_D$, and remarking that $T_D(\underline{u}_D^n, \underline{w}_D(\underline{u}_D^n), \underline{1}_D) = 0$ for all $n \ge 1$, we immediately infer the following discrete mass conservation property.

**Proposition 4 (Mass conservation).** If $(\underline{u}_D^n \in \underline{V}_D)_{n \ge 1}$ is a (positive) solution to (1.45), then

$$\forall n \in \mathbb{N}^*, \quad \int_{\Omega} u_M^n = \int_{\Omega} u_M^0 = \int_{\Omega} u^{in} = M.$$

Following Proposition 4, a discrete steady-state $\underline{u}_D^\infty \in \underline{V}_D$ of (1.45) shall satisfy

$$T_D(\underline{u}_D^\infty, \underline{w}_D(\underline{u}_D^\infty), \underline{v}_D) = 0 \quad \forall \underline{v}_D \in \underline{V}_D, \qquad (1.46)$$

and $\int_{\Omega} u_M^\infty = M$. Letting $\underline{w}_D^\infty = \underline{w}_D(\underline{u}_D^\infty)$, and testing (1.46) with $\underline{v}_D = \underline{w}_D^\infty$, by (1.39) and (1.44), since $r_K(\underline{u}_K^\infty) > 0$ for all $K \in \mathcal{M}$, we necessarily have $a_D^\Lambda(\underline{w}_D^\infty, \underline{w}_D^\infty) = 0$, which yields $|\underline{w}_D^\infty|_{1,D} = 0$ by the coercivity property (1.14). Hence, $\underline{w}_D^\infty = c \underline{1}_D$ for some constant $c \in \mathbb{R}$, and since $\int_{\Omega} u_M^\infty = M$, by (1.38), we necessarily have $e^c = \rho^M$, that is $c = \log(\rho^M)$ and $\underline{w}_D^\infty = \log(\rho^M) \underline{1}_D$. As a consequence, again by (1.38), $\underline{u}_D^\infty = \rho^M \underline{w}_D$, i.e., $\underline{u}_D^\infty$ is the HFV interpolate of $u_{th}^\infty$. Thus, just like the exponential
fitting scheme (cf. Remark 6), the nonlinear scheme preserves the thermal equilibrium. We notice that $w_D$, first defined by (1.38), can actually be modified up to an additive constant without any impact on the scheme (1.45). Hence, we can redefine $w_D$ as

$$w_K = \log\left(\frac{u_K}{u_K^\infty}\right) \quad \forall K \in \mathcal{M}, \quad w_\sigma = \log\left(\frac{u_\sigma}{u_\sigma^\infty}\right) \quad \forall \sigma \in \mathcal{E}. \qquad (1.47)$$

Another important consequence of the fact that $u_D^\infty$ is the HFV interpolate of $u_{th}^\infty$ is the following. Letting $u_b^\infty = \frac{M\omega_b}{|\Omega|\omega_\#} > 0$ and $u_\#^\infty = \frac{M\omega_\#}{|\Omega|\omega_b} > 0$ (recall that $\omega_b$ and $\omega_\#$ only depend on $\phi$ and $\Omega$), we have $u_b^\infty \le u_{th}^\infty \le u_\#^\infty$, but we also have that

$$u_b^\infty \underline{1}_D \le u_D^\infty \le u_\#^\infty \underline{1}_D, \qquad (1.48)$$

where the inequalities shall be understood coordinate-wise. In other words, the continuous bounds on the steady-state are transferred to the discrete level.

Given a (positive) solution $(\underline{u}_D^n \in \underline{V}_D)_{n \ge 1}$ to (1.45), we define the following discrete versions of the relative entropy and dissipation introduced in (1.37): for all $n \ge 1$,

$$\mathbb{E}^n = \int_{\Omega} u_M^\infty \Phi_1 \left( \frac{u_M^n}{u_M^\infty} \right) \quad \text{and} \quad \mathbb{D}^n = T_D(\underline{u}_D^n, \underline{w}_D^n, \underline{w}_D^n), \qquad (1.49)$$

where we let $\underline{w}_D^n = \underline{w}_D(\underline{u}_D^n)$, and we recall that $\Phi_1(s) = s \log(s) - s + 1$ for all $s > 0$. Notice that $\mathbb{E}^n \ge 0$ for all $n \ge 1$ since $\Phi_1 \ge 0$. For further use, we extend the function $\Phi_1$ by continuity to 0, letting $\Phi_1(0) = 1$. We can then define, in case there exists $K \in \mathcal{M}$ such that $u_K^0 = 0$, $\mathbb{E}^0 \ge 0$ according to (1.49). As far as $\mathbb{D}^n$ is concerned, by (1.39) and (1.44), for all $n \ge 1$, we have

$$\mathbb{D}^n = \sum_{K \in \mathcal{M}} r_K(\underline{u}_K^n) a_K^\Lambda(\underline{w}_K^n, \underline{w}_K^n) \ge 0.$$

We can now establish the following discrete counterpart of (1.36).

**Proposition 5 (Entropy dissipation).** If $(\underline{u}_D^n \in \underline{V}_D)_{n \ge 1}$ is a (positive) solution to (1.45), then

$$\forall n \in \mathbb{N}, \quad \frac{\mathbb{E}^{n+1} - \mathbb{E}^n}{\Delta t} + \mathbb{D}^{n+1} \le 0. \qquad (1.50)$$

*Proof.* Let $n \in \mathbb{N}$. By the expression (1.49) of the discrete relative entropy, and the convexity of $\Phi_1$, we have

$$\mathbb{E}^{n+1} - \mathbb{E}^n \le \sum_{K \in \mathcal{M}} |K| u_K^\infty \Phi_1' \left( \frac{u_K^{n+1}}{u_K^\infty} \right) \left( \frac{u_K^{n+1} - u_K^n}{u_K^\infty} \right).$$

Thus, by (1.47), we get

$$\mathbb{E}^{n+1} - \mathbb{E}^n \le \int_{\Omega} (u_M^{n+1} - u_M^n) \log \left( \frac{u_M^{n+1}}{u_M^n} \right) = (u_M^{n+1} - u_M^n, w_M^{n+1})_{\Omega}.$$

By (1.45a) and (1.49), we finally infer

$$\mathbb{E}^{n+1} - \mathbb{E}^n \le -\Delta t T_D(\underline{u}_D^{n+1}, \underline{w}_D^{n+1}, \underline{w}_D^{n+1}) = -\Delta t \mathbb{D}^{n+1},$$