# A posteriori error estimates and adaptivity for locally conservative methods. Inexpensive implementation and evaluation, polytopal meshes, iterative linearization and algebraic solvers, and applications to complex porous media flows

Martin Vohralík, Soleiman Yousef

► To cite this version:

Martin Vohralík, Soleiman Yousef. A posteriori error estimates and adaptivity for locally conservative methods. Inexpensive implementation and evaluation, polytopal meshes, iterative linearization and algebraic solvers, and applications to complex porous media flows. 2025. hal-05081807v3

HAL Id: hal-05081807

https://inria.hal.science/hal-05081807v3

Preprint submitted on 14 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
2.3 Numerical approximation: spatial discretization, temporal discretization, iterative linearization, and iterative linear algebraic resolution 6

2.4 Total error and its components 6

2.4.1 Intrinsic error measure 6

2.4.2 Spatial discretization error 6

2.4.3 Iterative linear algebraic solver error 7

2.4.4 Iterative linearization error 7

2.4.5 Temporal discretization error 7

2.4.6 Error components 7

2.5 A posteriori error estimate 8

2.5.1 Guaranteed error upper bound (reliability) 8

2.5.2 Error lower bound (efficiency) 8

2.5.3 Robustness with respect to the spatial and temporal domains as well as physical and numerical parameters 9

2.5.4 Asymptotic exactness 9

2.5.5 Inexpensive evaluation 9

2.5.6 Estimating the error components 9

2.6 Adaptivity: balancing the error component estimates 9

2.7 Mass balance recovery at any step of numerical resolution 10

2.8 Bibliographic resources 10

3 Sobolev spaces and their piecewise polynomial subspaces, meshes, and finite volume methods 11

3.1 Basic notation 11

3.1.1 Domain Ω 11

3.1.2 Lebesgue spaces 11

3.1.3 Measures, sizes, and cardinalities 11

3.2 Infinite-dimensional Sobolev spaces $H_0^1(\Omega)$ and $H(\operatorname{div}, \Omega)$ 11

3.2.1 Sobolev space $H_0^1(\Omega)$ 11

3.2.2 Sobolev space $H(\operatorname{div}, \Omega)$ 12

3.3 Simplicial meshes 13

3.4 Broken Sobolev spaces $H^1(\mathcal{T}_h)$ and $H(\operatorname{div}, \mathcal{T}_h)$ 14

3.5 Finite-dimensional subspaces of $H_0^1(\Omega)$ and $H(\operatorname{div}, \Omega)$ 14

3.5.1 Piecewise polynomial subspaces of $L^2(\Omega)$ 14

3.5.2 Lagrange piecewise polynomial subspaces of $H_0^1(\Omega)$ 15

3.5.3 Raviart-Thomas piecewise polynomial subspaces of $H(\operatorname{div}, \Omega)$ 15

3.6 Polytopal meshes with virtual simplicial submeshes 17

3.6.1 Polytopal mesh with a virtual simplicial submesh 17

3.6.2 Mesh faces and mesh vertices 17

3.7 Basic principle of finite volume and related lowest-order locally conservative methods 18

3.8 Bibliographic resources 19

4 The Poisson equation. Simplicial meshes and spatial discretization error 19

4.1 The Poisson equation 20

4.2 Weak solution, potential, and flux 20

4.3 Properties of the exact potential and flux 20

4.4 Approximate solution and spatial discretization error 21

4.5 Prager-Synge equality 21

4.6 Equivalence of the spatial discretization error to the flux nonconformity (dual norm of the residual) plus potential nonconformity 22

4.7 Cell-centered finite volume discretizations on a simplicial mesh 24

4.8 Flux reconstruction by lifting the face normal fluxes 24

4.9 Potential reconstruction by elementwise postprocessing and averaging 25

4.10 A guaranteed a posteriori error estimate 26

4.11 Numerical experiments 27

4.11.1 Regular solution 27

4.11.2 Singular solution 28

4.12 Bibliographic resources 28
5 Steady linear pure diffusion problems. Polytopal meshes and inexpensive implemen-
tation and evaluation 29

<table>
  <tbody>
    <tr>
      <td>5.1</td>
      <td>Singlephase steady linear Darcy flow</td>
      <td>30</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Weak solution and its properties</td>
      <td>30</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>Spatial discretization error</td>
      <td>30</td>
    </tr>
    <tr>
      <td>5.4</td>
      <td>Generic discretizations on a polytopal mesh</td>
      <td>30</td>
    </tr>
    <tr>
      <td>5.5</td>
      <td>Face normal fluxes</td>
      <td>31</td>
    </tr>
    <tr>
      <td>5.6</td>
      <td>Fictitious flux reconstruction by lifting the face normal fluxes</td>
      <td>32</td>
    </tr>
    <tr>
      <td>5.6.1</td>
      <td>A (fictitious) flux reconstruction under Assumption 5.1</td>
      <td>32</td>
    </tr>
    <tr>
      <td>5.6.2</td>
      <td>A fictitious flux reconstruction under Assumption 5.2</td>
      <td>34</td>
    </tr>
    <tr>
      <td>5.7</td>
      <td>Potential point values</td>
      <td>35</td>
    </tr>
    <tr>
      <td>5.7.1</td>
      <td>Potential point values under Assumption 5.1</td>
      <td>35</td>
    </tr>
    <tr>
      <td>5.7.2</td>
      <td>Potential point values under Assumption 5.2</td>
      <td>36</td>
    </tr>
    <tr>
      <td>5.8</td>
      <td>Fictitious potential reconstruction by respecting the point values</td>
      <td>37</td>
    </tr>
    <tr>
      <td>5.9</td>
      <td>A guaranteed a posteriori error estimate with inexpensive implementation and evaluation</td>
      <td>38</td>
    </tr>
    <tr>
      <td>5.10</td>
      <td>Numerical experiments</td>
      <td>40</td>
    </tr>
    <tr>
      <td>5.11</td>
      <td>Bibliographic resources</td>
      <td>42</td>
    </tr>
  </tbody>
</table>

6 Steady nonlinear pure diffusion problems. Iterative linearization and algebraic errors 42

<table>
  <tbody>
    <tr>
      <td>6.1</td>
      <td>Singlephase steady nonlinear Darcy flow</td>
      <td>42</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Weak solution and its properties</td>
      <td>43</td>
    </tr>
    <tr>
      <td>6.3</td>
      <td>Total error and its spatial discretization, linearization, and algebraic components</td>
      <td>43</td>
    </tr>
    <tr>
      <td>6.4</td>
      <td>Generic discretizations on a polytopal mesh, iterative linearization, and iterative algebraic system solution</td>
      <td>44</td>
    </tr>
    <tr>
      <td>6.5</td>
      <td>Face normal fluxes and potential point values</td>
      <td>45</td>
    </tr>
    <tr>
      <td>6.6</td>
      <td>Fictitious flux and potential reconstructions</td>
      <td>46</td>
    </tr>
    <tr>
      <td>6.7</td>
      <td>A guaranteed a posteriori error estimate distinguishing the error components with inex- pensive implementation and evaluation</td>
      <td>47</td>
    </tr>
    <tr>
      <td>6.8</td>
      <td>Numerical experiments</td>
      <td>50</td>
    </tr>
    <tr>
      <td>6.8.1</td>
      <td>Regular solution</td>
      <td>50</td>
    </tr>
    <tr>
      <td>6.8.2</td>
      <td>Singular solution</td>
      <td>52</td>
    </tr>
    <tr>
      <td>6.9</td>
      <td>Bibliographic resources</td>
      <td>53</td>
    </tr>
  </tbody>
</table>

7 Unsteady nonlinear coupled degenerate advection–diffusion–reaction problems. Temporal discretization error 54

<table>
  <tbody>
    <tr>
      <td>7.1</td>
      <td>Multiphase compositional unsteady Darcy flow</td>
      <td>55</td>
    </tr>
    <tr>
      <td>7.2</td>
      <td>Space-time mesh</td>
      <td>56</td>
    </tr>
    <tr>
      <td>7.3</td>
      <td>Weak solution and its properties</td>
      <td>56</td>
    </tr>
    <tr>
      <td>7.4</td>
      <td>Total error and its spatial discretization, temporal discretization, linearization, and algebraic components</td>
      <td>57</td>
    </tr>
    <tr>
      <td>7.5</td>
      <td>Backward Euler time discretization and locally conservative space discretization on polytopal meshes</td>
      <td>57</td>
    </tr>
    <tr>
      <td>7.6</td>
      <td>Iterative linearization</td>
      <td>58</td>
    </tr>
    <tr>
      <td>7.7</td>
      <td>Iterative algebraic system solution</td>
      <td>59</td>
    </tr>
    <tr>
      <td>7.8</td>
      <td>Face normal fluxes and potential point values</td>
      <td>59</td>
    </tr>
    <tr>
      <td>7.9</td>
      <td>Fictitious flux and potential reconstructions</td>
      <td>60</td>
    </tr>
    <tr>
      <td>7.10</td>
      <td>A guaranteed a posteriori error estimate distinguishing the error components with inexpensive implementation and evaluation</td>
      <td>61</td>
    </tr>
    <tr>
      <td>7.11</td>
      <td>Numerical experiments</td>
      <td>63</td>
    </tr>
    <tr>
      <td>7.11.1</td>
      <td>Setup</td>
      <td>64</td>
    </tr>
    <tr>
      <td>7.11.2</td>
      <td>Two-phase Darcy flow</td>
      <td>64</td>
    </tr>
    <tr>
      <td>7.11.3</td>
      <td>Three-phases, three-components Darcy flow</td>
      <td>67</td>
    </tr>
    <tr>
      <td>7.12</td>
      <td>Bibliographic resources</td>
      <td>70</td>
    </tr>
  </tbody>
</table>
# 1 Introduction

The purpose of this contribution is to present the theory of a posteriori error estimates and adaptivity in numerical approximation of partial differential equations (PDEs) by locally conservative methods of the finite volume type. We focus on estimates that certify the error and can be easily coded, cheaply evaluated, and efficiently used, even in complex practical simulations. We consider polytopal meshes, i.e., meshes formed by general polygonal or polyhedral elements, appearing in various applications. We cover the overall chain of computational practice, including iterative linearization of nonlinear problems and algebraic resolution of systems of linear algebraic equations.

An optimal a posteriori estimate can be described by a set of requested properties, including namely the bound on the error between the unknown solution of the partial differential equation and the available numerical approximation that is guaranteed and fully computable from the approximate solution. We introduce this concept in Section 2 at an abstract level. We start from the notion of a PDE, discuss some intrinsic physical and mathematical properties of the weak solutions of PDEs, the notion of a total error and its components, and of adaptivity which aims at balancing these error components. One particular property of the methodology we develop here is that it leads to (exact or approximate) mass balance recovery at any step of numerical resolution, namely during the iteration of the nonlinear and linear solvers.

This contribution is made as self-standing as possible. For this purpose, in Section 3, we set up the notation, recall the intrinsic function spaces related to the considered partial differential equations, namely the infinite-dimensional Sobolev spaces $H_0^1(\Omega)$ and $H(\operatorname{div}, \Omega)$, and recall their finite-dimensional subspaces composed of piecewise polynomials with respect to a simplicial mesh. We also describe there the polytopal meshes we consider and introduce the basic principle of the finite volume and related lowest-order locally conservative methods.

The heart of our exposition starts in Section 4, where we treat in detail the model steady linear diffusion problem with homogeneous and isotropic diffusion tensor (the Poisson equation). Only the basic cell-centered finite volume discretization on a simplicial mesh and only the spatial discretization error are investigated here, to expose the construction principles as clearly as possible. We namely motivate and present the methodology based on $H^1(\Omega)$-conforming potential and $H(\operatorname{div}, \Omega)$-conforming equilibrated flux reconstructions, which yields a guaranteed and fully computable error upper bound. These reconstructions are practically obtained here in piecewise polynomial spaces. This is perfectly fine for academic purposes and has an asymptotically optimal evaluation cost, since it is in particular explicit and no local problems need to be solved. A cheaper implementation and evaluation of a posteriori error estimates is, however, possible and suitable for more complicated model or real-life problems, which we treat later. Numerical experiments illustrate these basic developments.

In Section 5, we extend the above analysis to tensor-valued steady linear diffusion problems (the singlephase steady linear Darcy flow), general polytopal meshes, and any lowest-order locally conservative method such as mixed finite elements, mixed and hybrid finite volumes, mimetic finite differences, mixed virtual finite elements, or hybrid high-order methods. In particular, though we still rely on the notions of potential and equilibrated flux reconstructions, these become virtual here, since they are not constructed in practice anymore. Similarly, though we rely on a notion of a simplicial submesh of a general polytopal mesh, this is also virtual, where the physical construction and computer implementation of a simplicial submesh is avoided. In particular, the evaluation of our estimates merely consists in some local matrix-vector multiplications, where, on each mesh element, the matrices are either directly inherited from the given numerical method, or easily constructed from the element geometry, while the vectors are the algebraic unknowns of the flux and potential approximations on the given element. This gives an easy and practically accessible application to polytopal meshes of the general methodology of $H^1$-conforming potential reconstruction and $H(\operatorname{div}, \Omega)$-conforming flux reconstruction. Numerical experiments focus on general polygonal meshes and illustrate the performance of this methodology.

In Section 6, we extend the previous developments to encompass two cornerstones of numerical approximation in applications: the nonlinear and linear solvers. Our problem here is the singlephase steady nonlinear Darcy flow, where an iterative linearization is applied and where the algebraic resolution of the system of arising linear algebraic equations is addressed. We still derive a guaranteed upper bound on the total error. Moreover, we identify the three arising error components, related to discretization, linearization, and algebraic solver, and design adaptive stopping criteria for the iterative solvers. Remarkably, the estimators still take the form of a simple matrix-vector multiplication, with the same local matrices as in the linear case. Numerical illustrations confirm the theoretical findings.

In the last part of this contribution, Section 7, we consider an unsteady nonlinear coupled degenerate
advection-diffusion-reaction problem, the multiphase compositional unsteady Darcy flow. We apply the entire methodology developed above to obtain estimates that are still guaranteed and fully computable for the intrinsic error measure, apply to any lowest-order locally conservative method on a polytopal mesh, and take a simple matrix-vector multiplication form. Moreover, they are valid on each stage of the overall solution algorithm: on each time step n, each linearization step k, and each linear solver step i. They also allow to distinguish different error components and design adaptive stopping criteria for the solvers, as well as adaptive choice of space and time meshes. A comprehensive numerical study on realistic porous media problems is included.

The present a posteriori error estimates can be readily implemented into numerical codes on general polytopal meshes, with a minimal overhead. These estimates allow for a very fast evaluation and, according to the main results summarized in Theorems 4.18, 5.18, 6.2, and 7.6, give a guaranteed control over the intrinsic error committed in the numerical approximation. Additionally, all the different error components (time and space discretizations, linearization, algebraic) are identified, leading to fully adaptive algorithms with all adaptive stopping criteria for linear and nonlinear solvers, adaptive time step management, and adaptive mesh refinement. Numerical experiments on real-life problems confirm important computational speed-ups that can be achieved with our methodology, in addition to the certification of the computed output.

Our presentation proceeds along the lines of Vohralík [249], Ern and Vohralík [124], and Vohralík and Yousef [255], using some central concepts for general polytopal meshes from Vohralík and Wohlmuth [254]. Ample bibliographic references are given in the dedicated sections below.

## 2 A posteriori error estimation and adaptivity – the principles

In this section, we set up the main ideas. We proceed informally; notation is in detail set in Section 3.

### 2.1 Partial differential equations

Let $\Omega \subset \mathbb{R}^d$, $d \ge 1$, be a polytopal domain where our problem is set. Consider an abstract (steady but possibly nonlinear) partial differential equation in the form: for a source term $f$, find a real-valued function $p: \Omega \to \mathbb{R}$ such that

$$ \mathcal{A}(p) = f \quad \text{in } \Omega, \qquad (2.1a) $$

$$ p = 0 \quad \text{on } \partial\Omega. \qquad (2.1b) $$

The prominent example that we will consider below is the Laplace operator

$$ \mathcal{A}(p) = -\nabla \cdot (\nabla p) = -\Delta p = -\sum_{i=1}^{d} \partial_{x_i}^2 p. \qquad (2.2) $$

In general, $\mathcal{A}(\cdot)$ is a (nonlinear) operator including partial derivatives of the unknown solution $p$.

More generally, for a final time $t_F > 0$, a source term $f$, and an initial condition $p_0$, consider an abstract (unsteady nonlinear) partial differential equation in the form: find a real-valued function $p: \Omega \times (0, t_F) \to \mathbb{R}$ such that

$$ \mathcal{A}(p) = f \quad \text{in } \Omega \times (0, t_F), \qquad (2.3a) $$

$$ p = 0 \quad \text{on } \partial\Omega \times (0, t_F), \qquad (2.3b) $$

$$ p(\cdot, 0) = p_0 \quad \text{in } \Omega. \qquad (2.3c) $$

Here, $\mathcal{A}(\cdot)$ is a (nonlinear) operator including partial derivatives of $p$ and namely a partial derivative with respect to the time, for example

$$ \mathcal{A}(p) = \partial_t p - \nabla \cdot (\underline{\mathbf{K}}(|\nabla p|)\nabla p), \qquad (2.4) $$

where $|\nabla p|$ is the Euclidean size of the vector $\nabla p$ and $\underline{\mathbf{K}}(\cdot)$ is a scalar- or matrix-valued nonlinear function of a real variable.

### 2.2 Physical and mathematical properties of the weak solutions of PDEs

In order to properly describe the solution $p$ of (2.1) or (2.3), adequate mathematical function spaces are necessary. We develop this task in details in Section 3.2 below, introducing namely the Sobolev spaces
$H_0^1(\Omega)$ and $\mathbf{H}(\mathrm{div}, \Omega)$. These spaces also reflect the physical principles behind problems (2.1) and (2.3). For instance, for the Laplace operator (2.2) in (2.1), it turns out that the solution $p$ satisfies

$$p \in H_0^1(\Omega), \quad \mathbf{u} := -\nabla p \in \mathbf{H}(\mathrm{div}, \Omega), \quad \nabla \cdot \mathbf{u} = f. \qquad (2.5)$$

This reflects that the scalar-valued field $p$ (pressure, potential, primal variable) is continuous in the appropriate trace sense (cannot jump across interfaces), the normal trace of the vector-valued field $\mathbf{u} := -\nabla p$ (Darcy velocity, flux or flow field, dual variable) is normal-component continuous in the appropriate trace sense (what flows out from a subdomain through an interface flows in through this interface in the neighboring subdomain), the Darcy velocity is related to the pressure by the constitutive law $\mathbf{u} := -\nabla p$, and the Darcy velocity is in equilibrium with the load $f$ as expressed by $\nabla \cdot \mathbf{u} = f$.

## 2.3 Numerical approximation: spatial discretization, temporal discretization, iterative linearization, and iterative linear algebraic resolution

It is in general not possible to find analytically the exact potential $p$ and the exact flux $\mathbf{u}$ of problems like (2.1) or (2.3). Then, we want to find an approximation obtained by a mathematical algorithm implemented on a computer. Typically, we will take a subdivision of $\Omega$ into computational cells, a spatial mesh $\mathcal{T}_h$ of $\Omega$, with a characteristic size $h$, and use some spatial discretization scheme. For (2.3), we additionally consider a subdivision of the time interval $(0, t_F)$ into time steps, a temporal mesh of $(0, t_F)$, with a characteristic size $\tau$, and proceed to a temporal discretization. When a nonlinear function as $\mathbf{K}(\cdot)$ in (2.4) is present, we typically need an iterative linearization, to, on each iteration step $k$, approximate a nonlinear operator by a linear one. Finally, for steady problems or implicit discretizations of unsteady problems, we need an iterative linear algebraic solver for an approximate solution of a (large sparse) system of linear algebraic equations, with an iteration index $i$. In this text, we focus on discretizations yielding numerical approximations to the exact flux $\mathbf{u}$ (rather than to the exact potential $p$) that we denote by $\mathbf{u}_h^{k,i}$ or $\mathbf{u}_{h\tau}^{k,i}$.

## 2.4 Total error and its components

In principle, the numerical approximation is not exact, so that $\mathbf{u}_h^{k,i} \neq \mathbf{u}$ for (2.1) and $\mathbf{u}_{h\tau}^{k,i} \neq \mathbf{u}$ for (2.3). Evaluating the difference $\mathbf{u} - \mathbf{u}_h^{k,i}$ or $\mathbf{u} - \mathbf{u}_{h\tau}^{k,i}$ in some norm then gives the notion of the total error. This is in our exposition composed of the spatial discretization error, temporal discretization error, iterative linearization error, and iterative linear algebraic solver error.<fnref n="1" />

### 2.4.1 Intrinsic error measure

In the choice of the norm to measure the distances $\mathbf{u} - \mathbf{u}_h^{k,i}$ or $\mathbf{u} - \mathbf{u}_{h\tau}^{k,i}$, we in this contribution stick to problem-induced norms, written as

$$\|\mathbf{u} - \mathbf{u}_h^{k,i}\| \quad \text{or} \quad \|\mathbf{u} - \mathbf{u}_{h\tau}^{k,i}\|. \qquad (2.6)$$

These intrinsic error measures arise from the problem considered, the mathematical function spaces adopted, and are defined so as to yield the most direct link to the sum of the residual and nonconformity in $\mathbf{u}_h^{k,i}$ or $\mathbf{u}_{h\tau}^{k,i}$, cf. Theorem 4.9 below. In the Laplace case (2.2), our choice of the intrinsic error measure in particular gives

$$\|\mathbf{u} - \mathbf{u}_h^{k,i}\| = \|\mathbf{u} - \mathbf{u}_h^{k,i}\|, \qquad (2.7)$$

where $\|\cdot\|$ is the $L^2(\Omega)$ norm on $\Omega$. This $L^2(\Omega)$ flux error $\|\mathbf{u} - \mathbf{u}_h^{k,i}\|$ is a physically crucialy important quantity in the applications.

### 2.4.2 Spatial discretization error

The spatial discretization error is the error arising from the discretization on the spatial mesh $\mathcal{T}_h$ of $\Omega$. For steady problems (2.1) and supposing that we have obtained the approximate flux $\mathbf{u}_h$ as specified “on paper”, i.e., in particular solving “exactly” the associated (non)linear systems (so that there is no

<footnote n="1">¹In this text, we suppose exact computer implementation, i.e., we neglect rounding errors coming from approximate computer work with real numbers via floating point arithmetic. Such errors can also be taken into account, as well as additional errors such as those from (iterative) regularization or model selection.</footnote>
linearization index $k$ and no algebraic solver index $i$), this is the only error component. In the Laplace case (2.2), the spatial discretization error writes

$$ \|u - u_h\|. \qquad (2.8) $$

In the setting of locally conservative methods that we consider, the approximate flux $u_h$ will satisfy

$$ u_h \in H(\operatorname{div}, \Omega) \quad \text{with} \quad \nabla \cdot u_h = f \qquad (2.9) $$

for a piecewise constant source term $f$. Property (2.9) is precisely the “local conservation” from “locally conservative” methods. Additionally, we will be able to find an approximate potential $\tilde{p}_h$ such that

$$ -\nabla \tilde{p}_h|_K = u_h|_K $$

for all elements $K$ of the mesh $T_h$, so that the constitutive law is satisfied. In comparison with the continuous level (2.5), we see that the only deficiency of the discrete level is that the approximate potential $\tilde{p}_h$ is nonconforming, does not belong to the correct space,

$$ \tilde{p}_h \notin H_0^1(\Omega), \qquad (2.10) $$

in contrast to the situation at the continuous level, where $p \in H_0^1(\Omega)$ as per (2.5).

### 2.4.3 Iterative linear algebraic solver error

In order to obtain the approximate solution $u_h$ as specified by some spatial discretization scheme, one typically needs to solve a (large, sparse) system of linear algebraic equations. In practice, we typically proceed iteratively, yielding an approximation $u_h^i$ to $u_h$ on each algebraic solver step $i \ge 0$. Then the iterative linear algebraic solver error is

$$ \|u_h - u_h^i\|. \qquad (2.11) $$

### 2.4.4 Iterative linearization error

In presence of nonlinearities, one typically employs an iterative approximation of a nonlinear operator $A(\cdot)$ by linear operators $A^k(\cdot)$ (on the continuous or on the discrete level). For steady problems, this yields an approximation $u_h^k$ to $u_h$ on each iterative linearization step $k \ge 0$. The corresponding iterative linearization error can then be structurally expressed as

$$ \|u_h - u_h^k\|. \qquad (2.12) $$

### 2.4.5 Temporal discretization error

The temporal discretization error is the error arising from the discretization on the temporal mesh of the time interval $(0, t_F)$. Naturally, it only appears for unsteady problems of the form (2.3), when some temporal discretization has been applied.

### 2.4.6 Error components

Wrapping up the previous sections, the total error is typically composed of several error components. For problems of the form (2.3)–(2.4) namely, we still have the potential and the Darcy velocity

$$ p \in H_0^1(\Omega), \quad u := -\underline{K}(|\nabla p|)\nabla p \in H(\operatorname{div}, \Omega), \quad \partial_t p + \nabla \cdot u = f, \qquad (2.13) $$

as in (2.5). Then, for a numerical approximation $u_{h\tau}^{k,i}$, we will have

$$ \text{total error} \qquad \|u - u_{h\tau}^{k,i}\|, \qquad (2.14a) $$

$$ \text{temporal discretization error} \qquad \|u - u_{\tau}\|, \qquad (2.14b) $$

$$ \text{spatial discretization error} \qquad \|u_{\tau} - u_{h\tau}\|, \qquad (2.14c) $$

$$ \text{iterative linearization error} \qquad \|u_{h\tau} - u_{h\tau}^k\|, \qquad (2.14d) $$

$$ \text{iterative linear algebraic solver error} \qquad \|u_{h\tau}^k - u_{h\tau}^{k,i}\|, \qquad (2.14e) $$

in extension of (2.8), (2.12), and (2.11). Details are developed in Sections 6 and 7.
## 2.5 A posteriori error estimate

An a posteriori error estimate is real number $\eta(\mathbf{u}_h^{k,i})$ or $\eta(\mathbf{u}_{h\tau}^{k,i})$ that bounds the error between the unknown exact flux $\mathbf{u}$ of the partial differential equation like (2.1) or (2.3) and the available numerical approximation $\mathbf{u}_h^{k,i}$ or $\mathbf{u}_{h\tau}^{k,i}$ of $\mathbf{u}$. Importantly, $\eta(\mathbf{u}_h^{k,i})$ or $\eta(\mathbf{u}_{h\tau}^{k,i})$ has to be fully computable from the approximate solution $\mathbf{u}_h^{k,i}$ or $\mathbf{u}_{h\tau}^{k,i}$. The following are some highly desirable properties.

### 2.5.1 Guaranteed error upper bound (reliability)

We say that an a posteriori error estimate gives a *guaranteed error upper bound* when the intrinsic error measure (2.6) is bounded by it as

$$ \|\mathbf{u} - \mathbf{u}_h^{k,i}\| \le \eta(\mathbf{u}_h^{k,i}) \quad \text{or} \quad \|\mathbf{u} - \mathbf{u}_{h\tau}^{k,i}\| \le \eta(\mathbf{u}_{h\tau}^{k,i}). \qquad (2.15) $$

This allows to *certify* the error. This may appear as a miracle, since the left-hand side of (2.15) is unknown, whereas the right-hand side of (2.15) is known, but is indeed possible, as we will see in the developments below. The key will be to follow at the discrete level the physical and mathematical properties of the weak solutions of PDEs as expressed by (2.5) and (2.13). In the setting of Section 2.4.2, namely, from (2.5) and since $\mathbf{u}_h \in \mathbf{H}(\mathrm{div}, \Omega)$ with $\nabla \cdot \mathbf{u}_h = f$, we easily show that

$$ \|\mathbf{u} - \mathbf{u}_h\| \le \eta(\mathbf{u}_h) := \|\mathbf{u}_h + \nabla \zeta_h\| \qquad (2.16) $$

for an arbitrary

$$ \zeta_h \in H_0^1(\Omega); \qquad (2.17) $$

we will see the details in Theorem 4.18 below. Here, $\zeta_h$ is discrete, piecewise polynomial, and obtained by some *local* modification from $\tilde{p}_h$ (postprocessing and averaging); we call it a *potential reconstruction*. Similarly, when $\mathbf{u}_h^{k,i} \notin \mathbf{H}(\mathrm{div}, \Omega)$ or $\nabla \cdot \mathbf{u}_h^{k,i} \neq f$ (namely in presence of linear and nonlinear solvers), an *equilibrated flux reconstruction* $\sigma_h^{k,i}$ such that

$$ \sigma_h^{k,i} \in \mathbf{H}(\mathrm{div}, \Omega) \quad \text{and} \quad \nabla \cdot \sigma_h^{k,i} = f \qquad (2.18) $$

(we will also consider merely $\nabla \cdot \sigma_h^{k,i} = f + \rho_h^{k,i}$ where $\rho_h^{k,i}$ is a remainder made as small as necessary) will be the key together with the potential reconstruction (2.17).

### 2.5.2 Error lower bound (efficiency)

Relation (2.15) is only one-sided, so that it might potentially happen that the left-hand side of (2.15) is zero (small), whereas the right-hand side is nonzero (big). A mathematical equivalence between the error and the a posteriori error estimator is ensured when, in addition to (2.15), there also holds

$$ \eta(\mathbf{u}_h^{k,i}) \le C \|\mathbf{u} - \mathbf{u}_h^{k,i}\| \quad \text{or} \quad \eta(\mathbf{u}_{h\tau}^{k,i}) \le C \|\mathbf{u} - \mathbf{u}_{h\tau}^{k,i}\|, \qquad (2.19) $$

which is called (global) *efficiency*. Here, $C$ is a generic constant. Actually, (2.19) may only hold when on the right-hand side, there are some additional terms without structural importance, i.e., of vanishing importance with decreasing the mesh size $h$ and the time step $\tau$ (so-called “data oscillation” terms).

Strengthening (2.19) to hold locally in space (and in time) is actually often possible, and highly desirable. This then means that when we predict the presence of error in some mesh element or time step, then it is indeed localized there, or in a small neighborhood. In the example (2.16), in particular, we easily have

$$ \|\mathbf{u}_h + \nabla \zeta_h\|^2 = \sum_{K \in \mathcal{T}_h} \|\mathbf{u}_h + \nabla \zeta_h\|_K^2, \qquad (2.20) $$

and it is possible to strengthen (2.19) to

$$ \eta_K(\mathbf{u}_h) \le C \|\mathbf{u} - \mathbf{u}_h\|_{\omega_K}, \qquad (2.21) $$

where $\omega_K$ are the simplices of the mesh $\mathcal{T}_h$ sharing a vertex with the simplex $K$. This is called *local efficiency*.
### 2.5.3 Robustness with respect to the spatial and temporal domains as well as physical and numerical parameters

There may be various ways how to simultaneously achieve (2.15) and (2.19) (and also (2.21)). Then, a crucial indicator of the quality of the derived estimate $\eta(\mathbf{u}_h^{k,i})$ or $\eta(\mathbf{u}_{h\tau}^{k,i})$ is the nature and behavior of the “generic” constant $C$ from (2.19) (and (2.21)) (recall that there is no unknown constant in (2.15)). Ideally, $C$ is truly generic, only depending on the space dimension $d$ and the shape regularity (smallest angle) of the spatial meshes $\mathcal{T}_h$ (say taking a value $2-10$ in practical applications when $1 \le d \le 3$). In any case, it is requested that $C$ is independent of the spatial domain $\Omega$ (namely its size) and the temporal domain $(0, t_F)$ (the final simulation time $t_F$) and of any numerical parameters, namely the mesh size $h$ and the time step $\tau$. It is also desirable to have $C$ independent of the physical parameters such us medium properties (the smallest and highest (eigen)values of diffusion tensor $\mathbf{K}$ or specifications of the nonlinearity in $\mathbf{K}(\cdot)$). We then speak of robustness.

### 2.5.4 Asymptotic exactness

One characterization of the quality of an a posteriori error estimate is the so-called effectivity index given by

$$I_{\text{eff}} := \frac{\eta(\mathbf{u}_h^{k,i})}{\|\mathbf{u} - \mathbf{u}_h^{k,i}\|} \quad \text{or} \quad I_{\text{eff}} := \frac{\eta(\mathbf{u}_{h\tau}^{k,i})}{\|\mathbf{u} - \mathbf{u}_{h\tau}^{k,i}\|} \qquad (2.22)$$

Recall that in our setting, from (2.15), $I_{\text{eff}}$ is greater than or equal to 1. If $I_{\text{eff}}$ goes to the optimal value of 1 with decreasing mesh size $h$ and time step $\tau$, we speak of asymptotic exactness. We will see such cases below.

### 2.5.5 Inexpensive evaluation

There exist very sharp estimates $\eta(\mathbf{u}_h^{k,i})$ or $\eta(\mathbf{u}_{h\tau}^{k,i})$ in (2.15), which, however, are expensive to compute and evaluate, comparatively as expensive as to compute the numerical approximations $\mathbf{u}_h^{k,i}$ or $\mathbf{u}_{h\tau}^{k,i}$ themselves. Indeed, from (2.16) and (2.5), as sharp as requested upper bound can be obtained upon a solution of some global problem with increasing size, see Theorem 4.9 below for more details. Such classes of estimates are, however, excluded from our considerations, as not suitable in practice (they may serve to certify the error but it is hard to think of their use in adaptivity because of their cost). Some other classes of estimates request a solution of some small, local, mutually independent problems in patches of mesh elements. These are termed implicit and can be labelled “inexpensive”. In this contribution, however, we do not consider such estimates either, and we only focus on estimates whose evaluation from the approximate solution $\mathbf{u}_h^{k,i}$ or $\mathbf{u}_{h\tau}^{k,i}$ is yet cheaper, explicit. For example, the potential reconstruction $\zeta_h$ from (2.17) is obtained from $\mathbf{u}_h$ (passing through $\tilde{p}_h$) by some local postprocessing and averaging. On general polytopal meshes $\mathcal{T}_H$, we will then develop estimates whose evaluation merely consists in multiplications of local vectors by local matrices, where the local matrices are pre-processed from the scheme or cell geometry at hand and the local vectors are the degrees of freedom representing the available values of the local fluxes and potentials. These estimate are then as inexpensive as possible.

### 2.5.6 Estimating the error components

We have seen in Section 2.4.6 that the total error in a numerical approximation of, say (2.3)–(2.4), has several components, namely due to spatial discretization, temporal discretization, iterative linearization, and iterative linear algebraic resolution. Congruently, an a posteriori error estimator should identify and estimate these error components, developing (2.15) structurally into

$$\|\mathbf{u} - \mathbf{u}_h^{k,i}\| \le \eta(\mathbf{u}_h^{k,i}) \le \eta_{\text{sp}}(\mathbf{u}_h^{k,i}) + \eta_{\text{lin}}(\mathbf{u}_h^{k,i}) + \eta_{\text{alg}}(\mathbf{u}_h^{k,i}) \qquad (2.23a)$$
or

$$\|\mathbf{u} - \mathbf{u}_{h\tau}^{k,i}\| \le \eta(\mathbf{u}_{h\tau}^{k,i}) \le \eta_{\text{sp}}(\mathbf{u}_{h\tau}^{k,i}) + \eta_{\text{tm}}(\mathbf{u}_{h\tau}^{k,i}) + \eta_{\text{lin}}(\mathbf{u}_{h\tau}^{k,i}) + \eta_{\text{alg}}(\mathbf{u}_{h\tau}^{k,i}), \qquad (2.23b)$$

in congruence with the example of error components in (2.14).

## 2.6 Adaptivity: balancing the error component estimates

Having approximately identified the total error and its components as per (2.23), one can think of how to use this information. We are in particular going to develop adaptivity consisting in balancing of the
estimated error components. For (2.23b), for example, this consists in making all the estimators $\eta_{sp}(\mathbf{u}_{h\tau}^{k,i})$, $\eta_{lm}(\mathbf{u}_{h\tau}^{k,i})$, $\eta_{lin}(\mathbf{u}_{h\tau}^{k,i})$, and $\eta_{alg}(\mathbf{u}_{h\tau}^{k,i})$ of comparable size, uniformly throughout the entire simulation. Indeed, shall for example the iterative linear algebraic solver error estimate $\eta_{alg}(\mathbf{u}_{h\tau}^{k,i})$ be orders of magnitude smaller than the spatial or temporal discretization error estimates, we may be wasting our computational resources (unfortunately, this is common in practice). On a more common note, say for a steady problem, we will develop

$$[\eta_{sp}(\mathbf{u}_h^{k,i})]^2 = \sum_{K \in \mathcal{T}_h} [\eta_{K,sp}(\mathbf{u}_h^{k,i})]^2,$$

in extension of (2.20). We then design algorithms ensuring that the elementwise a posteriori error estimators $\eta_{K,sp}(\mathbf{u}_h^{k,i})$ have all comparable values. This is termed *mesh adaptivity* (adaptive mesh refinement).

## 2.7 Mass balance recovery at any step of numerical resolution

A locally conservative method produces an approximation to the Darcy velocity $\mathbf{u}$ from (2.5) in a way that

$$\mathbf{u}_h \in \mathbf{H}(\text{div}, \Omega) \quad \text{with} \quad \nabla \cdot \mathbf{u}_h = f$$

for a piecewise constant source term $f$, see (2.9). Unfortunately, during an iterative linearization or during iterative linear algebraic resolution, this may be lost in that

$$\mathbf{u}_h^{k,i} \notin \mathbf{H}(\text{div}, \Omega) \quad \text{and/or} \quad \nabla \cdot \mathbf{u}_h^{k,i} \neq f.$$

One of the specificities of our approach is that in the construction of our a posteriori error estimates giving a guaranteed error upper bound and distinguishing the total error components, we actually recover a locally conservative flux field (Darcy velocity) $\sigma_h^{k,i}$ such that

$$\sigma_h^{k,i} \in \mathbf{H}(\text{div}, \Omega) \quad \text{and} \quad \nabla \cdot \sigma_h^{k,i} = f,$$

see (2.18). In practice, we may go for a cheaper variant where $\nabla \cdot \sigma_h^{k,i} = f + \rho_h^{k,i}$ with $\rho_h^{k,i}$ is a remainder made as small as necessary.

## 2.8 Bibliographic resources

Focusing on the spatial discretization error, it has been rather soon understood how to structurally obtain a guaranteed error upper bound as in (2.16): this follows from the Prager-Synge equality [216], relying on the physical and mathematical properties of the weak solutions of PDEs as discussed in Section 2.2. In particular, in finite element methods, where one obtains an approximate solution $p_h \in H_0^1(\Omega)$, one was in quest of an equilibrated flux reconstruction $\sigma_h \in \mathbf{H}(\text{div}, \Omega)$ with $\nabla \cdot \sigma_h = f$ as in (2.18). A local construction, is however, a bit involved in finite elements, so that its establishment was rather labyrinthine, see Ladevèze [174], Ladevèze and Leguillon [175], Bossavit [56], Destuynder and Métivet [97], Larson and Niklasson [177], Luce and Wohlmuth [183], Vohralík [248], Braess and Schöberl [57], Cottereau et al. [96], Ern and Vohralík [124, 125], and the references therein.

In contrast, in methods locally conservative by construction, we already have $\mathbf{u}_h \in \mathbf{H}(\text{div}, \Omega)$ with $\nabla \cdot \mathbf{u}_h = f$: directly from mixed finite elements and by local postprocessing for finite volume-type methods, as detailed below. And since a potential reconstruction $\zeta_h \in H_0^1(\Omega)$ as in (2.17) is rather easy to design, a posteriori error estimates of the form (2.16) were put in place in Achdou et al. [3], Ainsworth [12, 13], Kim [162], and Vohralík [247, 249], see also the references therein. In discontinuous Galerkin and related methods, $\mathbf{u}_h \in \mathbf{H}(\text{div}, \Omega)$ with $\nabla \cdot \mathbf{u}_h = f$ is not directly available but can be readily reconstructed, see Bastian and Rivière [34], Kim [162], and Ern et al. [119]. More references specifically for lowest-order locally conservative methods/finite volumes are discussed in Sections 4.12 and 5.11 below.

Unsteady problems were analyzed in Picasso [214], Verfürth [242], Bergam et al. [47], Ern and Vohralík [123], and Ern et al. [120], see also the references therein. Iterative linear algebraic solver error and iterative linearization error were in particular addressed in Becker et al. [39], Jiránek et al. [154], Arioli et al. [25], Ern and Vohralík [124], Heid and Wihler [147], Gantner et al. [136], Haberl et al. [142], Mitra and Vohralík [190], and the references therein. In-depth presentation of a posteriori error estimates and adaptivity can then be found in the monographs by Synge [236], Hlaváček et al. [152], Zeidler [263], Ainsworth and Oden [14], Babuška and Strouboulis [29], Bangerth and Rannacher [31], Han [144], Ladevèze and Pelle [176], Repin [219], Nochetto et al. [198], Deuflhard and Weiser [99], and Verfürth [245], see also Chamoin and Legoll [76], Bartels and Kaltenbach [33], Mghazli [186, 187], and Smears [232].
# 3 Sobolev spaces and their piecewise polynomial subspaces, meshes, and finite volume methods

In this section, we describe the infinite-dimensional Sobolev spaces $H_0^1(\Omega)$ and $\mathcal{H}(\text{div}, \Omega)$ and their finite-dimensional subspaces formed by piecewise polynomials. We also set up the notation, define a simplicial mesh, and describe the considered polytopal meshes. We finally introduce the basic principle of the finite volume and related lowest-order locally conservative methods.

## 3.1 Basic notation

We introduce here the basic notation used throughout the manuscript.

### 3.1.1 Domain $\Omega$

We let $\Omega \subset \mathbb{R}^d$ be an open bounded connected set in a form of an interval or polygon or polyhedron for respectively $d = 1, 2, 3$ space dimensions, or a polytope in general for $d \ge 1$, with a Lipschitz-continuous boundary $\partial\Omega$. The notation $\omega$ is reserved for a subdomain $\omega \subset \Omega$.

### 3.1.2 Lebesgue spaces

For $\omega \subset \Omega$, we let $L^2(\omega)$ stand for square-integrable scalar-valued functions $v, w : \omega \to \mathbb{R}$ with the scalar product

$$ (v, w)_\omega := \int_\omega v(x)w(x) \, dx. \qquad (3.1a) $$

Square-integrable vector-valued functions $v, w : \omega \to \mathbb{R}^d$ are then collected in $L^2(\omega) := [L^2(\omega)]^d$, with the scalar product

$$ (v, w)_\omega := \int_\omega v(x) \cdot w(x) \, dx. \qquad (3.1b) $$

The associated norms are denoted by $\|\cdot\|_\omega$. When $\omega = \Omega$, we drop the subscripts.

### 3.1.3 Measures, sizes, and cardinalities

We write $|\omega|$ for the Lebesgue measure of a domain $\omega$, $|v|$ for the Euclidean size of a vector $v \in \mathbb{R}^d$, and $|S|$ for the cardinality of a set $S$. We denote by $n_\omega$ the exterior unit normal of a domain $\omega$.

## 3.2 Infinite-dimensional Sobolev spaces $H_0^1(\Omega)$ and $\mathcal{H}(\text{div}, \Omega)$

In the mathematical sense, the solutions of partial differential equations live in specific infinite-dimensional spaces that admit generalized (weak) partial derivatives and divergences. This is in particular the Sobolev space $H_0^1(\Omega)$ for scalar-valued functions and the Sobolev space $\mathcal{H}(\text{div}, \Omega)$ for vector-valued functions. We recapitulate some basics about these spaces here. We follow Allaire [15], Ern and Guermond [118], and Vohralík [251, 252].

### 3.2.1 Sobolev space $H_0^1(\Omega)$

We start by the notion of a weak partial derivative. Let $\mathcal{D}(\Omega)$ be the space of functions from $C^\infty(\Omega)$ with a compact support in $\Omega$.

**Definition 3.1** (Weak partial derivative). Let a scalar-valued function $v : \Omega \to \mathbb{R}$ be given. We say that $v$ admits a weak $i$-th partial derivative, $1 \le i \le d$, if

1. $v \in L^2(\Omega)$;

2. there exists a function $w_i : \Omega \to \mathbb{R}$ such that

(a) $w_i \in L^2(\Omega)$;

(b) $(v, \partial_{x_i}\varphi) = -(w_i, \varphi) \quad \forall \varphi \in \mathcal{D}(\Omega)$.

The function $w_i$ is called the weak $i$-th partial derivative of $v$. We use the notation $\partial_{x_i}v = w_i$.

We then define a weak gradient:
<figure><img src="image_2.png" /><figcaption>Figure 1: Example of a function belonging to the Sobolev space $H_0^1(\Omega)$ (left; the function actually lies in $C^0(\bar{\Omega})$, but not in $C^1(\bar{\Omega})$). Example of a function not belonging to the Sobolev space $H_0^1(\Omega)$ but belonging to the broken Sobolev space $H^1(\mathcal{T}_h)$ for a mesh $\mathcal{T}_h$ composed of two triangles (right; note that there is no trace-continuity on the common face of the two mesh elements)</figcaption></figure>

**Definition 3.2 (Weak gradient).** Let a scalar-valued function $v: \Omega \to \mathbb{R}$ be given. We say that $v$ admits a weak gradient if $v$ admits the weak i-th partial derivative for all $1 \le i \le d$. We set

$$ \nabla v := (\partial_{x_1} v, \dots, \partial_{x_d} v)^t. \qquad (3.2) $$

**Definition 3.3 (The space $H^1(\Omega)$).** The space $H^1(\Omega)$ is the space of all the functions that admit the weak gradient.

Let us recall from [15, 118] that $H^1(\Omega)$ is a Hilbert space for the scalar product given by $(u, v) + (\nabla u, \nabla v)$.

**Definition 3.4 (The space $H_0^1(\Omega)$).** The space $H_0^1(\Omega)$ is the space of functions $v \in H^1(\Omega)$ such that $v|_{\partial\Omega} = 0$.

Recalling again from [15, 118], $H_0^1(\Omega)$ is a Hilbert space for the scalar product given by $(\nabla u, \nabla v)$.

**Remark 3.5 (Trace and trace continuity).** In the above definition, we have used the notation $v|_{\partial\Omega}$. As the functions from $H^1(\Omega)$ are taken from the Lebesgue space $L^2(\Omega)$, $v|_{\partial\Omega}$ in the sense of a restriction is not defined since $\partial\Omega$ is a set of measure zero. In our notation, by $v|_{\partial\Omega}$, we mean the trace of $v \in H^1(\Omega)$ on the boundary of $\Omega$, a key notion for the Sobolev space $H^1(\Omega)$. Crucially, $v|_{\partial\Omega}$ in the sense of traces and $v|_{\partial\Omega}$ as a restriction coincide for any function $v \in H^1(\Omega)$ that is continuous, $v \in H^1(\Omega) \cap C^0(\bar{\Omega})$. Let us also recall that $H^1(\Omega) \subset C^0(\bar{\Omega})$ in one space dimension but $H^1(\Omega) \not\subset C^0(\bar{\Omega})$ in multiple space dimensions. Crucially, functions $v$ in $H^1(\Omega)$, though not necessarily continuous, are continuous in the trace sense. An illustration of a function $v$ belonging to $H^1(\Omega)$ is given in Figure 1, left, and of a function $v$ not belonging to $H^1(\Omega)$ in Figure 1, right. We refer for details to [15, 118, 251].

### 3.2.2 Sobolev space $H(\operatorname{div}, \Omega)$

We now turn to weak a divergence.

**Definition 3.6 (Weak divergence).** Let a vector-valued function $v: \Omega \to \mathbb{R}^d$ be given. We say that $v$ admits a weak divergence if

1. $v \in L^2(\Omega)$;

2. there exists a function $w: \Omega \to \mathbb{R}$ such that

(a) $w \in L^2(\Omega)$;

(b) $(v, \nabla w) = -(w, \varphi) \quad \forall \varphi \in \mathcal{D}(\Omega)$.

The function $w$ is called the weak divergence of $v$. We use the notation $\nabla \cdot v = w$.

**Definition 3.7 (The space $H(\operatorname{div}, \Omega)$).** The space $H(\operatorname{div}, \Omega)$ is the space of all the functions that admit the weak divergence.

Let us recall from [15, 118] that $H(\operatorname{div}, \Omega)$ is a Hilbert space for the scalar product given by $(u, v) + (\nabla u, \nabla v)$. Similarly to Remark 3.5, there holds:
<figure><img src="image_3.png" /><figcaption>Figure 2: Example of a function belonging to the Sobolev space $H(\text{div}, \Omega)$ (left; on the interface $x = 0$ between the two subdomains $\Omega_1 = (-1, 0) \times (0, 1)$ and $\Omega_2 = (0, 1) \times (0, 1)$, $v|_{\Omega_1} \cdot n = v|_{\Omega_2} \cdot n$, with $n = (1, 0)^t$; $v$ is not continuous for each component (the $y$ component of $v$ is discontinuous, as it passes from value 0 in $\Omega_1$ to a nonzero value in $\Omega_2$) but $v$ is normal-trace continuous). Example of a function not belonging to the Sobolev space $H(\text{div}, \Omega)$ but belonging to the broken Sobolev space $H(\text{div}, \mathcal{T}_h)$ (right, there is no normal trace-continuity across the mesh face at $x = 0$).</figcaption></figure>

**Remark 3.8** (Normal trace and normal-trace continuity). On $\partial\Omega$, functions $v \in H(\text{div}, \Omega)$ admit a normal component in appropriate sense, called normal trace and denoted as $(v \cdot n)|_{\partial\Omega}$. Namely, if $v \in H(\text{div}, \Omega) \cap [C^0(\bar{\Omega})^d]$, then $v \cdot n$ is the usual normal component of $v$ on the boundary $\partial\Omega$. Functions $v$ in $H(\text{div}, \Omega)$ are not necessarily continuous in each component, not necessarily normal-component continuous, but are normal-trace continuous. An illustration of a function $v$ belonging to $H(\text{div}, \Omega)$ is given in Figure 2, left, and of a function $v$ not belonging to $H(\text{div}, \Omega)$ in Figure 2, right. We refer for details to [15, 118].

The Sobolev spaces $H_0^1(\Omega)$ from Definition 3.4 and $H(\text{div}, \Omega)$ from Definition 3.7 form exactly the right setting for the following Green theorem:

**Theorem 3.9** (Green theorem). Let $v \in H_0^1(\Omega)$ and let $w \in H(\text{div}, \Omega)$. Then

$$ (w, \nabla v) + (\nabla \cdot w, v) = 0. \qquad (3.3) $$

In the following, we will also employ the above notation and developments for subdomains $\omega \subset \Omega$.

## 3.3 Simplicial meshes

In order to form finite-dimensional subspaces of $H_0^1(\Omega)$ and $H(\text{div}, \Omega)$, we will need a simplicial mesh $\mathcal{T}_h$ of $\Omega$. This a finite collection of closed $d$-simplices (intervals for $d = 1$, triangles for $d = 2$, tetrahedra for $d = 3$) such that $\cup_{K \in \mathcal{T}_h} K = \bar{\Omega}$ and such that the intersection of two different simplices is either empty or their entire $l$-dimensional face, $0 \le l \le (d-1)$. Here, for example for space dimension $d = 3$, a vertex is a 0-dimensional face, an edge is a 1-dimensional face, and a face is a 2-dimensional face. Henceforth, by face, we mean $(d-1)$-dimensional face. Figure 3 gives illustrations for $d = 2$ and $d = 3$. We denote by

$$ \theta_{\mathcal{T}_h} := \max_{K \in \mathcal{T}_h} \frac{h_K}{\iota_K} \qquad (3.4) $$

the shape-regularity parameter of the mesh $\mathcal{T}_h$, where $h_K$ is the diameter of the simplex $K \in \mathcal{T}_h$ and $\iota_K$ is the diameter of largest ball inscribed in $K$. This parameter expresses how much the mesh is “distorted”: it is close to 1 for meshes of equilateral simplices and increases in presence of small angles. For each mesh $\mathcal{T}_h$, $\theta_{\mathcal{T}_h}$ is a bounded number. Our theoretical constants such as $C$ from (2.19) and (2.21) will depend on $\theta_{\mathcal{T}_h}$. If we consider a sequence of meshes, then we need $\theta_{\mathcal{T}_h}$ uniformly bounded.

The set of the (($d-1$)-dimensional) faces of $\mathcal{T}_h$ is denoted as $\mathcal{F}_h$. It is divided into interior faces $\mathcal{F}_h^{\text{int}}$ and boundary faces $\mathcal{F}_h^{\text{ext}}$: $\sigma \in \mathcal{F}_h^{\text{int}}$ if there exist $K, L \in \mathcal{T}_h$, $K \neq L$, such that $\sigma = K \cap L$, and $\sigma \in \mathcal{F}_h^{\text{ext}}$ if there exist $K \in \mathcal{T}_h$ such that $\sigma \subset K \cap \partial\Omega$. We will also need the notation $\mathcal{V}_h$ for the vertices of the mesh $\mathcal{T}_h$ and $\mathcal{V}_\sigma$ for the vertices of the face $\sigma \in \mathcal{F}_h$. For an interior face $\sigma \in \mathcal{F}_h^{\text{int}}$, we fix an arbitrary orientation and denote the corresponding unit normal vector by $n_\sigma$. For a boundary face $\sigma \in \mathcal{F}_h^{\text{ext}}$, we make $n_\sigma$ coincide with the exterior unit normal $n_\Omega$ of $\Omega$. Finally, for all $K \in \mathcal{T}_h$, we denote by $n_K$ the unit normal vector to $\sigma$ pointing out of $K$ and by $\mathcal{F}_K$ the set of all faces $\sigma$ of $K$.
<figure><img src="image_4.png" /><figcaption>Figure 3: Simplicial mesh $\mathcal{T}_h$ for $d = 2$ (left) and $d = 3$ (right)</figcaption></figure>

## 3.4 Broken Sobolev spaces $H^1(\mathcal{T}_h)$ and $H(\operatorname{div}, \mathcal{T}_h)$

Before we pass to finite-dimensional subspaces of $H_0^1(\Omega)$ and $H(\operatorname{div}, \Omega)$, it is useful to define the broken Sobolev space $H^1(\mathcal{T}_h)$ as

$$H^1(\mathcal{T}_h) := \{v \in L^2(\Omega); v|_{\mathcal{T}_h} \in H^1(K) \quad \forall K \in \mathcal{T}_h\}, \quad (3.5)$$

and the broken Sobolev space $H(\operatorname{div}, \mathcal{T}_h)$ as

$$H(\operatorname{div}, \mathcal{T}_h) := \{v \in L^2(\Omega); v|_{\mathcal{T}_h} \in H(\operatorname{div}, K) \quad \forall K \in \mathcal{T}_h\}. \quad (3.6)$$

An illustration of a function $v \in H^1(\mathcal{T}_h)$ is given in Figure 1, right, and an illustration of a function $v \in H(\operatorname{div}, \mathcal{T}_h)$ is given in Figure 2, right. In this setting, Remarks 3.5 and 3.8 can respectively be developed into, see [100, Lemma 1.23] or [118, Theorem 18.8]:

**Lemma 3.10** (Trace continuity over mesh faces). There holds

$$v \in H^1(\mathcal{T}_h) \text{ and } (v|_{\mathcal{T}_h})|_{\sigma} = (v|_L)|_{\sigma} \quad \forall \sigma = K \cap L \in \mathcal{F}_h^{\text{int}} \iff v \in H^1(\Omega), \quad (3.7a)$$

$$v \in H^1(\mathcal{T}_h) \text{ and } (v|_{\mathcal{T}_h})|_{\sigma} = (v|_L)|_{\sigma} \quad \forall \sigma = K \cap L \in \mathcal{F}_h^{\text{int}}$$

$$\text{and } (v|_{\mathcal{T}_h})|_{\sigma} = 0 \quad \forall \sigma = K \cap \partial\Omega \in \mathcal{F}_h^{\text{ext}} \iff v \in H_0^1(\Omega). \quad (3.7b)$$

In words, functions from the broken Sobolev space $H^1(\mathcal{T}_h)$ belong to the Sobolev space $H^1(\Omega)$ if and only if their traces are uniquely defined on all mesh interior faces, with the same trace from both elements sharing the given interior face. We say that the functions from $H^1(\Omega)$ are trace continuous. Functions from $H_0^1(\Omega)$ then additionally have the traces zero on boundary faces.

Similarly, following [100, Lemma 1.24] or [118, Theorem 18.10], we have:

**Lemma 3.11** (Normal trace continuity over mesh faces). Functions from $H(\operatorname{div}, \mathcal{T}_h)$ belong to $H(\operatorname{div}, \Omega)$ if and only if their normal traces are continuous over mesh interior faces in appropriate normal-trace sense.

## 3.5 Finite-dimensional subspaces of $H_0^1(\Omega)$ and $H(\operatorname{div}, \Omega)$

We now introduce standard finite-dimensional subspaces of $H_0^1(\Omega)$ and $H(\operatorname{div}, \Omega)$, again following Al- laire [15] and Ern and Guermond [118].

### 3.5.1 Piecewise polynomial subspaces of $L^2(\Omega)$

Let $K \in \mathcal{T}_h$ be a simplex from the mesh $\mathcal{T}_h$ and let $k \ge 0$ be an integer. We denote by $\mathcal{P}_k(K)$ the space of scalar-valued polynomials on $K$ of total degree at most $k$. In particular, $\mathcal{P}_0(K)$ stands for constants on $K$, $\mathcal{P}_1(K)$ for affine functions on $K$, and $\mathcal{P}_2(K)$ for quadratic functions on $K$. The spaces for $k = 1, 2$ and $d = 2, 3$ are schematically visualized in Figure 4. In particular, fixing the point values as per Figure 4 uniquely determines a function in $\mathcal{P}_k(K)$. We will often use the space of piecewise $k$-degree polynomials on the simplicial mesh $\mathcal{T}_h$,

$$\mathcal{P}_k(\mathcal{T}_h) := \{v_h \in L^2(\Omega); v_h|_K \in \mathcal{P}_k(K) \quad \forall K \in \mathcal{T}_h\}. \quad (3.8)$$
<figure><img src="image_5.png" /><figcaption>Figure 4: Spaces $\mathcal{P}_k(\mathcal{T}_h)$ and point values uniquely fixing a function $v_h \in \mathcal{P}_k(\mathcal{T}_h)$ (the so-called Lagrange nodes) for a mesh element $K \in \mathcal{T}_h$. Polynomial degrees $k = 1$ and $k = 2$, space dimensions $d = 2$ (left) and $d = 3$ (right)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Spaces $\mathcal{RT}_0(K)$ and normal fluxes uniquely fixing a function $v_h \in \mathcal{RT}_0(K)$ for a simplex $K \in \mathcal{T}_h$. Space dimensions $d = 2$ (left) and $d = 3$ (right)</figcaption></figure>

### 3.5.2 Lagrange piecewise polynomial subspaces of $H_0^1(\Omega)$

The space $\mathcal{P}_k(\mathcal{T}_h)$ from (3.8) is a finite-dimensional subspace of the space $H^1(\mathcal{T}_h)$ from (3.5) but not of the Sobolev space $H_0^1(\Omega)$. To create a subspace of $H_0^1(\Omega)$ of the so-called Lagrange piecewise polynomials, we write, in several equivalent ways,

$$ 
\begin{align} 
\mathcal{P}_k(\mathcal{T}_h) \cap H_0^1(\Omega) &= \{v_h \in H_0^1(\Omega); v_h|_K \in \mathcal{P}_k(K) \quad \forall K \in \mathcal{T}_h\} \nonumber \\
&= \{v_h \in \mathcal{P}_k(\mathcal{T}_h); (v_h|_K)|_\sigma = (v_h|_L)|_\sigma \quad \forall \sigma = K \cap L \in \mathcal{F}_h^{\text{int}} \nonumber \\
&\quad \text{and } (v_h|_K)|_\sigma = 0 \quad \forall \sigma \subset K \cap \partial\Omega \in \mathcal{F}_h^{\text{ext}}\} \tag{3.9} \\
&= \{v_h \in \mathcal{P}_k(\mathcal{T}_h) \cap C^0(\bar{\Omega}); v_h = 0 \text{ on } \partial\Omega\}. \nonumber
\end{align} 
$$

Indeed, Lemma 3.10 holds with $H^1(\mathcal{T}_h)$ replaced by $\mathcal{P}_k(\mathcal{T}_h)$ and $H_0^1(\Omega)$ replaced by $\mathcal{P}_k(\mathcal{T}_h) \cap H_0^1(\Omega)$, where the traces on faces $\sigma$ now equal the usual restrictions to $\sigma$. Congruently, fixing the point values in each simplex $K$ of the mesh $\mathcal{T}_h$ as schematically visualized in Figure 4 by a unique bullet value yields the $H_0^1(\Omega)$- or $C^0(\bar{\Omega})$-conformity for piecewise $k$-degree polynomials. We again refer for details to [15, 118].

### 3.5.3 Raviart-Thomas piecewise polynomial subspaces of $H(\text{div}, \Omega)$

Recall that $\mathcal{P}_0(K)$ stands for constants on the simplex $K \in \mathcal{T}_h$. Following Whitney [258], Raviart and Thomas [218], and Nédélec [195], see alternatively Boffi et al. [53] or Ern and Guermond [118], we define the space of vector-valued polynomials

$$ 
\begin{align} 
\mathcal{RT}_0(K) &:= [\mathcal{P}_0(K)]^d + \mathbf{x}\mathcal{P}_0(K) \nonumber \\
&= \{v_h(\mathbf{x}) + w_h(\mathbf{x})\mathbf{x}; v_h \in [\mathcal{P}_0(K)]^d, w_h \in \mathcal{P}_0(K)\} \nonumber \\
&= \left\{ \begin{pmatrix} a \\ b \end{pmatrix} + c \begin{pmatrix} x \\ y \end{pmatrix} \quad a, b, c \in \mathbb{R}, \mathbf{x} = (x, y)^{\text{t}} \right\} \quad \text{for } d=2, \tag{3.10} \\
&= \left\{ \begin{pmatrix} a \\ b \\ c \end{pmatrix} + d \begin{pmatrix} x \\ y \\ z \end{pmatrix} \quad a, b, c, d \in \mathbb{R}, \mathbf{x} = (x, y, z)^{\text{t}} \right\} \quad \text{for } d=3. \nonumber
\end{align} 
$$
<figure><img src="image_7.png" /><figcaption>Figure 6: A basis function $v_\sigma$ of $RT_0(K)$ (supported on the simplex $K$) with nonzero normal flux only through the face $\sigma \in \mathcal{F}_K$ (left). A basis function $v_\sigma$ of $RT_0(\mathcal{T}_h) \cap \mathcal{H}(\text{div}, \Omega)$ (supported on all (two here) simplices sharing the face $\sigma$) with nonzero normal flux only through the face $\sigma \in \mathcal{F}_h$ (right).</figcaption></figure>

This space is slightly bigger than just constant vector-valued polynomials $[\mathcal{P}_0(K)]^d$ but slightly smaller than affine vector-valued polynomials $[\mathcal{P}_1(K)]^d$. Its dimension equals $d+1$, i.e., the number of faces of the simplex $K$, which means 3 in two space dimensions and 4 in three space dimensions. Its distinctive properties are

$$ \nabla \cdot v_h \in \mathcal{P}_0(K) \quad \text{and} \quad (v_h \cdot n_K)|_{\sigma} \in \mathcal{P}_0(\sigma) \quad \forall \sigma \in \mathcal{F}_K \quad \text{for } v_h \in \mathcal{RT}_0(K), \qquad (3.11) $$

i.e., the divergence is constant and the normal trace is constant on each face of the simplex $K$ for any function $v_h \in \mathcal{RT}_0(K)$. Let us recall here that for a polynomial $v_h$, the normal trace on a face $\sigma$ equals the restriction to $\sigma$ of the normal component $v_h \cdot n_K$. Moreover, a function $v_h$ in $\mathcal{RT}_0(K)$ is uniquely defined by prescribing the face normal fluxes, i.e.,

$$ \langle v_h \cdot n_K, 1 \rangle_{\sigma} := \int_{\sigma} v_h \cdot n_K \quad \forall \sigma \in \mathcal{F}_K, \qquad (3.12) $$

which is a part of the schematic visualization of the spaces $RT_0(K)$ for $d=2,3$ in Figure 5. We choose the $d+1$ basis functions $v_\sigma$ as associated with the faces $\sigma \in \mathcal{F}_K$. Setting $\langle v_\sigma \cdot n_K, 1 \rangle_\sigma = 1$ and $\langle v_\sigma \cdot n_K, 1 \rangle_{\sigma'} = 0$ for $\sigma' \neq \sigma$, we obtain the explicit expression

$$ v_{\sigma}(x) = \frac{1}{d|K|} (x - a_{K,\sigma}), \quad x \in K, \qquad (3.13) $$

cf. the illustration in Figure 6, left, where $a_{K,\sigma}$ is the vertex of $K$ opposite to $\sigma$. On the simplicial mesh $\mathcal{T}_h$, we now define

$$ \mathcal{RT}_0(\mathcal{T}_h) := \{v_h \in \mathcal{L}^2(\Omega); v_h|_{K} \in \mathcal{RT}_0(K) \quad \forall K \in \mathcal{T}_h\}. \qquad (3.14) $$

This is a finite-dimensional subspace of the space $\mathcal{H}(\text{div}, \mathcal{T}_h)$ from (3.6) but not of the Sobolev space $\mathcal{H}(\text{div}, \Omega)$. To create a subspace of $\mathcal{H}(\text{div}, \Omega)$, we set

$$ \begin{aligned} \mathcal{RT}_0(\mathcal{T}_h) \cap \mathcal{H}(\text{div}, \Omega) &= \{v_h \in \mathcal{H}(\text{div}, \Omega); v_h|_{K} \in \mathcal{RT}_0(K) \quad \forall K \in \mathcal{T}_h\} \\ &= \{v_h \in \mathcal{RT}_0(\mathcal{T}_h); (v_h|_{K} \cdot n_K)|_{\sigma} = -(v_h|_{L} \cdot n_L)|_{\sigma} \quad \forall \sigma = K \cap L \in \mathcal{F}_h^{\text{int}}\}; \end{aligned} \qquad (3.15) $$

it follows by Lemma 3.11 that the above characterizations coincide. Congruently, fixing the face normal values/integrals in each simplex $K \in \mathcal{T}_h$ as schematically visualized in Figure 5 by the arrows in particular yields the $\mathcal{H}(\text{div}, \Omega)$-conformity for piecewise $\mathcal{RT}_0(K)$-polynomials. The space (3.15) of normal-component-continuous piecewise vector-valued polynomials is called the Raviart-Thomas(-Nédélec) space. Its dimension is the number of mesh faces in the mesh $\mathcal{T}_h$. There is one basis function $v_\sigma$ for each face $\sigma \in \mathcal{F}_h$, prescribed such that $\langle v_\sigma \cdot n_\sigma, 1 \rangle_\sigma = 1$ and $\langle v_\sigma \cdot n_{\sigma'}, 1 \rangle_{\sigma'} = 0$ for all $\sigma' \neq \sigma$. An illustration is given in Figure 6, right. We again refer for details to [15, 118].

**Remark 3.12** (Meshes consisting of rectangular parallelepipeds). Above, we only consider simplicial meshes since these are our key to treat general polytopal meshes below. In practice, meshes consisting
<figure><img src="image_8.png" /><figcaption>Figure 7: Polygonal element $K \in \mathcal{T}_H$ (left). Virtual simplicial submesh $\mathcal{T}_K$ of $K$ (middle and right). One point inside $K \in \mathcal{T}_H$ shared by all the simplices in $\mathcal{T}_K$ and faces of $\mathcal{T}_h$ do not subdividing the faces of $\mathcal{T}_H$ (middle). General situation (right).</figcaption></figure>

of rectangles, blocks, or in general rectangular parallelepipeds are often used. On rectangles or blocks, in place of (3.8), piecewise polynomial subspaces of $L^2(\Omega)$ are typically chosen as

$$Q_k(\mathcal{T}_h) := \{v_h \in L^2(\Omega); v_h|_K \in Q_k(K) \quad \forall K \in \mathcal{T}_h\}, \quad (3.16)$$

where $Q_k(K)$ is the space of scalar-valued polynomials on $K$ of degree at most $k$ in each space coordinate. Note that $P_k(\mathcal{T}_h)$ and $Q_k(\mathcal{T}_h)$ coincide for the lowest polynomial degree $k = 0$, being simply piecewise constants with respect to the mesh $\mathcal{T}_h$. The Raviart–Thomas spaces $RT_0(K)$ then take the form

$$RT_0(K) := \left\{ \begin{pmatrix} a+cx \\ b+dy \end{pmatrix} \quad a,b,c,d \in \mathbb{R}, \mathbf{x} = (x,y)^t \right\} \quad \text{for } d=2, \quad (3.17a)$$

and

$$RT_0(K) := \left\{ \begin{pmatrix} a+dx \\ b+ey \\ c+ fz \end{pmatrix} \quad a,b,c,d,e,f \in \mathbb{R}, \mathbf{x} = (x,y,z)^t \right\} \quad \text{for } d=3, \quad (3.17b)$$

in place of (3.10). Subspaces of $H(\text{div}, \Omega)$ are then still created by (3.14)–(3.15).

In the following, we will sometimes employ the above notation and developments for subdomains $\omega \subset \Omega$ and the corresponding submeshes $\mathcal{T}_\omega$ of $\omega$.

## 3.6 Polytopal meshes with virtual simplicial submeshes

Recall that $\Omega \subset \mathbb{R}^d$, $d \ge 1$, is an open bounded connected polytope with a Lipschitz-continuous boundary $\partial \Omega$. Meshes with mesh elements formed by polygons (in two space dimensions), polyhedrons (in three space dimensions), or in general polytopes (in $d$ space dimensions) are useful in practice. An example of a general polygonal element ($d=2$) is given in Figure 7 (left); some examples of polygonal meshes can be found in Figures 8 (left) and 32 and of polyhedral meshes in Figures 34 and 36. These meshes stand in contrast to regular meshes of rectangular parallelepipeds by their generality and, in comparison with meshes formed by simplicial elements, they may contain (much) fewer elements, cf. Figure 8 (right).

### 3.6.1 Polytopal mesh with a virtual simplicial submesh

We call a polytopal mesh a partition $\mathcal{T}_H$ of the domain $\Omega$ into a finite collection of closed $d$-dimensional polytopes $K$ homotopic to a ball and with Lipschitz-continuous boundary $\partial K$ such that $\cup_{K \in \mathcal{T}_H} K = \bar{\Omega}$. We suppose that the intersection of two different polytopes is either empty or their entire $l$-dimensional face, $0 \le l \le (d-1)$, or a collection of their entire $l$-dimensional faces, $0 \le l \le (d-1)$. The elements $K$ can in particular be nonconvex and non star-shaped. The structural assumption we impose is that there exists a virtual, not to be constructed in practice, simplicial mesh $\mathcal{T}_h$ of $\Omega$ in the sense of Section 3.3 such that $\mathcal{T}_K := \mathcal{T}_h|_K$ is a simplicial mesh of the polytope $K$. An illustration is provided in Figures 7 and 8. Altogether, each element $\kappa \in \mathcal{T}_h$ is a $d$-simplex, $\cup_{\kappa \in \mathcal{T}_h} \kappa = \bar{\Omega}$, $\cup_{\kappa \in \mathcal{T}_K} \kappa = K$, and $\cup_{K \in \mathcal{T}_H} \mathcal{T}_K = \mathcal{T}_h$.

### 3.6.2 Mesh faces and mesh vertices

Let $\mathcal{F}_H$ be the set of the $(d-1)$-dimensional faces of $\mathcal{T}_H$. We divide it into interior faces $\mathcal{F}_H^{\text{int}}$ and boundary faces $\mathcal{F}_H^{\text{ext}}$: $\sigma \in \mathcal{F}_H^{\text{int}}$ if there exist $K, L \in \mathcal{T}_H$, $K \ne L$, such that $\sigma \subset K \cap L$, and $\sigma \in \mathcal{F}_H^{\text{ext}}$ if there exist
<figure><img src="image_9.png" /><figcaption>Figure 8: A polygonal mesh $\mathcal{T}_H$ and the corresponding triangular submesh $\mathcal{T}_h$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Example of a polygonal element $K$ with its faces $\mathcal{F}_K$, corresponding face fluxes $U_K^{\text{ext}}$, and pressure head $(P)_K$ (left); virtual simplicial submesh $\mathcal{T}_K$ of $K$ (right)</figcaption></figure>

$K \in \mathcal{T}_H$ such that $\sigma \subset K \cap \partial\Omega$. Recall from Section 3.3 that $\mathcal{F}_h$ is the set of the $(d-1)$-dimensional faces of the simplicial mesh $\mathcal{T}_h$. We denote by $\mathcal{F}_{H,h}$ the set of such faces from $\mathcal{F}_h$ that lie in some polytopal face of $\mathcal{F}_H$. Let $\mathcal{F}_K \subset \mathcal{F}_H$ be the set of the $(d-1)$-dimensional faces of the polytope $K \in \mathcal{T}_H$, and let $\mathcal{F}_{K,h}^{\text{ext}}$ collect those faces of $\mathcal{F}_{H,h}$ that lie on the boundary of the element $K \in \mathcal{T}_H$, whereas $\mathcal{F}_{K,h}^{\text{int}}$ those faces of $\mathcal{F}_h$ that lie inside the element $K \in \mathcal{T}_H$. We also set $\mathcal{F}_{K,h} := \mathcal{F}_{K,h}^{\text{ext}} \cup \mathcal{F}_{K,h}^{\text{int}}$ see Figure 9. If the faces of $\mathcal{T}_h$ do not subdivide the faces of $\mathcal{T}_H$, then $\mathcal{F}_{H,h} = \mathcal{F}_H$ and also $\mathcal{F}_{K,h}^{\text{ext}} = \mathcal{F}_K$. For an interior face $\sigma \in \mathcal{F}_H^{\text{int}}$, we fix an arbitrary orientation and denote the corresponding unit normal vector by $\mathbf{n}_\sigma$. For a boundary face $\sigma \in \mathcal{F}_H^{\text{ext}}$, $\mathbf{n}_\sigma$ coincides with the exterior unit normal $\mathbf{n}_\Omega$ of $\Omega$. Finally, for all $K \in \mathcal{T}_H$, we denote by $\mathbf{n}_K$ the unit normal vector to $\sigma$ pointing out of $K$. To complete the notation, we let $\mathcal{V}_{K,h}$ to be the set of vertices of all the elements $\kappa \in \mathcal{T}_K$, and, as for simplicial meshes, $\mathcal{V}_\sigma$ the set of vertices of the given face $\sigma$.

## 3.7 Basic principle of finite volume and related lowest-order locally conservative methods

We now introduce the basic principle of the finite volume and related lowest-order locally conservative methods considered in this contribution. In Sections 4–7 below, for each considered partial differential equation, a general lowest-order locally conservative discretization is introduced in detail.

Consider a steady diffusion partial differential equation of the form (2.1)–(2.2) and recall (2.5). A general lowest-order locally conservative method looks for real values $U_{K,\sigma}$ which for each polytopal mesh element $K \in \mathcal{T}_H$ and each of its faces $\sigma \in \mathcal{F}_K$ approximate the normal outflux $\langle \mathbf{u} \cdot \mathbf{n}_K, 1 \rangle_\sigma = \int_\sigma \mathbf{u} \cdot \mathbf{n}_K$ from $K$ through the face $\sigma$ such that

$$ \sum_{\sigma \in \mathcal{F}_K} U_{K,\sigma} = (f, 1)_K \quad \forall K \in \mathcal{T}_H. \qquad (3.18a) $$

Here the sum runs through all the faces of the polytope $K$ and, recalling the notation (3.1a), $(f, 1)_K = \int_K f(\mathbf{x}) \, d\mathbf{x}$ is the integral of the source term. From a neighboring element $L \in \mathcal{T}_H$ sharing the same face
$\sigma$, we then impose,

$$U_{L,\sigma} = -U_{K,\sigma}. \qquad (3.18b)$$

Requirement (3.18b) ensures the normal-trace continuity, as in the Raviart-Thomas space $RT_0(\mathcal{T}_h) \cap H(\text{div}, \Omega)$ of (3.15). Then (3.18b) expresses the conservation of mass across mesh faces and (3.18a) the equilibrium with the load, respectively the discrete versions of $u \in H(\text{div}, \Omega)$ and $\nabla \cdot u = f$. Mimicking the remaining items from (2.5), namely $p \in H_0^1(\Omega)$ and $u = -\nabla p$, then has to be built in the construction of the discrete normal fluxes $U_{K,\sigma}$ from elementwise approximations of the primal variable $p_K$.

For an unsteady problem of the form (2.3)-(2.4), we consider an implicit Euler time discretization on discrete times $0 = t^0 < \dots < t^n < \dots t^N = t_F$ with time steps $\tau^n := t^n - t^{n-1}$. We then get, for each discrete time $t^n$, $1 \le n \le N$,

$$\frac{p_K^n - p_K^{n-1}}{\tau^n} |K| + \sum_{\sigma \in \mathcal{F}_K} U_{K,\sigma}^n = (f^n, 1)_K \quad \forall K \in \mathcal{T}_H \qquad (3.19a)$$

and

$$U_{L,\sigma}^n = -U_{K,\sigma}^n, \qquad (3.19b)$$

where $f^n := \int_{t^{n-1}}^{t^n} f dt/\tau^n$. Again, the heart of each method is the construction of the discrete normal fluxes $U_{K,\sigma}^n$ for each polytopal mesh element $K \in \mathcal{T}_H$ and each face $\sigma \in \mathcal{F}_K$ from the elementwise approximations of the primal variable $p_K^n$.

## 3.8 Bibliographic resources

Details on the material rapidly exposed in Sections 3.2–3.6 can be found in Adams [5], Thomas [237], Cia- rlet [92], Allaire [15], Brenner and Scott [61], Boffi et al. [53], Ern and Guermond [118], and Vohralík [251, 252].

Finite volume schemes like (3.18) and (3.19) are very natural and used traditionally in the engineer- ing practice. Amongst first mathematical analyses for diffusion problems, one can cite Manteuffel and White [184], Rose [225], Morton and Süli [191], Faille [132], Herbin [149], and the overview in Eymard et al. [126]. For the closely related mixed finite elements, the construction of the spaces in the lowest-order case can be traced back to Whitney [258], whereas higher-order extensions and analysis have been under- taken in Raviart and Thomas [218] and Nédélec [195], see Brezzi and Fortin [63] for an early overview. The ample connections between finite volumes and mixed finite elements are discussed in, e.g., Russell and Wheeler [227], Agouzal et al. [7], Baranger et al. [32], Younès et al. [260, 259], Chavent et al. [79], Klausen and Russell [164], and Vohralík [246].

Numerous extensions of the initial approaches together with generalizations to polytopal meshes, including a priori analysis (existence, uniqueness, convergence, a priori error estimates), have been un- dertaken recently. Let us cite multi-point finite volumes by Aavatsmark et al. [1, 2], Edwards [113], and Breil and Maire [58], mimetic finite differences by Brezzi et al. [64, 65], mixed and hybrid finite volumes by Eymard et al. [108, 128], multipoint flux mixed finite elements by Wheeler and Yotov et al. [257, 256], compatible discrete operator schemes by Bonelle and Ern [55], mixed virtual elements by Brezzi et al. [62] and Beirão da Veiga et al. [40], hybrid high-order methods by Di Pietro, Ern, and Lemaire [102, 101], or mixed finite elements in [172, 173, 230, 88, 254], see also [52, 150, 87, 107]; unifying frameworks can be found in Droniou et al. [110, 111], Vohralík and Wohlmuth [254], Cockburn et al. [95], Boffi and Di Pietro [54], and the references therein. Amongst recent books, let us cite Jovanović and Süli [155] for finite differences, Beirão da Veiga et al. [41] for mimetic finite differences, Droniou et al. [109] for a unifying approach, and Nordbotten and Keilegavlen [199] for multi-point finite volumes.

# 4 The Poisson equation. Simplicial meshes and spatial discretization error

In this section, we start with the simplest partial differential equation for diffusion processes: the Poisson equation. We first recall the definition of a weak solution based on the Sobolev space $H_0^1(\Omega)$ of Section 3.2.1 and show that the associated flux lies in the Sobolev space $H(\text{div}, \Omega)$ of Section 3.2.2. This is very instrumental and will be repeated in all generalizations of Sections 5–7 below. We then define the energy error and explore its structure for a completely generic approximate solution. We subsequently consider a basic cell-centered finite volume discretization incarnating the principles of Section 3.7. We
only employ in this section a simplicial mesh as per Section 3.3 and define the approximate solution by lifting the face fluxes into a Raviart-Thomas piecewise polynomial space of Section 3.5.3. Next, we introduce and discuss two central notions of this work: the flux and potential reconstructions. These in particular lead to a posteriori error estimates fulfilling all the building principles of Section 2. Illustrative numerical experiments close this section.

## 4.1 The Poisson equation

In this section, for the domain $\Omega$ as specified in Section 3.1.1, we consider the problem of finding $p: \Omega \to \mathbb{R}$ such that

$$ 
\begin{array}{ll}
-\nabla \cdot (\nabla p) = f & \text{in } \Omega, \\
p = 0 & \text{on } \partial\Omega.
\end{array}
\qquad (4.1a) \\
(4.1b)
$$

Here $f \in L^2(\Omega)$ is a source term.

## 4.2 Weak solution, potential, and flux

The Sobolev space $H_0^1(\Omega)$ from Definition 3.4 gives a proper mathematical setting to define a unique solution $p$ of (4.1). Indeed, there is one and only one $p \in H_0^1(\Omega)$ such that

$$ (\nabla p, \nabla v) = (f, v) \quad \forall v \in H_0^1(\Omega); \qquad (4.2) $$

note that this follows trivially from the Riesz representation theorem, since $H_0^1(\Omega)$ is a Hilbert space for the scalar product given by $(\nabla p, \nabla v)$ and since $(f, v)$ is a bounded linear form over $H_0^1(\Omega)$. We call (4.2) the weak formulation of (4.1). From the weak solution $p$, we can define the dual variable

$$ u := -\nabla p. \qquad (4.3) $$

In relation to the applications below, $p$ may be a pressure and $u$ a Darcy velocity, but we rather generically refer to $p$ as potential and to $u$ as flux.

## 4.3 Properties of the exact potential and flux

Using Definitions 3.4 and 3.7 of the spaces $H_0^1(\Omega)$ and $H(\text{div}, \Omega)$, we have the following important result, detailing what we have discovered in (2.5):

**Proposition 4.1** (Properties of the weak solution (4.2)). Let $p$ be the solution of (4.2). Let $u$ be given by (4.3). Then

$$ p \in H_0^1(\Omega), \quad u \in H(\text{div}, \Omega), \quad \nabla \cdot u = f. \qquad (4.4) $$

*Proof.* The weak solution $p$ belongs to $H_0^1(\Omega)$ by definition. In order to verify that $u \in H(\text{div}, \Omega)$ and $\nabla \cdot u = f$, we need to check the three conditions of Definition 3.6. Condition 1 is obvious, as $p \in H_0^1(\Omega)$ and thus $-\nabla p = u$ is square-integrable in view of Definition 3.1, condition 2a. For the function $w$ of condition 2a in Definition 3.6, choose $w := f$ and note that $f \in L^2(\Omega)$ by assumption. Then condition 2b follows immediately from (4.2) and the fact that $D(\Omega) \subset H_0^1(\Omega)$. $\square$

**Remark 4.2** (Mathematical properties of the weak solution). It is remarkable that the flux $u$ prescribed from the potential $p$ by (4.3) belongs to $H(\text{div}, \Omega)$ and has a divergence in equilibrium with the load $f$, since this is not explicitly stipulated in the weak formulation (4.2). Congruently, it immediately follows (cf. Proposition 4.4 below) that the three properties (4.4) together with the constitutive relation (4.3) actually define the weak solution $p$ of (4.2) in a unique way. This mathematically describes the intrinsic properties of the exact potential $p$ and flux $u$: $p$ has to be trace continuous and $u$ has to be related by the constitutive relation (4.3) to it, normal-trace continuous, and in equilibrium with the load.

**Remark 4.3** (Physical properties of the weak solution). Scalar-valued physical variables such as the potential (pressure) $p$ are naturally trace continuous, whereas vector-valued physical variables such as the flux (Darcy velocity) $u = -\nabla p$ are naturally normal-trace continuous, cf. Figures 1 and 2. Moreover, $u$ has to be in equilibrium with the load (locally mass conservative) which amounts to $\nabla \cdot u = f$. This restates that the three properties (4.4) together with the constitutive relation (4.3) characterize the exact solution.
## 4.4 Approximate solution and spatial discretization error

Let $u_h \in L^2(\Omega)$ be an arbitrary function that we think of a (numerical) approximation of the exact flux $u$ from (4.2)–(4.3). We will give a specific example for the cell-centered finite volume discretizations on a simplicial mesh in Section 4.8 below, but we want to proceed abstractly for the moment, since there is no link to a mesh or to the discrete world of piecewise polynomials at this stage. We will call the $L^2(\Omega)$-distance of $u_h$ to $u$,

$$ \|u - u_h\|, \qquad (4.5) $$

the spatial discretization error.

## 4.5 Prager–Synge equality

We now state the following important result dating back to Prager and Synge [216]:

**Proposition 4.4 (Prager–Synge equality).** Let $p \in H_0^1(\Omega)$ be the weak solution of (4.2) and let $\zeta_h \in H_0^1(\Omega)$ and $\sigma_h \in H(\text{div}, \Omega)$ with $\nabla \cdot \sigma_h = f$ be arbitrary. Then

$$ \|\nabla(p - \zeta_h)\|^2 + \|\nabla p + \sigma_h\|^2 = \|\nabla \zeta_h + \sigma_h\|^2. \qquad (4.6) $$

*Proof.* Adding and subtracting $\nabla p$, we develop

$$ \begin{aligned} \|\nabla \zeta_h + \sigma_h\|^2 &= \|\nabla(\zeta_h - p) + \nabla p + \sigma_h\|^2 \\ &= \|\nabla(\zeta_h - p)\|^2 + \|\nabla p + \sigma_h\|^2 + 2(\nabla(\zeta_h - p), \nabla p + \sigma_h). \end{aligned} $$

Note from Proposition 4.1 that $\nabla p \in H(\text{div}, \Omega)$ with $\nabla \cdot (\nabla p) = -f$. Thus $(\nabla p + \sigma_h) \in H(\text{div}, \Omega)$ and in particular $\nabla \cdot (\nabla p + \sigma_h) = 0$. Thus, using that $\zeta_h - p \in H_0^1(\Omega)$, the Green theorem (3.3) gives

$$ (\underbrace{\nabla p + \sigma_h}_{\in H(\text{div}, \Omega)}, \underbrace{\nabla(\zeta_h - p)}_{\in H_0^1(\Omega)}) \stackrel{(3.3)}{=} -(\underbrace{\nabla \cdot (\nabla p + \sigma_h)}_{=0}, \zeta_h - p) = 0, $$

whence the assertion follows.

□

**Remark 4.5 (Prager–Synge equality (4.6)).** We stress that (4.6) is an equality, where on the right-hand side, there are, in practice, two known, discrete objects from the piecewise polynomial subspaces of $H_0^1(\Omega)$ and $H(\text{div}, \Omega)$, whereas the left-hand side are their spatial discretization errors with respect to the weak solution of (4.2). In other words, if we want to trace a simultaneous potential and flux error, then we can compute it, there is no (a posteriori) estimate.

From (4.6), we immediately observe (cf. Theorem 6.1 and Corollary 6.6 in [250] and the references therein):

**Corollary 4.6 (Prager–Synge (in)equality).** Let $p \in H_0^1(\Omega)$ be the solution of (4.2), let $u = -\nabla p$ from (4.3), and let $u_h \in H(\text{div}, \Omega)$ with $\nabla \cdot u_h = f$ be arbitrary. Then

$$ \|u - u_h\| = \min_{v \in H_0^1(\Omega)} \|u_h + \nabla v\|, \qquad (4.7) $$

so that in particular

$$ \|u - u_h\| \le \|u_h + \nabla \zeta_h\| \qquad (4.8) $$

for an arbitrary $\zeta_h \in H_0^1(\Omega)$.

*Proof.* The second assertion (4.8) follows from (4.6) upon taking $\sigma_h = u_h$ and noticing that

$$ \|\nabla p + u_h\|^2 = \|\nabla \zeta_h + u_h\|^2 - \|\nabla(p - \zeta_h)\|^2 \le \|\nabla \zeta_h + u_h\|^2 $$

for any $\zeta_h \in H_0^1(\Omega)$, using $u = -\nabla p$ from (4.3). This immediately gives (4.7) with infimum in place of minimum and the $\le$ sign in place of the equality. Since we can take $\zeta_h = p$ and since $u = -\nabla p$, the minimum in (4.7) is actually attained and the right-hand side equals the left-hand side. □

**Remark 4.7 (Prager–Synge equality (4.7)).** According to (4.3) and (4.4), to be equal to the exact flux $-\nabla p$, any $u_h \in H(\text{div}, \Omega)$ with $\nabla \cdot u_h = f$ only misses being a minus weak gradient of some function from $H_0^1(\Omega)$. This is the meaning of the characterization (4.7). Notice that (4.7) expresses the distance of $u_h$ to $\nabla H_0^1(\Omega)$.
**Remark 4.8** (Prager-Synge inequality (4.8)). Having $\mathbf{u}_h \in \mathbf{H}(\text{div}, \Omega)$ with $\nabla \cdot \mathbf{u}_h = f$ at disposal, which we will see in Section 4.8 below is immediate in finite volume methods for the Poisson equation (4.1), when $f$ is piecewise constant, the Prager-Synge inequality (4.8) can be readily used to produce a guaranteed error upper bound also satisfying the other requirements formulated in Section 2.5. This will still be the situation for the singlephase steady linear Darcy flow in Section 5 below, but when we want to take into account iterative algebraic solvers in Section 6 and component fluxes together with nonlinear zero-order terms in Section 7 for nonlinear compositional Darcy flows, it may not be feasible to locally construct $\mathbf{u}_h \in \mathbf{H}(\text{div}, \Omega)$ with $\nabla \cdot \mathbf{u}_h = f$ as an approximate solution. There, we will in general only have $\mathbf{u}_h \in L^2(\Omega)$. In this case, we will reconstruct $\sigma_h^{k,i} \in \mathbf{H}(\text{div}, \Omega)$ with $\nabla \cdot \sigma_h^{k,i} = f$ or $\nabla \cdot \sigma_h^{k,i} = f + \rho_h^{k,i}$ with $\rho_h^{k,i}$ a remainder made as small as necessary.

## 4.6 Equivalence of the spatial discretization error to the flux nonconformity (dual norm of the residual) plus potential nonconformity

Following [213, Theorem 3.1], [125, Theorem 3.3], and the references therein, we now extend Corollary 4.6 to the general case $\mathbf{u}_h \in L^2(\Omega)$.

**Theorem 4.9** (Error characterization for $\mathbf{u}_h \in L^2(\Omega)$). Let $p \in H_0^1(\Omega)$ be the weak solution of (4.2), let $\mathbf{u} = -\nabla p$ from (4.3), and let $\mathbf{u}_h \in L^2(\Omega)$ be arbitrary. Then

$$ \|\mathbf{u} - \mathbf{u}_h\|^2 = \underbrace{\min_{\substack{\mathbf{v} \in \mathbf{H}(\text{div}, \Omega) \\ \nabla \cdot \mathbf{v} = f}} \|\mathbf{u}_h - \mathbf{v}\|^2}_{\max_{\substack{\varphi \in H_0^1(\Omega); \|\nabla \varphi\| = 1 \\ \text{dual norm of the residual of } \mathbf{u}_h}} \{(f, \varphi) + (\mathbf{u}_h, \nabla \varphi)\}^2}, + \underbrace{\min_{\substack{\mathbf{v} \in H_0^1(\Omega) \\ \text{nonconformity} \\ \text{(distance of } \mathbf{u}_h \text{ to } \nabla H_0^1(\Omega)\}} \|\mathbf{u}_h + \nabla \mathbf{v}\|^2}_{\text{distance of } \mathbf{u}_h \text{ to } \nabla H_0^1(\Omega)} . \quad (4.9) $$

**Remark 4.10** (Error characterization (4.9)). Relation (4.9) extends (4.7) from $\mathbf{u}_h \in \mathbf{H}(\text{div}, \Omega)$ with $\nabla \cdot \mathbf{u}_h = f$ to any $\mathbf{u}_h \in L^2(\Omega)$. The supplementary arising term is naturally the distance of $\mathbf{u}_h$ to $\mathbf{H}(\text{div}, \Omega)$ under the constraint of the divergence being equal to $f$.

**Remark 4.11** (Dual norm of the residual). The first term in (4.9) has an equivalent expression investigating “how much, for a test function $\varphi \in H_0^1(\Omega)$, $\mathbf{u}_h$ misses to fulfill the weak formulation (4.2) when put in place of $-\nabla p$”, which is related to the residual of (4.2) for $\varphi$. The overall maximum over all $\varphi$ from $H_0^1(\Omega)$ with $\|\nabla \varphi\| = 1$ is then the dual norm of the residual.

**Remark 4.12** (Potential nonconformity). The second term in (4.9) expresses the distance of $\mathbf{u}_h$ to $\nabla H_0^1(\Omega)$, cf. Remark 4.7.

*Proof of Theorem 4.9.* Let us define a function $\zeta \in H_0^1(\Omega)$ by

$$ (\nabla \zeta, \nabla v) = -(\mathbf{u}_h, \nabla v) \quad \forall v \in H_0^1(\Omega). \quad (4.10) $$

There exists one and only one $\zeta$ by the Riesz representation theorem; indeed, recall from Section 3.2.1 that the left-hand side of (4.10) is a scalar product on $H_0^1(\Omega)$ and notice that from $|-(\mathbf{u}_h, \nabla v)| \le \|\mathbf{u}_h\|\|\nabla v\|$, the right-hand side of (4.10) is a bounded linear form on $H_0^1(\Omega)$ ($\mathbf{u}_h$ is here a fixed datum). The function $\zeta$ can be seen as the orthogonal projection of the approximate solution $\mathbf{u}_h$ onto $\nabla H_0^1(\Omega)$. With the aid of $\zeta$, we can thus write the Pythagorean equality

$$ \|\nabla p + \mathbf{u}_h\|^2 = \|\nabla(p - \zeta)\|^2 + \|\nabla \zeta + \mathbf{u}_h\|^2. \quad (4.11) $$

Indeed,

$$ \|\nabla p + \mathbf{u}_h\|^2 = \|\nabla p - \nabla \zeta + \nabla \zeta + \mathbf{u}_h\|^2 = \|\nabla(p - \zeta)\|^2 + \|\nabla \zeta + \mathbf{u}_h\|^2 + 2(\nabla(p - \zeta), \nabla \zeta + \mathbf{u}_h), $$

and the last term in the above expression vanishes in view of the orthogonality (4.10), since $p - \zeta$ can be taken as a test function $v \in H_0^1(\Omega)$ in (4.10). We continue in three steps.

1) Since $\zeta$ can be seen as a projection of $\mathbf{u}_h$,

$$ \|\nabla \zeta + \mathbf{u}_h\|^2 = \min_{v \in H_0^1(\Omega)} \|\nabla v + \mathbf{u}_h\|^2. \quad (4.12) $$

Indeed, from (4.11) used for an arbitrary function $v \in H_0^1(\Omega)$ in place of $p \in H_0^1(\Omega)$, we see

$$ \|\nabla v + \mathbf{u}_h\|^2 = \|\nabla(v - \zeta)\|^2 + \|\nabla \zeta + \mathbf{u}_h\|^2. $$
Therefrom, we get

$$ \|\nabla \zeta + \mathbf{u}_h\|^2 = \|\nabla v + \mathbf{u}_h\|^2 - \|\nabla(v - \zeta)\|^2 \le \|\nabla v + \mathbf{u}_h\|^2 \quad \forall v \in H_0^1(\Omega), $$

and there is an equality for $v = \zeta$. This shows that the second term in (4.11) indeed takes the form needed in (4.9).

2) For the first term in (4.11), we first use that $p - \zeta \in H_0^1(\Omega)$. Thus, the dual norm characterization and (4.10) give

$$ \|\nabla(p - \zeta)\| = \max_{\substack{\varphi \in H_0^1(\Omega) \\ \|\nabla\varphi\| = 1}} (\nabla(p - \zeta), \nabla\varphi) \stackrel{(4.10)}{=} \max_{\substack{\varphi \in H_0^1(\Omega) \\ \|\nabla\varphi\| = 1}} (\nabla p + \mathbf{u}_h, \nabla\varphi). \quad (4.13) $$

Let now $\varphi \in H_0^1(\Omega)$ with $\|\nabla\varphi\| = 1$ be fixed. Using the characterization (4.2) of the weak solution, we have

$$ (\nabla p + \mathbf{u}_h, \nabla\varphi) \stackrel{(4.2)}{=} (f, \varphi) + (\mathbf{u}_h, \nabla\varphi). \quad (4.14) $$

Thus, we have (4.9) with the first term in the form expressed in the underbrace.

3) We are left to show the equivalence

$$ \min_{\substack{\mathbf{v} \in \mathcal{H}(\text{div}, \Omega) \\ \nabla \cdot \mathbf{v} = f}} \|\mathbf{u}_h - \mathbf{v}\| = \max_{\substack{\varphi \in H_0^1(\Omega); \\ \|\nabla\varphi\| = 1}} \{(f, \varphi) + (\mathbf{u}_h, \nabla\varphi)\}. \quad (4.15) $$

3a) First, to motivate, we show that the right-hand side of (4.15) is bounded by the left-hand side one. Let $\mathbf{v} \in \mathcal{H}(\text{div}, \Omega)$ such that $\nabla \cdot \mathbf{v} = f$ be arbitrary. Then the Green theorem (3.3) gives

$$ (f, \varphi) + (\mathbf{u}_h, \nabla\varphi) = (\nabla \cdot \mathbf{v}, \varphi) + (\mathbf{u}_h, \nabla\varphi) = (\mathbf{u}_h - \mathbf{v}, \nabla\varphi). $$

Consequently, by the Cauchy–Schwarz inequality,

$$ \max_{\substack{\varphi \in H_0^1(\Omega); \\ \|\nabla\varphi\| = 1}} \{(f, \varphi) + (\mathbf{u}_h, \nabla\varphi)\} \le \min_{\substack{\mathbf{v} \in \mathcal{H}(\text{div}, \Omega) \\ \nabla \cdot \mathbf{v} = f}} \|\mathbf{u}_h - \mathbf{v}\|. \quad (4.16) $$

3b) We now show the equality (4.15). The argument of the constrained minimization in (4.15) is

$$ \boldsymbol{\sigma} := \arg \min_{\substack{\mathbf{v} \in \mathcal{H}(\text{div}, \Omega) \\ \nabla \cdot \mathbf{v} = f}} \|\mathbf{u}_h - \mathbf{v}\| $$

and is characterized by the Euler–Lagrange conditions as a function $\boldsymbol{\sigma} \in \mathcal{H}(\text{div}, \Omega)$ with $\nabla \cdot \boldsymbol{\sigma} = f$ such that

$$ (\boldsymbol{\sigma}, \mathbf{v}) = (\mathbf{u}_h, \mathbf{v}) \quad \forall \mathbf{v} \in \mathcal{H}(\text{div}, \Omega) \text{ with } \nabla \cdot \mathbf{v} = 0. $$

Imposing the constraint $\nabla \cdot \boldsymbol{\sigma} = f$ equivalently with test functions $q \in L^2(\Omega)$ and introducing the Lagrange multiplier $r \in L^2(\Omega)$, using that $\nabla \cdot \mathcal{H}(\text{div}, \Omega) = L^2(\Omega)$, the above problem is equivalent to finding $\boldsymbol{\sigma} \in \mathcal{H}(\text{div}, \Omega)$ and $r \in L^2(\Omega)$ such that

$$ (\boldsymbol{\sigma}, \mathbf{v}) - (r, \nabla \cdot \mathbf{v}) = (\mathbf{u}_h, \mathbf{v}) \quad \forall \mathbf{v} \in \mathcal{H}(\text{div}, \Omega), \quad (4.17a) $$

$$ (\nabla \cdot \boldsymbol{\sigma}, q) = (f, q) \quad \forall q \in L^2(\Omega). \quad (4.17b) $$

Now, (4.17a) implies by Definition 3.1 of the weak partial derivative that $r \in H_0^1(\Omega)$ with $\nabla r = -(\boldsymbol{\sigma} - \mathbf{u}_h)$. Consequently, by the Green theorem (3.3),

$$ \begin{aligned} \min_{\substack{\mathbf{v} \in \mathcal{H}(\text{div}, \Omega) \\ \nabla \cdot \mathbf{v} = f}} \|\mathbf{u}_h - \mathbf{v}\| &= \|\mathbf{u}_h - \boldsymbol{\sigma}\| = \|\nabla r\| = \max_{\substack{\varphi \in H_0^1(\Omega) \\ \|\nabla\varphi\| = 1}} (\nabla r, \nabla\varphi) \\ &= \max_{\substack{\varphi \in H_0^1(\Omega) \\ \|\nabla\varphi\| = 1}} (-\boldsymbol{\sigma} + \mathbf{u}_h, \nabla\varphi) = \max_{\substack{\varphi \in H_0^1(\Omega) \\ \|\nabla\varphi\| = 1}} \{(f, \varphi) + (\mathbf{u}_h, \nabla\varphi)\}, \end{aligned} \quad (4.18) $$

which is (4.15). □
<figure><img src="image_11.png" /><figcaption>Figure 10: Admissible triangular mesh $\mathcal{T}_h$ and notation for the cell-centered finite volume scheme</figcaption></figure>

## 4.7 Cell-centered finite volume discretizations on a simplicial mesh

Consider now a simplicial mesh $\mathcal{T}_h$ of the domain $\Omega$ as defined in Section 3.3. We additionally assume that 1) there exists a point $x_K$ associated with each simplex $K \in \mathcal{T}_h$; 2) there exists a point $x_{K,\sigma}$ associated with each face $\sigma$ lying in the boundary $\partial\Omega$, $\sigma \in \mathcal{F}_K \cap \mathcal{F}_h^{\text{ext}}$, $K \in \mathcal{T}_h$, lying in the interior of the face $\sigma$; 3) all the points $x_K$ and $x_{K,\sigma}$ are distinct; 4) the straight line connecting $x_K$ and $x_L$ for two neighboring elements $K$ and $L$ is orthogonal to their common face $\sigma_{K,L} = \partial K \cap \partial L \in \mathcal{F}_h^{\text{int}}$; 5) the straight line connecting $x_K$ and $x_{K,\sigma}$ is orthogonal to $\sigma$ for $\sigma \in \mathcal{F}_K \cap \mathcal{F}_h^{\text{ext}}$. The points $x_K$ typically, but not necessarily, lie in the interior of the element $K \in \mathcal{T}_h$. In particular, Delaunay meshes satisfy this requirement of “admissibility” in the sense of [126, Definition 9.1]. An illustration is provided in Figure 10.

The cell-centered finite volume scheme for the Poisson equation (4.1) on the simplicial mesh $\mathcal{T}_h$ reads: find the real values $p_K$, $K \in \mathcal{T}_h$, the approximations to the mean values of $p$ in the mesh elements $K$, such that

$$ \sum_{\sigma \in \mathcal{F}_K} U_{K,\sigma} = (f, 1)_K \quad \forall K \in \mathcal{T}_h, \qquad (4.19a) $$

where $U_{K,\sigma} \in \mathbb{R}$ for each face $\sigma \in \mathcal{F}_K$ approximates the normal (out)flux $\langle u \cdot n_K, 1 \rangle_{\sigma}$ from $K$ over the face $\sigma$ (recall the notation (3.12)) by

$$ U_{K,\sigma} := -\frac{|\sigma_{K,L}|}{d_{K,L}}(p_L - p_K) \quad \sigma = \sigma_{K,L} \in \mathcal{F}_h^{\text{int}}, \qquad (4.19b) $$

$$ U_{K,\sigma} := -\frac{|\sigma|}{d_{K,\sigma}}(0 - p_K) \quad \sigma \in \mathcal{F}_K \cap \mathcal{F}_h^{\text{ext}}. \qquad (4.19c) $$

Recall that $|\sigma|$ is the $(d-1)$-dimensional measure of the face $\sigma$. Moreover, $d_{K,L} = |x_K - x_L|$ is the Euclidean distance of the points $x_K$ and $x_L$ for an interior face $\sigma = \sigma_{K,L} \in \mathcal{F}_h^{\text{int}}$, unless the positions of $x_K$ and $x_L$ are interchanged with respect to the positions of $K$ and $L$, where we take $d_{K,L} = -|x_K - x_L|$. Since we suppose $x_K$ and $x_L$ distinct, $d_{K,L} \neq 0$. Similarly, $d_{K,\sigma} = |x_K - x_{K,\sigma}|$ is the Euclidean distance of the points $x_K$ and $x_{K,\sigma}$ for a boundary face $\sigma \in \mathcal{F}_K \cap \mathcal{F}_h^{\text{ext}}$, unless $x_K$ lies outside of $K$ past $\sigma$, where we take $d_{K,\sigma} = -|x_K - x_{K,\sigma}|$. Figure 10 provides an illustration.

**Remark 4.13** (Flux conservation and equilibrium in finite volumes on simplicial meshes). From (4.19b), there holds

$$ U_{K,\sigma} = -U_{L,\sigma} \quad \forall \sigma = \sigma_{K,L} \in \mathcal{F}_h^{\text{int}}. \qquad (4.20) $$

This is the flux conservation in finite volumes: what flows from the element $K$ out through the face $\sigma_{K,L}$ flows in the neighbor element $L$ through the same face $\sigma_{K,L}$, cf. Lemma 3.11 and definition (3.15). The requirement (4.19a) then stipulates the equilibrium of these face fluxes with the load $f$, mimicking $\nabla \cdot u = f$ from Proposition 4.1.

## 4.8 Flux reconstruction by lifting the face normal fluxes

Following [127], we now lift the normal fluxes $U_{K,\sigma}$ from (4.19b)-(4.19c) to form a vector-valued approximate flux field $u_h$ that is $H(\text{div}, \Omega)$-conforming and in equilibrium with the load $f$. More precisely, we construct $u_h$ in the piecewise polynomial Raviart-Thomas space $\mathcal{RT}_0(\mathcal{T}_h) \cap H(\text{div}, \Omega)$ from (3.15).
<figure><img src="image_12.png" /><figcaption>Figure 11: Finite volume face fluxes $U_{K,\sigma}$ (dashed black arrows) and flux reconstruction $u_h$ of Definition 4.14 (blue arrows)</figcaption></figure>

Recalling the definition (3.10) of the space $RT_0(K)$ for a simplex $K$ and the notation from Figure 6, we start with:

**Definition 4.14 (Flux reconstruction).** Let the cell-centered finite volume discretization on a simplicial mesh be given by (4.19). Then, on each simplex $K \in \mathcal{T}_h$, define

$$ u_h|_K \in RT_0(K), \quad \langle u_h \cdot n_K, 1 \rangle_\sigma = U_{K,\sigma} \quad \forall \sigma \in \mathcal{F}_K, \qquad (4.21) $$

or, equivalently, using the Raviart–Thomas basis functions $v_\sigma$ from (3.13),

$$ (u_h|_K)(\mathbf{x}) := \sum_{\sigma \in \mathcal{F}_K} U_{K,\sigma} v_\sigma(\mathbf{x}) = \sum_{\sigma \in \mathcal{F}_K} U_{K,\sigma} \frac{1}{d|K|} (\mathbf{x} - \mathbf{a}_{K,\sigma}), \quad \mathbf{x} \in K. \qquad (4.22) $$

An illustration of the flux reconstruction $u_h$ from the finite volume face fluxes $U_{K,\sigma}$ is provided in Figure 11. Crucially, we have:

**Lemma 4.15 (Flux reconstruction).** Let $u_h$ be given by Definition 4.14. Then

$$ u_h \in RT_0(\mathcal{T}_h) \cap H(\mathrm{div}, \Omega), \quad (\nabla \cdot u_h)|_K = \frac{(f, 1)_K}{|K|} \quad \forall K \in \mathcal{T}_h. \qquad (4.23) $$

*Proof.* From the finite volume conservation (4.20) and definition (4.21), we have the normal trace continuity needed in (3.15); this yields the first property in (4.23). As for the second property in (4.23), we write, using the Green theorem, definition (4.21), and the finite volume scheme (4.19a)

$$ (\nabla \cdot u_h, 1)_K = \sum_{\sigma \in \mathcal{F}_K} \langle u_h \cdot n_K, 1 \rangle_\sigma \stackrel{(4.21)}{=} \sum_{\sigma \in \mathcal{F}_K} U_{K,\sigma} \stackrel{(4.19a)}{=} (f, 1)_K. $$

The conclusion follows since $\nabla \cdot u_h$ is constant by virtue of (3.11). $\square$

## 4.9 Potential reconstruction by elementwise postprocessing and averaging

Supposing for the moment that $f$ is piecewise constant, so that $u_h \in H(\mathrm{div}, \Omega)$ with $\nabla \cdot u_h = f$ from (4.23), we will be in position to use Corollary 4.6 to estimate the error between the finite volume flux reconstruction $u_h$ of Definition 4.14 and the exact flux $u = -\nabla p$. Inequality (4.8) tells us that we will need a suitable scalar-valued object in the primal space $H_0^1(\Omega)$ whose negative gradient is as close as possible to $u_h$. Since the piecewise constant function given by the finite volume values $p_K$ on each $K \in \mathcal{T}_h$ cannot be used directly, as it does not lie in the space $H_0^1(\Omega)$, we proceed in two steps.

Following [247, Section 4.1] and [249, Section 3.2], we first create a scalar-valued piecewise polynomial whose negative gradient equals $u_h$ and whose mean values are given by $p_K$ on each mesh element $K \in \mathcal{T}_h$. Recalling the piecewise polynomial space $P_k(\mathcal{T}_h)$ from (3.8), we define:

**Definition 4.16 (Potential postprocessing).** Let $p_K$, $K \in \mathcal{T}_h$, be the finite volume approximate solution of (4.19). Let $u_h$ be given by Definition 4.14. Define the potential postprocessing $\tilde{p}_h$ as a piecewise quadratic polynomial on the simplicial mesh $\mathcal{T}_h$, $\tilde{p}_h \in P_2(\mathcal{T}_h)$, given by

$$ -\nabla \tilde{p}_h|_K = u_h|_K, \quad \frac{(\tilde{p}_h, 1)_K}{|K|} = p_K \quad \forall K \in \mathcal{T}_h. \qquad (4.24) $$
<figure><img src="image_13.png" /><figcaption>Figure 12: Potential postprocessing $\tilde{p}_h \in \mathcal{P}_2(\mathcal{T}_h)$ (left) and potential reconstruction by averaging $\zeta_h \in \mathcal{P}_2(\mathcal{T}_h) \cap H^1(\Omega)$ (right)</figcaption></figure>

It is easy to verify that any object from $\mathcal{RT}_0(K)$ for a simplex $K$ is a gradient of some polynomial which is at most second-order, since fields in $\mathcal{RT}_0(K)$ are curl-free and at most affine. Thus (4.24) is well-posed.

Unfortunately, in general $\tilde{p}_h$ from Definition 4.16 does not lie in $H_0^1(\Omega)$; if this was the case, then, actually, by Proposition 4.4 or Theorem 4.9, $\tilde{p}_h$ would be the weak solution $p$. An illustration of a prototypical situation is provided in Figure 12, left. For this reason, we continue with a second step, in the spirit of averaging operators in Durán and Padra [112], Achdou et al. [3], Karakashian and Pascal [157], Ainsworth [12], or Burman and Ern [68]. Recall the Lagrange spaces and nodes of Section 3.5.2 and Figure 4.

**Definition 4.17** (Potential reconstruction by averaging). Let $\tilde{p}_h \in \mathcal{P}_2(\mathcal{T}_h)$ be given by Definition 4.16. We call a potential reconstruction by averaging a piecewise quadratic and continuous polynomial $\zeta_h \in \mathcal{P}_2(\mathcal{T}_h) \cap H_0^1(\Omega)$ that is prescribed by its values in the Lagrange nodes of the space $\mathcal{P}_2(\mathcal{T}_h) \cap H_0^1(\Omega)$ by

$$ \zeta_h(\mathbf{x}) := \frac{1}{|\mathcal{T}_h|} \sum_{K \in \mathcal{T}_h} \tilde{p}_h|_K(\mathbf{x}) \quad \mathbf{x} \text{ is a Lagrange node of } \mathcal{P}_2(\mathcal{T}_h) \cap H_0^1(\Omega) \text{ included in } \Omega, \quad (4.25a) $$

$$ \zeta_h(\mathbf{x}) := 0 \quad \mathbf{x} \text{ is a Lagrange node of } \mathcal{P}_2(\mathcal{T}_h) \cap H_0^1(\Omega) \text{ included in } \partial\Omega, \quad (4.25b) $$

where $\mathcal{T}_h$ denotes the set of elements of the mesh $\mathcal{T}_h$ that contain the point $\mathbf{x}$ and $|\mathcal{T}_h|$ is the cardinality (number of elements) of this set.

We simply average all the (typically) different values that $\zeta_h$ takes in the Lagrange nodes. An illustration is given in Figure 12, right (where the 0 values at the boundary $\partial\Omega$ are not imposed).

## 4.10 A guaranteed a posteriori error estimate

With the above developments, we are now in a position to present a guaranteed a posteriori error estimate for the cell-centered finite volume discretization (4.19) of the Poisson equation (4.1) on a simplicial mesh $\mathcal{T}_h$. We follow [249, Theorem 4.1].

**Theorem 4.18** (A guaranteed a posteriori error estimate).

Let $p \in H_0^1(\Omega)$ be the exact potential given by (4.2) and let $\mathbf{u}$ be the exact flux given by (4.3). Let the cell-centered finite volume discretization be given by (4.19) and let the approximate flux $\mathbf{u}_h \in \mathcal{RT}_0(\mathcal{T}_h) \cap \mathbf{H}(\text{div}, \Omega)$ be constructed following Definition 4.14. Let the potential postprocessing $\tilde{p}_h \in \mathcal{P}_2(\mathcal{T}_h)$ be given by Definition 4.16 and the potential reconstruction by averaging $\zeta_h \in \mathcal{P}_2(\mathcal{T}_h) \cap H_0^1(\Omega)$ by Definition 4.17. Then there holds

$$ \|\mathbf{u} - \mathbf{u}_h\| \le \eta := \left\{ \sum_{K \in \mathcal{T}_h} \eta_K^2 \right\}^{1/2}, \quad (4.26a) $$
where

$$ \eta_K^2 := \|\mathbf{u}_h + \nabla \zeta_h\|_K^2 + \frac{h_K^2}{\pi^2} \left\| f - \frac{(f, 1)_K}{|K|} \right\|_K^2. \quad (4.26b) $$

*Proof.* If $f$ is piecewise constant, then $\mathbf{u}_h \in \mathbf{H}(\mathrm{div}, \Omega)$ with $\nabla \cdot \mathbf{u}_h = f$ from (4.23). In this case, we merely employ inequality (4.8) from Corollary 4.6.

In the general case $f \in L^2(\Omega)$, we rather use equality (4.9) from Theorem 4.9. The second term on the right-hand side of (4.9) gives rise to $\|\mathbf{u}_h + \nabla \zeta_h\|^2$. As for the first one, let $\varphi \in H_0^1(\Omega)$ with $\|\nabla \varphi\| = 1$ be fixed. By the Green theorem (3.3), the flux reconstruction property (4.23), the Cauchy-Schwarz inequality, and the Poincaré inequality

$$ \left\| \varphi - \frac{(\varphi, 1)_K}{|K|} \right\|_K \le \frac{h_K}{\pi} \|\nabla \varphi\|_K, \quad \varphi \in H^1(K), \quad (4.27) $$

cf. [212, 37], we see

$$ \begin{aligned} (f, \varphi) + (\mathbf{u}_h, \nabla \varphi) &\stackrel{(3.3)}{=} (f - \nabla \cdot \mathbf{u}_h, \varphi) = \sum_{K \in \mathcal{T}_h} (f - \nabla \cdot \mathbf{u}_h, \varphi)_K \stackrel{(4.23)}{=} \sum_{K \in \mathcal{T}_h} \left( f - \frac{(f, 1)_K}{|K|}, \varphi \right)_K \\ &= \sum_{K \in \mathcal{T}_h} \left( f - \frac{(f, 1)_K}{|K|}, \varphi - \frac{(\varphi, 1)_K}{|K|} \right)_K \stackrel{(4.27)}{\le} \sum_{K \in \mathcal{T}_h} \left( \left\| f - \frac{(f, 1)_K}{|K|} \right\|_K \frac{h_K}{\pi} \|\nabla \varphi\|_K \right) \\ &\le \left\{ \sum_{K \in \mathcal{T}_h} \left( \frac{h_K^2}{\pi^2} \left\| f - \frac{(f, 1)_K}{|K|} \right\|_K^2 \right) \right\}^{1/2} \underbrace{\|\nabla \varphi\|}_{=1} = \left\{ \sum_{K \in \mathcal{T}_h} \left( \frac{h_K^2}{\pi^2} \left\| f - \frac{(f, 1)_K}{|K|} \right\|_K^2 \right) \right\}^{1/2}. \end{aligned} $$

Thus (4.26) follows. $\square$

**Remark 4.19** (Local efficiency). Local efficiency of the estimate of (4.26) as per Section 2.5.2 holds true, see [249, Theorem 4.2].

**Remark 4.20** (Rectangular meshes). All the above results are in [249] also presented on rectangular meshes.

## 4.11 Numerical experiments

In this Section, we illustrate the application of the a posteriori error estimates of the form of Theorem 4.18 to the finite volume method on two test cases in two space dimensions.

### 4.11.1 Regular solution

In this test, we consider a model problem from [192] of form (4.1) with $\Omega = (0, 1) \times (0, 1)$ and the load term $f$ given such that the exact solution reads

$$ p(x, y) = \frac{10^2}{4} x(1-x)y(1-y) \exp(-100 \left( (x-0.75)^2 + (y-0.75)^2 \right)). \quad (4.28) $$

This solution is regular but contains a localized exponential peak. The numerical tests are performed on a sequence of uniformly refined rectangular meshes $\mathcal{T}_0, \mathcal{T}_1, \dots, \mathcal{T}_J, J = 3$, cf. Remarks 3.12 and 4.20. Figure 13 illustrates the approximate solution given by the values $p_K, K \in \mathcal{T}_h$, the elementwise errors $\|\mathbf{u} - \mathbf{u}_h\|_K$, and the corresponding a posteriori error estimators $\eta_K$ arising from Theorem 4.18 on mesh $\mathcal{T}_3$. We observe that the error is localized in a circular zone around the peak and that the actual and predicted error distributions match very closely.

In the left part of Figure 14, we plot the relative exact errors $\|\mathbf{u} - \mathbf{u}_h\|/\|\mathbf{u}_h\|$ on each of the meshes $\mathcal{T}_0-\mathcal{T}_3$ together with their relative a posteriori error estimates $\eta/\|\mathbf{u}_h\|$ from Theorem 4.18, with respect to the number of mesh elements/number of unknowns. We indeed observe that the estimators give a guaranteed upper bound on the errors. In the right part of Figure 14 we then plot the effectivity indices defined as per (2.22) by the ratio of the estimator to the error,

$$ I_{\text{eff}} := \frac{\eta}{\|\mathbf{u} - \mathbf{u}_h\|}. \quad (4.29) $$

We observe that they are remarkably close to one and tend to the optimal value of one with mesh refinement, showing asymptotic exactness as per Section 2.5.4.
<figure><img src="image_14.png" /><figcaption>Figure 13: [Section 4.11.1, mesh $\mathcal{T}_3$] The approximate solution given by the values $p_K$, $K \in \mathcal{T}_h$, (left), the exact errors $\|u - u_h\|_K$ (middle), and the error estimators $\eta_K$ from Theorem 4.18 (right).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14: [Section 4.11.1, uniformly refined meshes $\mathcal{T}_0-\mathcal{T}_3$] Relative errors $\|u - u_h\|/\|u_h\|$ and relative estimates $\eta/\|u_h\|$ from Theorem 4.18 (left), effectivity indices $I_{\text{eff}}$ from (4.29) (right)</figcaption></figure>

## 4.11.2 Singular solution

We consider a benchmark test case defined on the L-shaped domain

$$ \Omega = (-1, 1) \times (-1, 1) \setminus (-1, 0] \times (-1, 0] $$

with the exact solution

$$ p(r, \theta) = r^{\frac{2}{3}} \sin \left( \frac{2}{3}\theta + \frac{3}{2}\pi \right), $$

where $(r, \theta)$ are the polar coordinates. The Dirichlet boundary condition is partly inhomogeneous here and given by the value of the exact solution on $\partial\Omega$. We illustrate in Figure 15 the approximate solution given by the values $p_K$, $K \in \mathcal{T}_h$, the elementwise errors $\|u - u_h\|_K$, and the corresponding a posteriori error estimators $\eta_K$ from Theorem 4.18. We observe that a posteriori error estimators $\eta_K$ detect perfectly the error related in the vicinity of the singularity at $(0, 0)$.

The convergence histories on the sequence of the uniformly-refined meshes $\mathcal{T}_0-\mathcal{T}_2$ are depicted in the left part of Figure 16. The right part of Figure 16 then displays the corresponding effectivity indices from (4.29). We see that, for this test case, the effectivity indices are still close to the optimal value of one, even with singular solution, though no asymptotic exactness is observed here.

## 4.12 Bibliographic resources

One of the first a posteriori error analyses for finite volume methods was performed by Angermann [20]. Later, Achdou et al. [3], Afif et al. [6], Lazarov and Tomov [178], Carstensen et al. [73], and Amaziane et
<figure><img src="image_16.png" /><figcaption>Figure 15: [Section 4.11.2, mesh $\mathcal{T}_2$] The approximate solution given by the values $p_K$, $K \in \mathcal{T}_h$, (left), the exact errors $\|u - u_h\|_K$ (middle), and the error estimators $\eta_K$ from Theorem 4.18 (right).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 16: [Section 4.11.2, uniformly refined meshes $\mathcal{T}_0-\mathcal{T}_2$] Relative errors $\|u - u_h\|/\|u_h\|$ and relative estimates $\eta/\|u_h\|$ from Theorem 4.18 (left), effectivity indices $I_{\text{eff}}$ from (4.29) (right)</figcaption></figure>

al. [18] also considered the so-called vertex-centered finite volume discretizations, closely related to lowest-order finite elements. Cell-centered finite volumes of the form (4.19a) were then addressed in Agouzal and Oudin [8], Achdou et al. [3], Nicaise [196], Vohralík [249], and Erath and Praetorius [116].

# 5 Steady linear pure diffusion problems. Polytopal meshes and inexpensive implementation and evaluation

We describe here first a model steady linear problem, the Darcy flow problem. This extends the Poisson equation (4.1) from Section 4 and forms the basis for complex multiphase compositional porous media flows as that of Section 7. We also focus on extending the finite volume discretizations on simplicial meshes from Section 4.7 to general polytopal meshes. Finally, we investigate how the guaranteed a posteriori error estimates of the form of Theorem 4.18 can be implemented and evaluated in the most inexpensive way, and this on polytopal meshes. In particular, the flux and potential reconstructions of Sections 4.8–4.9 will only be performed here virtually, avoiding their factual construction. Instead, the a posteriori error estimates will only need *one flux value per polytopal face* and *one potential value per vertex of a (virtual) simplicial submesh*. The former are directly available from the numerical scheme (these are the generalizations of the face normal fluxes $U_{K,\sigma}$ in (4.19)) and the latter will be trivially obtained by some postprocessing and averaging of the elementwise potential values generalizing $p_K$ in (4.19).
## 5.1 Singlephase steady linear Darcy flow

Consider the problem of finding $p: \Omega \to \mathbb{R}$ such that

$$ 
\begin{aligned}
-\nabla \cdot (\underline{\mathbf{K}} \nabla p) &= f && \text{in } \Omega, \\
p &= 0 && \text{on } \partial\Omega.
\end{aligned}
\qquad (5.1a) \\
(5.1b)
$$

Here $f$ is a source term and $\underline{\mathbf{K}}$ is a symmetric, positive definite, and bounded diffusion tensor with values in $\mathbb{R}^{d \times d}$; we suppose for simplicity here that both $f$ and $\underline{\mathbf{K}}$ are piecewise constant with respect to the polytopal mesh $\mathcal{T}_H$.

## 5.2 Weak solution and its properties

The Sobolev space $H_0^1(\Omega)$ from Definition 3.4 still gives a proper mathematical setting to define a unique solution $p$ of problem (5.1), which is $p \in H_0^1(\Omega)$ such that

$$ (\underline{\mathbf{K}}\nabla p, \nabla v) = (f, v) \quad \forall v \in H_0^1(\Omega). \qquad (5.2) $$

From the pressure head $p$, we can define the Darcy velocity

$$ \mathbf{u} := -\underline{\mathbf{K}}\nabla p; \qquad (5.3) $$

recall that we refer to $p$ as potential and to $\mathbf{u}$ as flux. As in Proposition 4.1, it follows from (5.2) that

$$ \mathbf{u} \in \mathbf{H}(\mathrm{div}, \Omega), \quad \nabla \cdot \mathbf{u} = f. \qquad (5.4) $$

Thus, just like in Remark 4.3, the scalar-valued pressure head $p$ is trace continuous, cf. Figure 1, whereas the vector-valued Darcy velocity $\mathbf{u}$ is normal-trace continuous, cf. Figure 2. Moreover, $\mathbf{u}$ is locally mass conservative, in equilibrium with the load, satisfying $\nabla \cdot \mathbf{u} = f$.

## 5.3 Spatial discretization error

Applying the abstract reflections of Section 2.4.2, we will evaluate the spatial discretization error as

$$ \|\mathbf{u} - \mathbf{u}_h\|_{\underline{\mathbf{K}}^{-\frac{1}{2}}} \qquad (5.5) $$

for an approximate Darcy velocity $\mathbf{u}_h$. Here, the energy norm on $\omega \subset \Omega$, for $v \in [L^2(\omega)]^d$, writes as

$$ \|v\|_{\underline{\mathbf{K}}^{-\frac{1}{2},\omega}} := \|\underline{\mathbf{K}}^{-1/2}v\|_{L^2(\omega)} = \left\{ \int_{\omega} |\underline{\mathbf{K}}^{-1/2}(x)v(x)|^2 dx \right\}^{\frac{1}{2}}; \qquad (5.6) $$

when $\omega = \Omega$, we as usual drop the subscript $\omega$.

## 5.4 Generic discretizations on a polytopal mesh

Let $\mathcal{T}_H$ be a polytopal mesh of $\Omega$ in the sense of Section 3.6. Define the load vector $\mathbf{F} := \{\mathbf{F}_K\}_{K \in \mathcal{T}_H} \in \mathbb{R}^{|\mathcal{T}_H|}$ by

$$ \mathbf{F}_K := (f, 1)_K \quad \forall K \in \mathcal{T}_H. \qquad (5.7) $$

We consider any numerical discretization of (5.1) that can be written under the following abstract lowest-order locally conservative form:

**Assumption 5.1** (Locally conservative discretization on a polytopal mesh). Find the algebraic vector $\mathbf{P} := \{\mathbf{P}_K\}_{K \in \mathcal{T}_H} \in \mathbb{R}^{|\mathcal{T}_H|}$ such that the flux balance

$$ \sum_{\sigma \in \mathcal{F}_K} \mathbf{U}_{\sigma} \mathbf{n}_K \cdot \mathbf{n}_{\sigma} = \mathbf{F}_K \quad \forall K \in \mathcal{T}_H \qquad (5.8) $$

is satisfied. Here $\mathbf{U} := \{\mathbf{U}_{\sigma}\}_{\sigma \in \mathcal{F}_H} \in \mathbb{R}^{|\mathcal{F}_H|}$, and $\mathbf{U}_{\sigma} \in \mathbb{R}$ for each face $\sigma \in \mathcal{F}_H$ approximates the normal flux $\langle \mathbf{u} \cdot \mathbf{n}_{\sigma}, 1 \rangle_{\sigma}$ over the face $\sigma$ and depends linearly on $\mathbf{P}$. The unknowns $\mathbf{P}_K \in \mathbb{R}$ for each element $K \in \mathcal{T}_H$ approximate the potential $p$ in the element $K$.
Assumption 5.1 is very generic (allows for any lowest-order locally conservative method). Note in particular that for the a posteriori error analysis performed below, we do not need to know how the flux unknowns U depend on the pressure head unknowns P. Recall that on admissible simplicial meshes, the easiest example for this relation is the two-point formula (4.19b)-(4.19c).

Many lowest-order locally conservative methods take a form comprised in Assumption 5.1 but more specific. We consider them separately under the following assumption, in the spirit of mimetic finite differences [64, Theorem 5.1], mixed finite elements [254, Theorems 7.2 and 7.3], and the unifying frame- works [110, 111]:

**Assumption 5.2** (Locally conservative saddle-point discretization on a polytopal mesh). Find the alge- braic vectors U := {U$_{σ}$}$_{σ}$$_{∈F}$$_{H}$ ∈ R |F$_{H}$ | and P := {P$_{K}$}$_{K}$$_{∈T}$$_{H}$ ∈ R |T$_{H}$ | such that

$$
\begin{pmatrix} \mathbf{A} & \mathbf{B}^{\mathrm{t}} \\ \mathbf{B} & 0 \end{pmatrix} \begin{pmatrix} \mathbf{U} \\ \mathbf{P} \end{pmatrix} = \begin{pmatrix} \mathbf{0} \\ -\mathbf{F} \end{pmatrix}, \qquad (5.9)
$$

where

1) the unknown U$_{σ}$ ∈ ℝ for each face σ ∈ F$_{H}$ approximates the normal flux 〈 u · n$_{σ}$ , 1 〉$_{σ}$ over the face σ and the unknown P$_{K}$ ∈ ℝ for each element K ∈ T$_{H}$ approximates the potential p in the element K ;

2) the matrix B ∈ ℝ |T$_{H}$|×|F$_{H}$ | has a full rank;

3) for each polytopal cell K ∈ T$_{H}$ and face σ ∈ F$_{H}$ , B$_{K,σ}$ = − n$_{K}$ · n$_{σ}$ if σ is a face of K , σ ∈ F$_{K}$ , and B$_{K,σ}$ = 0 otherwise;

4) the matrix A ∈ ℝ |F$_{H}$|×|F$_{H}$ | is invertible;

5) for faces σ, σ' ∈ F$_{H}$ , A$_{σ,σ'}$ = 0 if σ and σ' are not faces of the same element K ∈ T$_{H}$ , { σ, σ' } ⊂ F$_{K}$ for some K ∈ T$_{H}$ , and

$$
\mathbb{A}_{\sigma, \sigma'} = \sum_{K \in \mathcal{T}_H, \{\sigma, \sigma'\} \subset \mathcal{F}_K} \mathbf{n}_K \cdot \mathbf{n}_\sigma \mathbf{n}_K \cdot \mathbf{n}_{\sigma'} (\hat{\mathbb{A}}_K)_{\sigma, \sigma'} \quad (5.10)
$$

otherwise, where ̂ A$_{K}$ ∈ ℝ |F$_{K}$|×|F$_{K}$ | are symmetric and positive definite (the “element matrices” of the given method for each K ∈ T$_{H}$ );

6) on each polytopal cell K ∈ T$_{H}$ , there exists a lifting ˜ u$_{h}$ |$_{K}$ ∈ H (div , K ) of the face normal fluxes U ext K = { U$_{σ}$ }$_{σ}$$_{∈F}$$_{K}$ to the interior of K such that

$$
\langle \tilde{u}_h |_K \cdot \mathbf{n}_\sigma, 1 \rangle_\sigma = U_\sigma \quad \forall \sigma \in \mathcal{F}_K, \qquad (5.11a)
$$

$$
\nabla \cdot \tilde{u}_h |_K = f |_K, \qquad (5.11b)
$$

$$
\| \tilde{u}_h \|_{K^{-\frac{1}{2}, K}}^2 = (U_K^{\text{ext}})^{\mathrm{t}} \hat{A}_K U_K^{\text{ext}}. \qquad (5.11c)
$$

5.5 Face normal fluxes

In Theorem 5.18 below, we will need a vector of real values with one value U$_{σ}$ per face σ ∈ F$_{H,h}$ (approximate face normal fluxes):

**Definition 5.3** (Face normal fluxes) . Let a discretization scheme of the form of Assumption 5.1 or 5.2 be given, leading to the face normal fluxes U := { U$_{σ}$ }$_{σ}$$_{∈F}$$_{H}$ ∈ R |F$_{H}$ | . For each polytopal face σ ∈ F$_{H}$ and each simplicial subface σ' ∈ F$_{H,h}$ as per Section 3.6, define the face normal flux

$$
U_{\sigma'} := \frac{U_{\sigma}}{|\sigma'|} |\sigma'| \quad \forall \sigma' \in \mathcal{F}_{H,h}, \sigma' \subset \sigma, \sigma \in \mathcal{F}_{H} \qquad (5.12)
$$

and collect them in the element vectors U ext K ∈ ℝ |F ext K,h |

$$
U_K^{\text{ext}} := \{U_\sigma\}_{\sigma \in \mathcal{F}_{K,h}^{\text{ext}}} \qquad (5.13)
$$

for each polytopal cell K ∈ T$_{H}$ .

Remark that if the faces of the simplicial mesh T$_{h}$ do not subdivide the faces of the polytopal mesh T$_{H}$ , then F$_{H,h}$ = F$_{H}$ and consequently U$_{σ}$ ′ = U$_{σ}$ for all σ ∈ F$_{K}$ and U ext K = { U$_{σ}$ }$_{σ}$$_{∈F}$$_{K}$ . Figures 7 and 9 give an illustration in two space dimensions.
## 5.6 Fictitious flux reconstruction by lifting the face normal fluxes

In Definition 4.14, on simplicial meshes, we have reconstructed a discrete $H(\text{div}, \Omega)$-conforming flux $u_h$ by lifting the face normal fluxes of a finite volume scheme. We extend here this procedure to general polytopal meshes. In contrast to Theorem 4.18, however, Theorem 5.18 below will not need $u_h$ to evaluate the posteriori error estimate, so one does not need to perform physically/on the computer this reconstruction. Congruently, we may not have a closed formula for the reconstruction, which on top may not be discrete. This is the reason for our naming fictitious.

### 5.6.1 A (fictitious) flux reconstruction under Assumption 5.1

We first treat the generic case of Assumption 5.1; the case of Assumption 5.2 is treated below. Following the concept of lifting operators used in, e.g., Eymard et al. [127, Section 1.2], Kuznetsov and Repin [172], Brezzi et al. [64, Theorem 5.1], Beirão da Veiga [42, Section 2.1], Beirão da Veiga and Manzini [43, Section 2.4], Kuznetsov [173], Vohralík [249, Section 3.2], or Sboui et al. [230], we will extend the fluxes $U_\sigma$, $\sigma \in \mathcal{F}_K$, to the interior of each polytopal element $K \in \mathcal{T}_H$. The obtained flux $u_h|_K$ will be the approximation of the exact flux $u$ from (5.3) for which we will estimate the energy error $\|u - u_h\|_{\mathcal{K}^{-\frac{1}{2}}}$. Let us first present a motivation.

**Motivation 5.4** (Motivation for the fictitious flux reconstruction). Following [249, equation (3.14)], consider in each element $K \in \mathcal{T}_H$ the infinite-dimensional Darcy problem: for the fluxes $\{U_\sigma\}_{\sigma \in \mathcal{F}_K}$ and the potential $P_K$ from Assumption 5.1, find $\tilde{p}_K: K \to \mathbb{R}$ such that

$$ -\nabla \cdot (\mathbf{K} \nabla \tilde{p}_K) = f|_K, \qquad (5.14a) $$

$$ \frac{(\tilde{p}_K, 1)_K}{|K|} = P_K, \qquad (5.14b) $$

$$ -\mathbf{K} \nabla \tilde{p}_K \cdot \mathbf{n}_\sigma = \frac{U_\sigma}{|\sigma|} \quad \forall \sigma \in \mathcal{F}_K. \qquad (5.14c) $$

We have actually already seen this problem, in Definition 4.16 of Section 4.9. There, on a simplicial mesh, the solution was easily practically computable and discrete: the second-order polynomial denoted as $\tilde{p}_h|_K$ (indeed, (4.24) gives a solution to (5.14)) on a simplicial mesh. Unfortunately, for a general polytopal mesh element $K \in \mathcal{T}_H$, (5.14) amounts to a partial differential equation on $K$, with the solution $\tilde{p}_K$ living in the infinite-dimensional Sobolev space $H^1(K)$ of Definition 3.3 and the associated flux $\tilde{u}_K := -\mathbf{K} \nabla \tilde{p}_K$ living in the infinite-dimensional Sobolev space $H(\text{div}, K)$ of Definition 3.7.

Based on Motivation 5.4, we would like to produce a discrete version of (5.14), which will in particular enable us to derive (easily) computable a posteriori error estimates. Following [249, Theorem 5.1] and [254, Theorems 7.2 and 7.3], we consider the mixed finite element method on the simplicial submesh $\mathcal{T}_K$ of $K$ to approximate the solution $\tilde{p}_K$ of (5.14). Note that this is to be done separately for each $K \in \mathcal{T}_H$, so it is an intrinsically parallel procedure. Recall Definition 5.3 and consider the following sets and spaces of piecewise polynomials over the simplicial submesh $\mathcal{T}_K$ of $K$:

$$ \begin{align} V_{h,N}^K &:= \{v_h \in \mathcal{R}\mathcal{T}_0(\mathcal{T}_K) \cap \mathbf{H}(\text{div}, K); \langle v_h \cdot \mathbf{n}_\sigma, 1 \rangle_\sigma = U_\sigma \quad \forall \sigma \in \mathcal{F}_{K,h}^{\text{ext}}\}, \nonumber \\ V_{h,0}^K &:= \{v_h \in \mathcal{R}\mathcal{T}_0(\mathcal{T}_K) \cap \mathbf{H}(\text{div}, K); \langle v_h \cdot \mathbf{n}_\sigma, 1 \rangle_\sigma = 0 \quad \forall \sigma \in \mathcal{F}_{K,h}^{\text{ext}}\}, \tag{5.15} \\ Q_{h,N}^K &:= \{q_h \in \mathcal{P}_0(\mathcal{T}_K); \frac{(q_h, 1)_K}{|K|} = P_K\}, \nonumber \\ Q_{h,0}^K &:= \{q_h \in \mathcal{P}_0(\mathcal{T}_K); (q_h, 1)_K = 0\}. \nonumber \end{align} $$

Here, $V_{h,N}^K$ and $V_{h,0}^K$ are subsets/subspaces of the lowest-order Raviart-Thomas spaces as defined in Section 3.5.3, whereas $Q_{h,N}^K$ and $Q_{h,0}^K$ are subsets/subspaces of the piecewise polynomial spaces as defined in Section 3.5.1. We then let:

**Definition 5.5** (Flux reconstruction on polytopal meshes under Assumption 5.1). Let $f$ be constant on each $K \in \mathcal{T}_H$ and let a discretization scheme of the form of Assumption 5.1 be given, leading to the normal face fluxes $U_K^{\text{ext}} = \{U_\sigma\}_{\sigma \in \mathcal{F}_K}$, $K \in \mathcal{T}_H$, and consequently to $\{U_\sigma\}_{\sigma \in \mathcal{F}_{K,h}^{\text{ext}}}$ as per Definition 5.3. For all $K \in \mathcal{T}_H$, define $u_h|_K$ by the mixed finite element discretization of (5.14), described by the constrained quadratic minimization problem

$$ u_h|_K := \arg \min_{v_h \in V_{h,N}^K, \nabla \cdot v_h = \text{constant} = f|_K} \|v_h\|_{\mathbf{K}^{-\frac{1}{2}}, K}^2. \qquad (5.16) $$
<figure><img src="image_18.png" /><figcaption>Figure 17: Face barycentres $\mathbf{x}_{K,\sigma}$ and vertices $\mathbf{a}_{K,\sigma}$ opposite to a face $\sigma \in \mathcal{F}_K$. Simplex $K \in \mathcal{T}_h$ in space dimensions $d=2$ (left) and $d=3$ (right)</figcaption></figure>

Equivalently, writing the Euler–Lagrange conditions of (5.16) and imposing the divergence constraint with a Lagrange multiplier, $\mathbf{u}_h|_K \in \mathbf{V}_{h,N}^K$ together with $p_h|_K \in Q_{h,0}^K$ are given by

$$ (\underline{\mathbf{K}}^{-1}\mathbf{u}_h, \mathbf{v}_h)_K - (p_h, \nabla \cdot \mathbf{v}_h)_K = 0 \qquad \forall \mathbf{v}_h \in \mathbf{V}_{h,0}^K, \qquad (5.17a) $$

$$ -(\nabla \cdot \mathbf{u}_h, q_h)_K = -(f, q_h)_K = 0 \qquad \forall q_h \in Q_{h,0}^K. \qquad (5.17b) $$

It is worth noting that the divergence constraint in (5.16) can be written with either a generic constant, or more specifically with $f|_K$; indeed, the assumption that $f$ is constant on each $K \in \mathcal{T}_H$, the definition (5.15) of $\mathbf{V}_{h,N}^K$, and the flux balance (5.8) give the equivalence. Congruently, (5.17b) can be equivalently written as

$$ (\nabla \cdot \mathbf{u}_h, q_h)_K = (f, q_h)_K \qquad \forall q_h \in \mathcal{P}_0(\mathcal{T}_K), $$

which is further equivalent to $\nabla \cdot \mathbf{u}_h|_K = f|_K$. Note also that since the normal fluxes $U_\sigma$ are univalued on the faces from $\mathcal{F}_{H,h}$, the resulting flux reconstruction $\mathbf{u}_h$ has the normal trace continuous over all interior mesh faces from $\mathcal{F}_{H,h}$ and thus belongs to the Raviart–Thomas space $\mathcal{RT}_0(\mathcal{T}_h) \cap \mathbf{H}(\mathrm{div}, \Omega)$ associated with the simplicial mesh $\mathcal{T}_h$ of the entire domain $\Omega$. In consequence, we have:

**Lemma 5.6** (Flux reconstruction on polytopal meshes under Assumption 5.1). Let $f$ be constant on each $K \in \mathcal{T}_H$. Then the flux reconstruction of Definition 5.5 satisfies

$$ \mathbf{u}_h \in \mathcal{R}\mathcal{T}_0(\mathcal{T}_h) \cap \mathbf{H}(\mathrm{div}, \Omega) \quad \text{with} \quad \nabla \cdot \mathbf{u}_h = f. $$

**Construction 5.7** (Matrix form of problem (5.17)). Let $\mathbf{v}_\sigma$, $\sigma \in \mathcal{F}_{K,h}^{\mathrm{int}}$, be the basis functions of the space $\mathbf{V}_{h,0}^K$, and $\mathbf{v}_\sigma$, $\sigma \in \mathcal{F}_{K,h}^{\mathrm{ext}}$, be the remaining basis functions of the space $\mathcal{R}\mathcal{T}_0(\mathcal{T}_K) \cap \mathbf{H}(\mathrm{div}, K)$, see (3.13) and Figures 6 and 17. Similarly, let $q_i$, $1 \le i \le |Q_{h,0}^K|$, be the basis functions of the space $Q_{h,0}^K$ (piecewise constants with mean value 0). Consider these basis functions of $\mathbf{V}_{h,0}^K$ and $Q_{h,0}^K$ as test functions $\mathbf{v}_h$ and $q_h$ in (5.17) and develop

$$ \mathbf{u}_h|_K = \sum_{\substack{\sigma \in \mathcal{F}_{K,h}^{\mathrm{ext}} \\ \mathbf{u}_{h,K}^{\mathrm{ext}}}} (U_K^{\mathrm{ext}})_\sigma \mathbf{v}_\sigma + \sum_{\substack{\sigma \in \mathcal{F}_{K,h}^{\mathrm{int}} \\ \mathbf{u}_{h,K}^{\mathrm{int}}}} (U_K^{\mathrm{int}})_\sigma \mathbf{v}_\sigma $$

and

$$ p_h|_K = P_K + \sum_{\substack{i=1 \\ p_{h,K}^0}}^{|Q_{h,0}^K|} (P_K^0)_i q_i, $$

where $\{U_\sigma\}_{\sigma \in \mathcal{F}_K}$ and $P_K$ are the given data of the problem (5.14) and the data $U_K^{\mathrm{ext}} := \{U_\sigma\}_{\sigma \in \mathcal{F}_{K,h}^{\mathrm{ext}}}$ appearing in problem (5.16) are given by Definition 5.3. Then (5.17) in matrix form corresponds to: find $U_K^{\mathrm{int}} \in \mathbb{R}^{|\mathbf{V}_{h,0}^K|}$ and $P_K^0 \in \mathbb{R}^{|Q_{h,0}^K|}$ such that

$$ \left( \begin{array}{cc} \mathbf{A}_K^{\mathrm{int},\mathrm{int}} & (\mathbb{B}_K^{0,\mathrm{int}})^{\mathrm{t}} \\ \mathbb{B}_K^{0,\mathrm{int}} & 0 \end{array} \right) \left( \begin{array}{c} \mathbf{U}_K^{\mathrm{int}} \\ \mathbf{P}_K^0 \end{array} \right) = \left( \begin{array}{c} -\mathbf{A}_K^{\mathrm{int},\mathrm{ext}} \mathbf{U}_K^{\mathrm{ext}} \\ -\mathbb{B}_K^{0,\mathrm{ext}} \mathbf{U}_K^{\mathrm{ext}} \end{array} \right), \qquad (5.18) $$
where

$$A_K^{\text{int,int}} \in \mathbb{R}^{|\mathcal{F}_{K,h}^{\text{int}}| \times |\mathcal{F}_{K,h}^{\text{ext}}|}, \quad (A_K^{\text{int,int}})_{\sigma, \sigma'} := (\underline{\mathbf{K}}^{-1} \mathbf{v}_{\sigma'}, \mathbf{v}_{\sigma})_K, \quad \sigma, \sigma' \in \mathcal{F}_{K,h}^{\text{int}}, \quad (5.19a)$$

$$A_K^{\text{int,ext}} \in \mathbb{R}^{|\mathcal{F}_{K,h}^{\text{int}}| \times |\mathcal{F}_{K,h}^{\text{ext}}|}, \quad (A_K^{\text{int,ext}})_{\sigma, \sigma'} := (\underline{\mathbf{K}}^{-1} \mathbf{v}_{\sigma'}, \mathbf{v}_{\sigma})_K, \quad \sigma \in \mathcal{F}_{K,h}^{\text{int}}, \sigma' \in \mathcal{F}_{K,h}^{\text{ext}}, \quad (5.19b)$$

$$\mathbb{B}_K^{0,\text{int}} \in \mathbb{R}^{|\mathcal{F}_{h,0}^K| \times |\mathcal{F}_{K,h}^{\text{int}}|}, \quad (\mathbb{B}_K^{0,\text{int}})_{i,\sigma'} := -(\nabla \cdot \mathbf{v}_{\sigma'}, q_i)_K, \quad 1 \le i \le |\mathcal{F}_{h,0}^K|, \sigma' \in \mathcal{F}_{K,h}^{\text{int}}, \quad (5.19c)$$

$$\mathbb{B}_K^{0,\text{ext}} \in \mathbb{R}^{|\mathcal{F}_{h,0}^K| \times |\mathcal{F}_{K,h}^{\text{ext}}|}, \quad (\mathbb{B}_K^{0,\text{ext}})_{i,\sigma'} := -(\nabla \cdot \mathbf{v}_{\sigma'}, q_i)_K, \quad 1 \le i \le |\mathcal{F}_{h,0}^K|, \sigma' \in \mathcal{F}_{K,h}^{\text{ext}}, \quad (5.19d)$$

It will also be useful to define

$$(A_K^{\text{ext,ext}})_{\sigma, \sigma'} := (\underline{\mathbf{K}}^{-1} \mathbf{v}_{\sigma'}, \mathbf{v}_{\sigma})_K, \quad \sigma, \sigma' \in \mathcal{F}_{K,h}^{\text{ext}}.$$

Definition 5.5 is practical and we can solve problem (5.17) or (5.18) to have $\mathbf{u}_h$ available. We, however, want to avoid this (unless we want to compute the energy error $\|\mathbf{u} - \mathbf{u}_h\|_{\underline{\mathbf{K}}^{-\frac{1}{2}},K}$ in model cases with known exact solution $\mathbf{u}$). In Theorem 5.18 below, we will principally only need the energy $\|\mathbf{u}_h\|_{\underline{\mathbf{K}}^{-\frac{1}{2}},K}$. In this respect, we observe that Definition 5.5 lifts the information from the boundary of each polytopal element $K \in \mathcal{T}_H$ given by the fluxes $\mathcal{U}_K^{\text{ext}}$ to the interior of the element $K$. It is thus clear that $\|\mathbf{u}_h\|_{\underline{\mathbf{K}}^{-\frac{1}{2}},K}$ must only depend on $\mathcal{U}_K^{\text{ext}}$ and can be obtained directly therefrom. It turns out that this link can be expressed by a single element matrix constructed solely from the geometry of the simplicial submesh $\mathcal{T}_K$ of $K$, using Construction 5.7. As in [254, proof of Theorem 7.3], we indeed have:

**Lemma 5.8 (Energy norm and MFE element matrix).** Let a discretization scheme of the form of Assumption 5.1 be given, leading to the normal face fluxes $\mathcal{U}_K^{\text{ext}} = \{\mathcal{U}_{\sigma}\}_{\sigma \in \mathcal{F}_K}$, $K \in \mathcal{T}_H$, and consequently to $\{\mathcal{U}_{\sigma}\}_{\sigma \in \mathcal{F}_{K,h}^{\text{ext}}}$ as per Definition 5.3. For each $K \in \mathcal{T}_H$, let $\mathbf{u}_h|_K \in \mathbf{V}_{h,N}^K$ be given by Definition 5.5. Relying on Construction 5.7, define the element matrix of the mixed finite element method

$$\hat{\mathbf{A}}_{\text{MFE},K} := A_K^{\text{ext,ext}} - \begin{pmatrix} A_K^{\text{int,ext}} \\ \mathbb{B}_K^{0,\text{ext}} \end{pmatrix}^t \begin{pmatrix} A_K^{\text{int,int}} & (\mathbb{B}_K^{0,\text{int}})^t \\ \mathbb{B}_K^{0,\text{int}} & 0 \end{pmatrix}^{-1} \begin{pmatrix} A_K^{\text{int,ext}} \\ \mathbb{B}_K^{0,\text{ext}} \end{pmatrix}, \quad (5.20)$$

where the building matrices are given by (5.19). Then

$$\|\mathbf{u}_h\|_{\underline{\mathbf{K}}^{-\frac{1}{2},K}}^2 = (\mathcal{U}_K^{\text{ext}})^t \hat{\mathbf{A}}_{\text{MFE},K} \mathcal{U}_K^{\text{ext}}. \quad (5.21)$$

*Proof.* Let $K \in \mathcal{T}_H$ be a fixed polytopal mesh element. Following Construction 5.7, decompose $\mathbf{u}_h|_K = \mathbf{u}_{h,K}^{\text{ext}} + \mathbf{u}_{h,K}^{\text{int}}$ and $p_h|_K = P_K + p_{h,K}^0$. Note that choosing $\mathbf{u}_{h,K}^{\text{int}} \in \mathbf{V}_{h,0}^K$ as the test function in (5.17a) and $p_{h,K}^0 \in Q_{h,0}^K$ as the test function in (5.17b), one has

$$(\underline{\mathbf{K}}^{-1} \mathbf{u}_h, \mathbf{u}_{h,K}^{\text{int}})_K \stackrel{(5.17a)}{=} (p_h, \nabla \cdot \mathbf{u}_{h,K}^{\text{int}})_K = (p_{h,K}^0, \nabla \cdot \mathbf{u}_{h,K}^{\text{int}})_K \stackrel{(5.17b)}{=} -(p_{h,K}^0, \nabla \cdot \mathbf{u}_{h,K}^{\text{ext}})_K,$$

where we have also used $(P_K, \nabla \cdot \mathbf{u}_{h,K}^{\text{int}})_K = 0$ in the middle equality. Consequently,

$$\begin{align*}
\|\mathbf{u}_h\|_{\underline{\mathbf{K}}^{-\frac{1}{2},K}}^2 &= (\underline{\mathbf{K}}^{-1} \mathbf{u}_h, \mathbf{u}_h)_K \\
&= (\underline{\mathbf{K}}^{-1} \mathbf{u}_h, \mathbf{u}_{h,K}^{\text{int}})_K + (\underline{\mathbf{K}}^{-1} \mathbf{u}_{h,K}^{\text{int}}, \mathbf{u}_{h,K}^{\text{ext}})_K + (\underline{\mathbf{K}}^{-1} \mathbf{u}_{h,K}^{\text{ext}}, \mathbf{u}_{h,K}^{\text{ext}})_K \\
&= -(p_{h,K}^0, \nabla \cdot \mathbf{u}_{h,K}^{\text{ext}})_K + (\underline{\mathbf{K}}^{-1} \mathbf{u}_{h,K}^{\text{int}}, \mathbf{u}_{h,K}^{\text{ext}})_K + (\underline{\mathbf{K}}^{-1} \mathbf{u}_{h,K}^{\text{ext}}, \mathbf{u}_{h,K}^{\text{ext}})_K \\
&= (\mathcal{U}_K^{\text{ext}})^t \begin{pmatrix} A_K^{\text{int,ext}} \\ \mathbb{B}_K^{0,\text{ext}} \end{pmatrix}^t \begin{pmatrix} \mathcal{U}_K^{\text{int}} \\ P_K^0 \end{pmatrix} + (\mathcal{U}_K^{\text{ext}})^t \hat{\mathbf{A}}_{K}^{\text{ext,ext}} \mathcal{U}_K^{\text{ext}}.
\end{align*}$$

Combining this with the MFE matrix form (5.18) and the MFE element matrix definition (5.20) yields the claim (5.21). $\square$

### 5.6.2 A fictitious flux reconstruction under Assumption 5.2

There exist other ways how to lift the face normal fluxes $\mathcal{U}_K^{\text{ext}}$ into the interior of each polytope $K \in \mathcal{T}_H$ than via Definition 5.5. For instance, one could use a finite volume discretization of (5.14) (like (4.19a)) in place of the mixed finite element one (5.17). Under Assumption 5.2, which defines $\hat{\mathbf{u}}_h$, a much more direct and simple procedure suggests itself: to use directly the lifting operator from (5.11). This is in particular the theoretical vehicle in mimetic finite differences according to, e.g., [64, Theorem 5.1]. In
contrast to (5.16), where we can solve for $\tilde{u}_h$ if needed, the reconstruction $\tilde{u}_h$ we now introduce may be truly fictitious in the sense that it may not be practically possible to construct it. Thus, we will not be able to compute the actual error $\|u - \tilde{u}_h\|_{\mathcal{K}^{-\frac{1}{2}}}$ in model cases with known exact flux $u$, but otherwise this not be an issue (we will still be able to compute the a posteriori error estimate of Theorem 5.18 below). In particular, thanks to (5.11c), we have the following trivial counterpart of Lemma 5.8:

**Lemma 5.9** (Energy norm and given discretization scheme element matrix). Let a discretization scheme of the form of Assumption 5.2 be given, leading to the element matrices $\hat{\mathbb{A}}_K$ and normal face fluxes $U_K^{\text{ext}} = \{U_\sigma\}_{\sigma \in \mathcal{F}_K}$, $K \in \mathcal{T}_H$. Then the lifting $\tilde{u}_h$ of (5.11) satisfies

$$ \|\tilde{u}_h\|_{\mathcal{K}^{-\frac{1}{2},K}}^2 = (U_K^{\text{ext}})^t \hat{\mathbb{A}}_K U_K^{\text{ext}}. \qquad (5.22) $$

Following [254, Theorem 7.3], the mixed finite element element matrix $\hat{\mathbb{A}}_{MFE,K}$ defined by (5.20) belongs to the same family as the element matrices $\hat{\mathbb{A}}_K$ of mimetic finite differences, mixed finite volumes, and hybrid finite volumes of Assumption 5.2 when the simplicial faces of $\mathcal{T}_h$ do not subdivide the faces of the polytopal mesh $\mathcal{T}_H$. This leads to the following strong interconnection of Lemmas 5.8 and 5.9 (this is an illuminating observation not needed in practice/thory):

**Lemma 5.10** (Approximate energy norm). Let the assumptions of Lemmas 5.8 and 5.9 and of [64, Theorem 5.1] be satisfied. Let the simplicial faces of $\mathcal{T}_h$ do not subdivide the faces of the polytopal mesh $\mathcal{T}_H$. Then there exist two constants $c, C > 0$ independent of the matrices $\hat{\mathbb{A}}_K, \hat{\mathbb{A}}_{MFE,K}$, and the mesh size but possibly depending on the shape-regularity of the simplicial submesh $\mathcal{T}_K$ such that

$$ c \|\tilde{u}_h\|_{\mathcal{K}^{-\frac{1}{2},K}}^2 \le \underbrace{\|u_h\|_{\mathcal{K}^{-\frac{1}{2},K}}^2}_{(U_K^{\text{ext}})^t \hat{\mathbb{A}}_K U_K^{\text{ext}}} \le C \underbrace{\|\tilde{u}_h\|_{\mathcal{K}^{-\frac{1}{2},K}}^2}_{(U_K^{\text{ext}})^t \hat{\mathbb{A}}_{MFE,K} U_K^{\text{ext}}} \qquad (5.23) $$

that we denote by

$$ \|u_h\|_{\mathcal{K}^{-\frac{1}{2},K}}^2 \approx (U_K^{\text{ext}})^t \hat{\mathbb{A}}_K U_K^{\text{ext}}. \qquad (5.24) $$

*Proof.* Any scheme of the form of Assumption 5.2 with the lifting operator according to [64, Theorem 5.1] satisfies the equality (5.7) from [64], i.e., (5.11c). As mixed finite elements on polytopal meshes belong to the mimetic finite difference family by [254, Theorem 7.3], this implies (5.23). $\square$

## 5.7 Potential point values

In Theorem 5.18 below, we will need a vector of real values with one value $Z_a$ per vertex $a \in V_h$ of the simplicial submesh $\mathcal{T}_h$ of $\mathcal{T}_H$ and one real value $Z_\sigma$ per face $\sigma \in \mathcal{F}_{H,h}$ of the simplicial submesh $\mathcal{T}_h$ of $\mathcal{T}_H$ (approximate point values of the pressure heads). We identify them here for the two Assumptions 5.1 and 5.2.

### 5.7.1 Potential point values under Assumption 5.1

Here, under the general Assumption 5.1, we define the vectors $\{Z_a\}_{a \in V_h}$ and $\{Z_\sigma\}_{\sigma \in \mathcal{F}_{H,h}}$ by average values of $P_K$ in the neighboring polytopal elements:

**Definition 5.11** (Potential point values on polytopal meshes under Assumption 5.1). Let a discretization scheme of the form of Assumption 5.1 be given, leading to the elementwise pressure heads $P := \{P_K\}_{K \in \mathcal{T}_H}$. For each vertex $a$ of the simplicial mesh $\mathcal{T}_h$ lying on $\partial K$ for some $K \in \mathcal{T}_H$ but not on $\partial \Omega$, let $\mathcal{T}_a$ denote the set of polytopal elements $K \in \mathcal{T}_H$ sharing $a$. We set

$$ Z_a := \frac{1}{|\mathcal{T}_a|} \sum_{K \in \mathcal{T}_a} P_K, \qquad (5.25a) $$

where, recall, $|\mathcal{T}_a|$ is the cardinality of the set $\mathcal{T}_a$, i.e., the number of polytopal elements $K \in \mathcal{T}_H$ sharing $a$. We also set

$$ Z_a := P_K \text{ for any vertex } a \text{ of } \mathcal{T}_h \text{ lying inside some } K \in \mathcal{T}_H, \qquad (5.25b) $$

$$ Z_a := 0 \text{ for any vertex } a \text{ lying simultaneously on } \partial K \text{ for some } K \in \mathcal{T}_H \text{ and on } \partial \Omega, \qquad (5.25c) $$
<figure><img src="image_19.png" /><figcaption>Figure 18: Example of element vectors $Z_K$ and $Z_K^{\text{ext}}$ of nodal and facial potential point values</figcaption></figure>

where the latter condition is set in accordance with the boundary condition (5.1b). We will also need the face values $\{Z_\sigma\}_{\sigma \in \mathcal{F}_{H,h}}$

$$Z_{\sigma} := \frac{1}{d} \sum_{\mathbf{a} \in \mathcal{V}_{\sigma}} Z_{\mathbf{a}} \qquad (5.26)$$

for any simplicial face $\sigma \in \mathcal{F}_{H,h}$, where, recall, $\mathcal{V}_\sigma$ collects the vertices of the given face $\sigma \in \mathcal{F}_{H,h}$ and $d$ is the space dimension. Then, for each polytopal cell $K \in \mathcal{T}_H$, we let the element vector $Z_K \in \mathbb{R}^{|\mathcal{V}_{K,h}|}$ collect the values $Z_\mathbf{a}$ for the vertices $\mathbf{a}$ of the simplicial mesh $\mathcal{T}_K$ lying on the boundary $\partial K$ and the value(s) inside $K$

$$Z_K := \{Z_\mathbf{a}\}_{\mathbf{a} \in \mathcal{V}_{K,h}} \qquad (5.27)$$

and the element vector $Z_K^{\text{ext}} \in \mathbb{R}^{|\mathcal{F}_{K,h}^{\text{ext}}|}$ collect the values associated with faces $\sigma \in \mathcal{F}_{K,h}^{\text{ext}}$

$$Z_K^{\text{ext}} := \{Z_\sigma\}_{\sigma \in \mathcal{F}_{K,h}^{\text{ext}}}. \qquad (5.28)$$

Figure 18 gives an illustration in two space dimensions.

### 5.7.2 Potential point values under Assumption 5.2

A more straightforward and precise potential reconstruction can be designed under Assumption 5.2. The main idea is that any scheme of the form (5.9) can be hybridized, following, e.g., Roberts and Thomas [224], Brezzi and Fortin [63], or Droniou et al. [110], giving rise to one Lagrange multiplier $\Lambda_\sigma$ per face $\sigma \in \mathcal{F}_H$:

**Definition 5.12** (Lagrange multipliers under Assumption 5.2). Let a discretization scheme of the form of Assumption 5.2 be given and fix a face $\sigma \in \mathcal{F}_H$. Consider the line associated with $\sigma$ in the first block equation of (5.9). For an interior face $\sigma$, let it be shared by two polytopes $K, L \in \mathcal{T}_H$ such that $\mathbf{n}_\sigma$ points from $K$ to $L$. From the structure of the matrices $\mathbb{A}$ and $\mathbb{B}$ supposed in Assumption 5.2, it follows that

$$\sum_{\sigma' \in \mathcal{F}_K} \mathbf{n}_K \cdot \mathbf{n}_{\sigma'} (\hat{\mathbb{A}}_K)_{\sigma, \sigma'} \mathbf{U}_{\sigma'} - \mathbf{P}_K = \sum_{\sigma' \in \mathcal{F}_L} \mathbf{n}_{L, \sigma'} \cdot \mathbf{n}_{\sigma'} (\hat{\mathbb{A}}_L)_{\sigma, \sigma'} \mathbf{U}_{\sigma'} - \mathbf{P}_L,$$

so that these expressions are univ alued from both elements $K$ and $L$. Such an expression is also clearly univ alued on boundary faces. It allows to define the Lagrange multipliers

$$\Lambda_{\sigma} := \mathbf{P}_{K} - \sum_{\sigma' \in \mathcal{F}_{K}} \mathbf{n}_{K} \cdot \mathbf{n}_{\sigma'} (\hat{\mathbb{A}}_{K})_{\sigma, \sigma'} \mathbf{U}_{\sigma'} \qquad (5.29)$$

for each element $K \in \mathcal{T}_H$ and each face $\sigma \in \mathcal{F}_K$.

**Definition 5.13** (Potential point values on polytopal meshes under Assumption 5.2). Let a discretization scheme of the form of Assumption 5.2 be given, leading to the normal face fluxes $\mathbf{U} = \{\mathbf{U}_\sigma\}_{\sigma \in \mathcal{F}_H}$ and the elementwise pressure heads $\mathbf{P} := \{\mathbf{P}_K\}_{K \in \mathcal{T}_H}$. Let the Lagrange multipliers $\Lambda_\sigma$ be given by Definition 5.12. For each vertex $\mathbf{a}$ of the simplicial mesh $\mathcal{T}_h$ lying on $\partial K$ for some $K \in \mathcal{T}_H$ but not on $\partial \Omega$, let $\mathcal{F}_\mathbf{a}$ denote the set of polytopal faces $\sigma \in \mathcal{F}_H$ sharing $\mathbf{a}$ and $|\mathcal{F}_\mathbf{a}|$ its cardinality. We set

$$Z_{\mathbf{a}} := \frac{1}{|\mathcal{F}_{\mathbf{a}}|} \sum_{\sigma \in \mathcal{F}_{\mathbf{a}}} \Lambda_{\sigma}. \qquad (5.30a)$$
We also set

$$Z_a := P_K \text{ for any vertex } \mathbf{a} \text{ of } \mathcal{T}_h \text{ lying inside some } K \in \mathcal{T}_H, \qquad (5.30b)$$

$$Z_a := 0 \text{ for any vertex } \mathbf{a} \text{ lying simultaneously on } \partial K \text{ for some } K \in \mathcal{T}_H \text{ and on } \partial \Omega, \qquad (5.30c)$$

where the latter condition is set in accordance with the boundary condition (5.1b). We then define the face values $\{Z_\sigma\}_{\sigma \in \mathcal{F}_{H,h}}$ and the element vectors $Z_K \in \mathbb{R}^{|\mathcal{V}_{K,h}|}$ and $Z_K^{\text{ext}} \in \mathbb{R}^{|\mathcal{F}_{K,h}^{\text{ext}}|}$ by respectively (5.26), (5.27), and (5.28) as in Definition 5.11.

## 5.8 Fictitious potential reconstruction by respecting the point values

In Definition 4.17, we have reconstructed a discrete $H_0^1(\Omega)$-conforming potential $\zeta_h$ from a finite volume approximation on a simplicial mesh. We extend here this procedure to general polytopal meshes. As in Section 5.6, in contrast to Theorem 4.18, however, Theorem 5.18 below will not need $\zeta_h$ to evaluate the posteriori error estimate, so one again does not need to perform physically/on the computer this reconstruction. This is the reason for our naming fictitious (though we will give a closed formula for $\zeta_h$ and $\zeta_h$ will be discrete and possible to construct if desired). Recall the notation from Section 3.6.

Given the point values $Z_a$, either prescribed by (5.25) from Definition 5.11 or by (5.30) of Definition 5.13, the potential reconstruction $\zeta_h \in P_1(\mathcal{T}_h) \cap H_0^1(\Omega)$ (piecewise affine with respect to the simplicial submesh $\mathcal{T}_h$ and $H_0^1(\Omega)$-conforming) is trivial:

**Definition 5.14** (Potential reconstruction on polytopal meshes). Let the potential point values $Z_a$ be given by either Definition 5.11 or Definition 5.13. Then the potential reconstruction $\zeta_h \in P_1(\mathcal{T}_h) \cap H_0^1(\Omega)$ is prescribed by

$$\zeta_h(\mathbf{a}) := Z_a \quad \forall \mathbf{a} \in \mathcal{V}_h.$$

Note that by (5.26), since $\zeta_h$ is affine on each face,

$$(Z_K^{\text{ext}})_\sigma = Z_\sigma = \zeta_h(\mathbf{x}_\sigma), \qquad (5.31)$$

i.e., the face values $Z_\sigma$ are the punctual values of the reconstruction $\zeta_h$ in the face barycenters $\mathbf{x}_\sigma$, $\sigma \in \mathcal{F}_{K,h}^{\text{ext}}$, cf. Figure 17.

Consider the usual hat basis functions $\psi_\mathbf{a}$ of the space $P_1(\mathcal{T}_h) \cap H^1(\Omega)$ on the simplicial mesh $\mathcal{T}_h$: $\psi_\mathbf{a}$ is piecewise affine on $\mathcal{T}_h$ and $H^1(\Omega)$-conforming, $\psi_\mathbf{a}(\mathbf{a}') = 1$ if $\mathbf{a} = \mathbf{a}'$ and 0 otherwise, where $\mathbf{a}, \mathbf{a}' \in \mathcal{V}_h$ are the vertices of $\mathcal{T}_h$. For any polytopal element $K \in \mathcal{T}_H$, we will need below the element stiffness matrix $\hat{S}_{\text{FE},K} \in \mathbb{R}^{|\mathcal{V}_{K,h}| \times |\mathcal{V}_{K,h}|}$ defined by

$$(\hat{S}_{\text{FE},K})_{\mathbf{a},\mathbf{a}'} := (\underline{\mathbf{K}}\nabla\psi_{\mathbf{a}'}, \nabla\psi_{\mathbf{a}})_K \quad \mathbf{a}, \mathbf{a}' \in \mathcal{V}_{K,h} \qquad (5.32)$$

and the element mass matrix $\hat{M}_{\text{FE},K} \in \mathbb{R}^{|\mathcal{V}_{K,h}| \times |\mathcal{V}_{K,h}|}$ defined by

$$(\hat{M}_{\text{FE},K})_{\mathbf{a},\mathbf{a}'} := (\psi_{\mathbf{a}'}, \psi_{\mathbf{a}})_K \quad \mathbf{a}, \mathbf{a}' \in \mathcal{V}_{K,h}. \qquad (5.33)$$

An immediate consequence is:

**Lemma 5.15** (Energy norm and integration using the stiffness and mass matrices). Let the potential reconstruction be given by Definition 5.14 and the stiffness and mass matrices respectively by (5.32) and (5.33). Then there holds

$$\| \underline{\mathbf{K}} \nabla \zeta_h \|^2_{\underline{\mathbf{K}}^{-\frac{1}{2}}, K} = Z_K^t \hat{S}_{\text{FE},K} Z_K \qquad (5.34)$$

and

$$(1, \zeta_h)_K = 1^t \hat{M}_{\text{FE},K} Z_K. \qquad (5.35)$$

**Remark 5.16** (Element matrices $\hat{S}_{\text{FE},K}$ and $\hat{M}_{\text{FE},K}$). To obtain the finite element stiffness and mass matrices $\hat{S}_{\text{FE},K}$ of (5.32) and $\hat{M}_{\text{FE},K}$ of (5.33), one could use a finite element assembly code linked to a simplicial mesh on the polytopal cell $K \in \mathcal{T}_H$. Note, however, that there exist simple analytical formulas for the basis functions $\psi_\mathbf{a}$ of piecewise affine Lagrange finite elements and for their gradients $\nabla\psi_\mathbf{a}$ that merely necessitate the position of the vertices of the polygonal element $K$ and of the central point(s). In this sense, we say that neither $\mathcal{T}_h$ nor $\mathcal{T}_K$ need not be constructed in practice in order to obtain $\hat{S}_{\text{FE},K}$ and $\hat{M}_{\text{FE},K}$.
We finish this section by a remark on an alternative potential reconstruction in the spirit of Definitions 4.16–4.17. This is still more precise, though slightly more computationally demanding; we will use it for comparison in numerical experiments below.

**Remark 5.17** (Piecewise quadratic potential reconstruction). Let $(u_h|_K, p_h|_K)$ be given by Definition 5.5 on all $K \in \mathcal{T}_H$. We let $\tilde{p}_h$ be a piecewise quadratic polynomial on the simplicial mesh $\mathcal{T}_h$ given by

$$-\underline{\mathbf{K}}\nabla \tilde{p}_h|_{\kappa} = \mathbf{u}_h|_{\kappa}, \quad \frac{(\tilde{p}_h, 1)_{\kappa}}{|\kappa|} = p_h|_{\kappa} \quad \forall \kappa \in \mathcal{T}_h.$$

Then a piecewise quadratic potential reconstruction $\zeta_h \in \mathcal{P}_2(\mathcal{T}_h) \cap H_0^1(\Omega)$ can be obtained by averaging the values of $\tilde{p}_h$ in all Lagrange nodes of freedom as per Definition 4.17.

## 5.9 A guaranteed a posteriori error estimate with inexpensive implementation and evaluation

We now present our a posteriori error estimate for the model problem (5.1). Our main result is given for the error between the exact Darcy velocity $\mathbf{u} \in \mathbf{H}(\text{div}, \Omega)$ of (5.2)–(5.3) and the reconstruction $\mathbf{u}_h \in \mathcal{RT}_0(\mathcal{T}_h) \cap \mathbf{H}(\text{div}, \Omega)$ of Definition 5.5 evaluated as per (5.5), under the general Assumption 5.1:

**Theorem 5.18** (A guaranteed a posteriori error estimate with inexpensive implementation and evaluation).

Let $f$ be constant on each $K \in \mathcal{T}_H$ and let the exact flux $\mathbf{u}$ be given by (5.2)–(5.3). For any polytopal discretization satisfying Assumption 5.1, let, for each polytopal element $K \in \mathcal{T}_H$, the element vectors of face normal fluxes $U_K^{\text{ext}}$ be given by Definition 5.3. Let the element vectors of potential point values $Z_K$ and $Z_K^{\text{ext}}$ be given by Definition 5.11. Let finally the matrices $\hat{\mathbf{A}}_{\text{MFE},K}$, $\hat{\mathbf{S}}_{\text{FE},K}$, and $\hat{\mathbf{M}}_{\text{FE},K}$ be respectively defined by (5.20), (5.32), and (5.33). Then there holds

$$\|\mathbf{u} - \mathbf{u}_h\|_{\underline{\mathbf{K}}^{-\frac{1}{2}}} \le \left\{ \sum_{K \in \mathcal{T}_H} \eta_K^2 \right\}^{1/2}, \quad (5.36)$$

where

$$\eta_K^2 := (U_K^{\text{ext}})^{\text{t}} \hat{\mathbf{A}}_{\text{MFE},K} U_K^{\text{ext}} + Z_K^{\text{t}} \hat{\mathbf{S}}_{\text{FE},K} Z_K + 2(U_K^{\text{ext}})^{\text{t}} Z_K^{\text{ext}} - 2\mathbf{F}_K |K|^{-1} \mathbf{1}^{\text{t}} \hat{\mathbf{M}}_{\text{FE},K} \bar{Z}_K. \quad (5.37)$$

Here the flux reconstruction $\mathbf{u}_h \in \mathcal{RT}_0(\mathcal{T}_h) \cap \mathbf{H}(\text{div}, \Omega)$ with $\nabla \cdot \mathbf{u}_h = f$ is obtained following Definition 5.5.

*Proof.* Recall from Lemma 5.6 that $\mathbf{u}_h \in \mathbf{H}(\text{div}, \Omega)$ with $\nabla \cdot \mathbf{u}_h = f$. Thus, considering the diffusion tensor $\underline{\mathbf{K}}$ in place of the identity matrix and proceeding as for the Prager–Synge equality of Corollary 4.6, we have

$$\|\mathbf{u} - \mathbf{u}_h\|_{\underline{\mathbf{K}}^{-\frac{1}{2}}} = \min_{v \in H_0^1(\Omega)} \|\mathbf{u}_h + \underline{\mathbf{K}}\nabla v\|_{\underline{\mathbf{K}}^{-\frac{1}{2}}}. \quad (5.38)$$

Consequently, for an arbitrary $\zeta_h \in H_0^1(\Omega)$,

$$\|\mathbf{u} - \mathbf{u}_h\|_{\underline{\mathbf{K}}^{-\frac{1}{2}}} \le \|\mathbf{u}_h + \underline{\mathbf{K}}\nabla \zeta_h\|_{\underline{\mathbf{K}}^{-\frac{1}{2}}}. \quad (5.39)$$

We now choose for $\zeta_h$ the fictitious potential reconstruction $\zeta_h \in \mathcal{P}_1(\mathcal{T}_h) \cap H_0^1(\Omega)$ of Definition 5.14, continuous and piecewise affine with respect to the simplicial submesh $\mathcal{T}_h$ and given by the nodal values of the vector $Z$ from Definition 5.11. We focus on the evaluation of $\|\mathbf{u}_h + \underline{\mathbf{K}}\nabla \zeta_h\|_{\underline{\mathbf{K}}^{-\frac{1}{2}}}$. Developing, we obtain, for each polytopal mesh element $K \in \mathcal{T}_H$,

$$\|\mathbf{u}_h + \underline{\mathbf{K}}\nabla \zeta_h\|_{\underline{\mathbf{K}}^{-\frac{1}{2},K}}^2 = \|\mathbf{u}_h\|_{\underline{\mathbf{K}}^{-\frac{1}{2},K}}^2 + 2(\mathbf{u}_h, \nabla \zeta_h)_K + \|\underline{\mathbf{K}}\nabla \zeta_h\|_{\underline{\mathbf{K}}^{-\frac{1}{2},K}}^2.$$

We now use Lemma 5.8 for the first term and (5.34) from Lemma 5.15 for the last one. For the middle term, recall first that the normal components of vector fields in $\mathcal{RT}_0(\mathcal{T}_h) \cap \mathbf{H}(\text{div}, \Omega)$ are constant on
each face. Thus the Green theorem together with (5.26) and its consequence (5.31), $\nabla \cdot \mathbf{u}_h = f$ and (5.7), and (5.35) from Lemma 5.15 give

$$ (\mathbf{u}_h, \nabla \zeta_h)_K = \langle \mathbf{u}_h \cdot \mathbf{n}, \zeta_h \rangle_{\partial K} - (\nabla \cdot \mathbf{u}_h, \zeta_h)_K = (\mathbf{U}_K^{\text{ext}})^{\text{t}} \mathbf{Z}_K^{\text{ext}} - \mathbf{F}_K |K|^{-1} \mathbf{1}^{\text{t}} \hat{\mathbf{M}}_{\text{FE},K} \mathbf{Z}_K. \quad (5.39) $$

Thus the proof is finished.

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

Note that the a posteriori error estimators $\eta_K$ of Theorem 5.18 take a simple form of matrix-vector multiplication on each polytopal mesh element $K \in \mathcal{T}_H$, only need the element vectors $\mathbf{U}_K^{\text{ext}}$, $\mathbf{Z}_K$, and $\mathbf{Z}_K^{\text{ext}}$ together with the element matrices $\hat{\mathbf{A}}_{\text{MFE},K}$, $\hat{\mathbf{S}}_{\text{FE},K}$, and $\hat{\mathbf{M}}_{\text{FE},K}$, but not the reconstructions $\mathbf{u}_h$ and $\zeta_h$ (which become fictitious), and yet (5.36) delivers a guaranteed upper bound on the Darcy velocity error.

For schemes of the form of Assumption 5.2, Definition 5.13 can alternatively be used in place of Definition 5.11 in Theorem 5.18. In this case, however, it is probably still more attractive to proceed as:

**Corollary 5.19** (A simple guaranteed estimate with the given element matrices $\hat{\mathbf{A}}_K$). Let $f$ be constant on each $K \in \mathcal{T}_H$ and let the exact flux $\mathbf{u}$ be given by (5.2)-(5.3). For any polytopal discretization satisfying Assumption 5.2, let, for each polytopal element $K \in \mathcal{T}_H$, the element vectors of face normal fluxes $\mathbf{U}_K^{\text{ext}}$ be given by Definition 5.3. Let the element vectors of potential point values $\mathbf{Z}_K$ and $\mathbf{Z}_K^{\text{ext}}$ be given by Definition 5.13. Let the matrices $\hat{\mathbf{A}}_K$ be the element matrices from Assumption 5.2. Let finally $\hat{\mathbf{S}}_{\text{FE},K}$ and $\hat{\mathbf{M}}_{\text{FE},K}$ be respectively defined by (5.32) and (5.33). Then there holds

$$ \| \mathbf{u} - \tilde{\mathbf{u}}_h \|_{\underline{\mathbf{K}}^{-\frac{1}{2}}} \le \left\{ \sum_{K \in \mathcal{T}_H} \eta_K^2 \right\}^{1/2}, \quad (5.40) $$

where

$$ \eta_K^2 := (\mathbf{U}_K^{\text{ext}})^{\text{t}} \hat{\mathbf{A}}_K \mathbf{U}_K^{\text{ext}} + \mathbf{Z}_K^{\text{t}} \hat{\mathbf{S}}_{\text{FE},K} \mathbf{Z}_K + 2(\mathbf{U}_K^{\text{ext}})^{\text{t}} \mathbf{Z}_K^{\text{ext}} - 2\mathbf{F}_K |K|^{-1} \mathbf{1}^{\text{t}} \hat{\mathbf{M}}_{\text{FE},K} \mathbf{Z}_K, \quad (5.41) $$

i.e., the same form as (5.37) but with $\hat{\mathbf{A}}_K$ in place of $\hat{\mathbf{A}}_{\text{MFE},K}$. Here the fictitious flux reconstruction $\tilde{\mathbf{u}}_h \in \mathcal{R}\mathcal{T}_0(\mathcal{T}_h) \cap \mathbf{H}(\text{div}, \Omega)$ with $\nabla \cdot \mathbf{u}_h = f$ is obtained following (5.11).

**Proof.** The proof is similar to that of Theorem 5.18, using the facts that $\tilde{\mathbf{u}}_h \in \mathbf{H}(\text{div}, \Omega)$, $\nabla \cdot \tilde{\mathbf{u}}_h = f$, and $\tilde{\mathbf{u}}_h \cdot \mathbf{n}_\sigma = \mathbf{U}_\sigma / |\sigma|$ for all $\sigma \in \mathcal{F}_K$ following (5.11a), and replacing Lemma 5.8 by Lemma 5.9. $\square$

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

Finally, relying on Lemma 5.10, we also have:

**Corollary 5.20** (A simple estimate with the given element matrices $\hat{\mathbf{A}}_K$). Let the assumptions of Corollary 5.19 hold, let $\mathcal{F}_{K,h}^{\text{ext}} = \mathcal{F}_K$, i.e., the simplicial faces of $\mathcal{T}_h$ do not subdivide the faces of the polytopal mesh $\mathcal{T}_H$, and let the flux $\mathbf{u}_h \in \mathcal{R}\mathcal{T}_0(\mathcal{T}_h) \cap \mathbf{H}(\text{div}, \Omega)$ be constructed by Definition 5.5. Then

$$ \| \mathbf{u} - \mathbf{u}_h \|_{\underline{\mathbf{K}}^{-\frac{1}{2}}} \lesssim \left\{ \sum_{K \in \mathcal{T}_H} \eta_K^2 \right\}^{1/2}, \quad (5.42) $$

where $\eta_K$ is given by (5.41) and the “approximately less than or equal to” denoted by $\lesssim$ stems from the approximation (5.24).

A few comments are in order:

**Remark 5.21** (Comparison of Theorem 5.18 with Corollaries 5.19 and 5.20). Theorem 5.18 is applicable to any scheme of the form of Assumption 5.1 and relies on the the mixed finite element matrices $\hat{\mathbf{A}}_{\text{MFE},K}$ that one needs to construct by (5.20). In contrast, Corollaries 5.19 and 5.20 are only applicable to schemes of the form of Assumption 5.2, but use directly the matrices $\hat{\mathbf{A}}_K$. The estimates (5.36) and (5.42) hold for the reconstruction $\mathbf{u}_h$ that one can compute in practice if needed by Definition 5.5, whereas (5.40) only holds for the generally unavailable lifting $\tilde{\mathbf{u}}_h$, as in [42, 43]. Finally, both (5.36) and (5.40) are guaranteed, with the $\le$ inequality, whereas (5.42) is not. In the numerical experiments in Section 5.10 below, though, it is hard to distinguish the two estimators (5.37) and (5.41).

**Remark 5.22** (Simplicial submeshes and solution of local problems). We can avoid the physical construction of the simplicial submeshes $\mathcal{T}_K$ under the conditions that 1) the given scheme takes the form of Assumption 5.2; 2) the element matrices $\hat{\mathbf{A}}_K$ are explicitly given; 3) the stiffness matrix $\hat{\mathbf{S}}_{\text{FE},K}$ of (5.32) and the mass matrix $\hat{\mathbf{M}}_{\text{FE},K}$ of (5.33) of the conforming finite element method on $\mathcal{T}_K$ for each polytopal
<figure><img src="image_22.png" /><figcaption>Figure 19: The exact solution (5.43)</figcaption></figure>

cell $K \in T_H$ are constructed from the geometry of $T_K$ only, in the sense of Remark 5.16; 4) we use the estimator $\eta_K$ (5.41) of Corollary 5.19 or 5.20. Note, however, that one may need a submesh of each polytopal cell $K$ already to produce the scheme element matrix $\hat{A}_K$. Similarly, a typical implementation of $\hat{A}_{MFE,K}$ (5.20) will also need the submesh $T_K$, though it can also be constructed from the geometry of $T_K$ only. Remark finally that behind (5.20), there is the local Neumann problem of Definition 5.5.

**Remark 5.23** (Higher-order locally conservative methods). Higher-order locally conservative methods can be treated similarly to the above exposition but are not subject of this work.

## 5.10 Numerical experiments

The purpose of this section is to numerically illustrate the performance of the estimators of Theorem 5.18, as well as the simplified estimate (5.42) of Corollary 5.20. We also compare this methodology to a reference estimate on the simplicial submesh in particular using the potential reconstruction of Remark 5.17. The test is taken from [188], where a collection of two-dimensional elliptic problems for testing adaptive grid refinement algorithms is proposed. We approximate $-\Delta p = f$ on the space domain $\Omega = (0, 1)^2$ where the analytical solution is

$$p(x, y) = 2^{4\alpha} x^{\alpha} (1-x)^{\alpha} y^{\alpha} (1-y)^{\alpha}, \quad \alpha = 200, \qquad (5.43)$$

see Figure 19. The source term $f$ is prescribed correspondingly and we neglect here that it is not piecewise constant, i.e., the data oscillation term as the second term in (4.26b). We use homogeneous Dirichlet boundary conditions.

The mesh $T_H$ consists of general polygonal elements; together with the fictitious triangular submesh $T_h$, it is shown in Figure 8. We consider the hybrid finite volume (HFV) discretization (Droniou et al. [128] or [110, Section 2.2]), taking the form (5.9) of Assumption 5.2 with the matrix $A$ formed by local element matrices $\hat{A}_K$. We compare three versions of a posteriori error estimates: a) Theorem 5.18, with the estimators $\eta_K$ evaluated via the matrices $\hat{A}_{MFE,K}$, $\hat{S}_{FE,K}$, and $\hat{M}_{FE,K}$ (polygonal MFE estimate); b) Corollary 5.20, using the element matrix $\hat{A}_K$ of the HFV scheme that is already available (polygonal HFV estimate); c) estimate (5.36) with the estimators $\eta_K$ replaced by the expression $\|K^{-1/2}u_h + K^{1/2}\nabla\zeta_h\|_K$, where $u_h|_K \in V_{h,N}^K$ is given in Definition 5.5 and $\zeta_h \in P_2(T_h) \cap H_0^1(\Omega)$ is described in Remark 5.17 (triangular estimate, sharpest-possible, for comparison); In the first two cases, we use Definition 5.11 to obtain the vectors of the potential point values $Z_K$ and $Z_K^{\text{ext}}$ by (5.25) and (5.26).

Note that if the mixed finite element method on polygonal meshes of [254, Theorem 7.2] was used instead of the HFV discretization, the third procedure would be, following [254, Remark 7.3], fully equivalent to solving the problem (5.1) directly on the simplicial mesh $T_h$ by the lowest-order Raviart-Thomas mixed finite element method and applying the a posteriori error estimates of Theorem 4.18. This is why the triangular estimate serves here as a reference a posteriori error estimate. The first procedure only uses the piecewise affine potential reconstruction $\zeta_h$ of Definitions 5.11 and 5.14, which allows for the simple matrix form (5.37) of the estimators. The second procedure is definitely the easiest choice in practice, where only the already available element matrices $\hat{A}_K$ are used and there is no need to construct the mixed finite element matrices $\hat{A}_{MFE,K}$ via (5.20).
<figure><img src="image_23.png" /><figcaption>Figure 20: Actual and estimated error distributions, entire domain (top) and center zoom (bottom)</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 21: Relative error $\frac{\|u-u_h\|}{\|u_h\|}$ and relative estimators and $\frac{\eta}{\|u_h\|}$ (left) and effectivity indices (right), uniform mesh refinement</figcaption></figure>

In Figure 20 we compare the actual and predicted error distributions. The energy error and the sharp triangular estimate distributions match perfectly. The polygonal MFE estimate and the polygonal HFV estimates give similar results and match also well with the energy error, though the error distribution around the peak is less well captured. We depict in Figure 21 the error and estimates as a function of the total number of unknowns and the corresponding effectivity indices for a uniform mesh refinement. This is in this test performed as follows: we refine uniformly the triangular submesh $T_h$, giving actually rise to a Delaunay triangulation on each step, and then we merge triangles into polygons. All the three estimators behave in a similar way, with a slight advantage for the triangular estimate. The graphs confirm in particular that replacing the mixed finite element matrix $\hat{A}_{MFE,K}$ in (5.37) by $\hat{A}_K$ in (5.41) has a very small influence.

Figure 22 shows the results for adaptive mesh refinement, achieved by using the local distribution of the predicted error as an indicator to refine only the cells of the mesh where the error is important. More precisely, we refine the cells $K \in \mathcal{T}_H$ such that $\eta_K \ge 0.7 \max_{L \in \mathcal{T}_H} \eta_L$. We observe quasi-identical values of the polygonal MFE and polygonal HFV estimates; consequently, we obtain the same number of unknowns at each step of adaptivity and almost identical final energy error for these two estimators. The left part of Figure 22 displays the results where we present the estimators and two energy errors, one resulting with adaptivity based on the polygonal estimates and another with adaptivity based on the triangular estimates. Almost undistinguishable results are obtained. Finally, adaptive mesh refinement
<figure><img src="image_25.png" /><figcaption>Figure 22: Relative error $\frac{\|u-u_h\|}{\|u_h\|}$ and relative estimators and $\frac{\eta}{\|u_h\|}$ (left) and effectivity indices (right), adaptive mesh refinement</figcaption></figure>

leads to both smaller error and better effectivity indices in comparison with uniform mesh refinement.

## 5.11 Bibliographic resources

Literature on a posteriori analysis on polytopal meshes is much less plentiful than that discussed in Section 4.12 for simplicial meshes. Beirão da Veiga [42] and Beirão da Veiga and Manzini [43] derive a posteriori error estimates for low-order mimetic finite difference methods; extensions are presented in Antonietti et al. [22]. Omnes et al. [206] considered the discrete duality finite volume method. Recent contributions on the subject include Munar et al. [193] or Li et al. [180], see also the references therein. Here we have followed Vohralík and Yousef [255], developing the idea of the local Neumann problems (5.14) from Vohralík [249, Section 3.2]. Numerous additional results were obtained for more complex model problems, which we discuss in Section 7.12 below.

# 6 Steady nonlinear pure diffusion problems. Iterative linearization and algebraic errors

As a passage between the model steady linear Darcy flow of Section 5 and the complex multiphase compositional Darcy flow of Section 7, we consider here a steady nonlinear Darcy flow. We showcase that our methodology applies here in a similar simple fashion as in Section 5. We in particular give all the details how inexact linearization and linear solvers can be taken into account, following the concept in Ern and Vohralík [124], giving rise to a simple *adaptive inexact Newton method* on polytopal meshes. Its particularity is that it gives a guaranteed upper bound on the total error in the fluxes on each linearization and each algebraic solver step, distinguishes the different error components, and is obtained by simple multiplications of the same element matrices as in the linear case by vectors of face normal fluxes and potential values immediately available in each lowest-order locally conservative polytopal scheme.

## 6.1 Singlephase steady nonlinear Darcy flow

For simplicity of exposition, let us consider the following quasi-linear version of (5.1): find $p: \Omega \to \mathbb{R}$ such that

$$-\nabla \cdot (\underline{\mathbf{K}}(|\nabla p|)\nabla p) = f \quad \text{in } \Omega, \qquad (6.1a)$$

$$p = 0 \quad \text{on } \partial\Omega. \qquad (6.1b)$$

We suppose that the nonlinearity can be inverted in the sense that

$$v = -\underline{\mathbf{K}}(|w|)w \iff w = -\tilde{\underline{\mathbf{K}}}(|v|)v \qquad (6.2)$$

for all $v, w \in \mathbb{R}^d$. We suppose strong monotonicity and Lipschitz-continuity, i.e., that there exist two positive constants $c_{\underline{\mathbf{K}}}, C_{\underline{\mathbf{K}}}$ so that for all $v, w \in \mathbb{R}^d$,

$$c_{\underline{\mathbf{K}}} |v - w|^2 \le (v - w) \cdot (\tilde{\underline{\mathbf{K}}}(|v|)v - \tilde{\underline{\mathbf{K}}}(|w|)w), \qquad (6.3a)$$

$$|\tilde{\underline{\mathbf{K}}}(|v|)v - \tilde{\underline{\mathbf{K}}}(|w|)w| \le C_{\underline{\mathbf{K}}} |v - w|. \qquad (6.3b)$$
<figure><img src="image_26.png" /><figcaption>Figure 23: Prototypical monotonous and Lipschitz-continuous function $\tilde{k}(r) := c_{\tilde{K}} + \frac{C_{\tilde{K}} - c_{\tilde{K}}}{\sqrt{1+r^2}}$ for $c_{\tilde{K}} = 1$ and different values of $C_{\tilde{K}}$ (left), $\tilde{k}(r)r$ (middle), and $(\tilde{k}(r)r)'$ (right)</figcaption></figure>

Moreover, we let $\tilde{K}$ and $\tilde{\tilde{K}}$ take symmetric values and, namely to arrive at a simple matrix-vector multiplication form of the estimates, we also suppose that, for all $v, w \in \mathbb{R}^d$,

$$c_{\tilde{K}}|v|^2 \le v \cdot \tilde{K}(|w|)v, \quad |\tilde{K}(|w|)v| \le C_{\tilde{K}}|v|. \qquad (6.4)$$

Note that for a linear problem, $\tilde{K} = K^{-1}$, $c_{\tilde{K}}$ is simply the smallest eigenvalue of $\tilde{K}$, i.e., the reciprocal of the largest eigenvalue of $K$, and $C_{\tilde{K}}$ is the largest eigenvalue of $\tilde{K}$, i.e., the reciprocal of the smallest eigenvalue of $K$.

A prototypical example satisfying (6.3)-(6.4) is, for two positive constant $c_{\tilde{K}}, C_{\tilde{K}}$,

$$\tilde{k}(r) := c_{\tilde{K}} + \frac{C_{\tilde{K}} - c_{\tilde{K}}}{\sqrt{1+r^2}}, \quad \tilde{K} := \tilde{k}I. \qquad (6.5)$$

In this case, it can be easily verified that

$$c_{\tilde{K}} \le \tilde{k}(r) \le C_{\tilde{K}}, \quad c_{\tilde{K}} \le (\tilde{k}(r)r)' \le C_{\tilde{K}}.$$

An illustration is provided in Figure 23. A much wider class of nonlinearities is considered below in Section 7.1.

## 6.2 Weak solution and its properties

In a generalization of Section 5.2, the weak formulation of (6.1) looks for $p \in H_0^1(\Omega)$ such that

$$(\underline{K}(|\nabla p|)\nabla p, \nabla v) = (f, v) \quad \forall v \in H_0^1(\Omega). \qquad (6.6)$$

From the pressure head $p$, the Darcy velocity is given by

$$u := -\underline{K}(|\nabla p|)\nabla p, \qquad (6.7)$$

and

$$\nabla p = -\tilde{K}(|u|)u \qquad (6.8)$$

from (6.2). Still as in Proposition 4.1, it follows from (6.6) that

$$u \in H(\mathrm{div}, \Omega), \quad \nabla \cdot u = f. \qquad (6.9)$$

## 6.3 Total error and its spatial discretization, linearization, and algebraic components

In the Poisson model problem in Section 4, we measure the error as $||u - u_h||$, see (4.5). In the Darcy problem in Section 5, the intrinsic error measure becomes $||u - u_h||_{K^{-1/2}}$, see (5.5). From the mathematical setting in Section 6.1, a natural extension is here to consider

$$c_{\tilde{K}}^{1/2} ||u - u_h^{k,i}|| \qquad (6.10)$$
as the error, where the weight $c_{\tilde{K}}$ in front of the $L^2(\Omega)$-norm is the monotonicity constant from (6.3a).

In (6.10), we have denoted the approximate solution by $u_h^{k,i}$, anticipating that our numerical approximation will here consist in spatial discretization (with characteristic mesh size $h$), iterative linearization (with iteration index $k$), and iterative algebraic resolution (with iteration index $i$). Our goal will be to distinguish the error components as in (2.14) in Section 2.4.6 and to design an a posteriori error estimate of the form (2.23a), following the general concepts outlined in Section 2.5.

## 6.4 Generic discretizations on a polytopal mesh, iterative linearization, and iterative algebraic system solution

Proceeding as generally as possible, we accommodate any discretization of the form of Assumption 5.1. Let $F_K := (f, 1)_K$ for all $K \in \mathcal{T}_H$, as in (5.7). For the model nonlinear problem (6.1), the flux balance (5.8) now takes the form: find the algebraic vector $P := \{P_K\}_{K \in \mathcal{T}_H} \in \mathbb{R}^{|\mathcal{T}_H|}$ such that

$$ \sum_{\sigma \in \mathcal{F}_K} (U(P))_{\sigma} n_K \cdot n_{\sigma} = F_K \quad \forall K \in \mathcal{T}_H. \qquad (6.11) $$

Here $(U(P))_{\sigma} \in \mathbb{R}$ for each face $\sigma \in \mathcal{F}_H$ again approximates the normal flux $\langle u \cdot n_{\sigma}, 1 \rangle_{\sigma}$ over the face $\sigma$, but it depends in a nonlinear way on the pressures $P$, which is denoted by $U(P)$. It follows that (6.11) is a system of $|\mathcal{T}_H|$ nonlinear algebraic equations for the $|\mathcal{T}_H|$ unknowns $P$.

The solution to (6.11) is typically sought by an iterative linearization such as the Newton or the fixed point method. Given an initial guess $P^0$, on each linearization step $k \ge 1$, this leads to a system of linear algebraic equations of the form: find the algebraic vector $P^k := \{P_K^k\}_{K \in \mathcal{T}_H} \in \mathbb{R}^{|\mathcal{T}_H|}$ such that

$$ \sum_{\sigma \in \mathcal{F}_K} (U^{k-1}(P^k))_{\sigma} n_K \cdot n_{\sigma} = F_K \quad \forall K \in \mathcal{T}_H. \qquad (6.12) $$

Here, $U^{k-1}(P^k)$ are linearized face normal fluxes. These are obtained from the previous iterate $P^{k-1}$ and $U(P^{k-1})$ but depend in a linear (affine) way on $P^k$. For the Newton linearization, in particular,

$$ (U^{k-1}(P^k))_{\sigma} := \sum_{K' \in \mathcal{T}_H} \frac{\partial(U(P^{k-1}))_{\sigma}}{\partial P_{K'}} \cdot (P_{K'}^k - P_{K'}^{k-1}) + (U(P^{k-1}))_{\sigma}. \qquad (6.13) $$

In Section 7 below, (7.27) provides a general Newton linearization for a more involved problem.

The solution to (6.12) is typically sought by an iterative algebraic solver. Given an initial guess $P^{k,0}$, typically the last algebraic solver iterate on the previous linearization step $k-1$ when $k > 1$, on each algebraic solver step $i \ge 1$, the constant

$$ (R)_K^{k,i} := F_K - \sum_{\sigma \in \mathcal{F}_K} (U^{k-1}(P^{k,i}))_{\sigma} n_K \cdot n_{\sigma} \qquad (6.14) $$

is the algebraic residual on the element $K$, the misfit of linearized fluxes $U^{k-1}(P^{k,i})$ not to satisfy (6.12). Consequently, the flux balance on linearization step $k$ and algebraic step $i$ takes the form

$$ \sum_{\sigma \in \mathcal{F}_K} (U^{k-1}(P^{k,i}))_{\sigma} n_K \cdot n_{\sigma} = F_K - (R)_K^{k,i} \quad \forall K \in \mathcal{T}_H, \qquad (6.15) $$

i.e., (6.12) up to the algebraic residual. In order to estimate the algebraic error, we will, following [140, 154, 124], also employ $j \ge 1$ additional algebraic solver steps in (6.15), giving rise to

$$ \sum_{\sigma \in \mathcal{F}_K} (U^{k-1}(P^{k,i+j}))_{\sigma} n_K \cdot n_{\sigma} = F_K - (R)_K^{k,i+j} \quad \forall K \in \mathcal{T}_H. \qquad (6.16) $$

The motivation is to take $j$ large enough to make the algebraic residual $(R)_K^{k,i+j}$ considerably smaller than $(R)_K^{k,i}$, see Remark 6.3 below.

**Example 6.1** (Discretization of problem (6.1) by the finite volume scheme on a simplicial mesh with fixed-point linearization and an arbitrary iterative algebraic solver). Consider the setting and notation
of Section 4.7. The cell-centered finite volume scheme for the nonlinear Darcy problem (6.1) on the simplicial mesh $T_h$ reads: find the real values $p_K$, $K \in T_h$, the approximations to the mean values of $p$ in the mesh elements $K$, such that

$$ \sum_{\sigma \in \mathcal{F}_K} \xi_{K,\sigma} U_{K,\sigma} = (f, 1)_K \quad \forall K \in T_h, \qquad (6.17a) $$

where the real number $\xi_{K,\sigma} U_{K,\sigma}$ for each face $\sigma \in \mathcal{F}_K$ approximates the normal (out)flux $\langle \mathbf{u} \cdot \mathbf{n}_K, 1 \rangle_{\sigma} = \langle (-\mathbf{K}(|\nabla p|)\nabla p) \cdot \mathbf{n}_K, 1 \rangle_{\sigma}$ from $K$ over the face $\sigma$. This is composed of the approximation of $\langle -\nabla p \cdot \mathbf{n}_K, 1 \rangle_{\sigma}$ given, as in (4.19b)-(4.19c), by

$$ U_{K,\sigma} := -\frac{|\sigma_{K,L}|}{d_{K,L}}(p_L - p_K) \quad \sigma = \sigma_{K,L} \in \mathcal{F}_h^{\text{int}}, \qquad (6.17b) $$

$$ U_{K,\sigma} := -\frac{|\sigma|}{d_{K,\sigma}}(0 - p_K) \quad \sigma \in \mathcal{F}_K \cap \mathcal{F}_h^{\text{ext}} \qquad (6.17c) $$

and of the approximation of $\mathbf{K}(|\nabla p|)$ on $\sigma$ given by

$$ \xi_{K,\sigma} := \underline{\mathbf{K}}\left(\left|\frac{1}{2}(\xi_h|_K(\mathbf{x}_{K,\sigma}) + \xi_h|_L(\mathbf{x}_{K,\sigma}))\right|\right) \quad \sigma = \sigma_{K,L} \in \mathcal{F}_h^{\text{int}}, \qquad (6.17d) $$

$$ \xi_{K,\sigma} := \underline{\mathbf{K}}(|\xi_h|_K(\mathbf{x}_{K,\sigma})|) \quad \sigma \in \mathcal{F}_K \cap \mathcal{F}_h^{\text{ext}}, \qquad (6.17e) $$

where the approximation of $-\nabla p$ is, as in (4.22), given by

$$ (\xi_h|_K)(\mathbf{x}) := \sum_{\sigma \in \mathcal{F}_K} U_{K,\sigma} \frac{1}{d|K|}(\mathbf{x} - \mathbf{a}_{K,\sigma}), \quad \mathbf{x} \in K, K \in T_h, \qquad (6.17f) $$

see the illustrations in Figure 17. Note that in contrast to (4.19b)-(4.19c), the face fluxes $\xi_{K,\sigma} U_{K,\sigma}$ are now nonlinear functions of the pressure unknowns $p_K$ and $p_L$ (unless $\mathbf{K}(\cdot)$ is a constant function). Also note that (6.17) writes in the form (6.11) with $(\mathbf{U}(\mathbf{P}))_{\sigma} = \xi_{K,\sigma} U_{K,\sigma} \mathbf{n}_K \cdot \mathbf{n}_{\sigma}$ (recall that $\mathbf{n}_K \cdot \mathbf{n}_{\sigma} = \pm 1$ only determines the sign).

Iterative linearization of (6.17) by the fixed point method then reads: given the initial values $p_K^0$, $K \in T_h$, on each linearization step $k \ge 1$, find the real values $p_K^k$, $K \in T_h$, such that

$$ \sum_{\sigma \in \mathcal{F}_K} \xi_{K,\sigma}^{k-1} U_{K,\sigma}^k = (f, 1)_K \quad \forall K \in T_h, \qquad (6.18a) $$

where

$$ U_{K,\sigma}^k := -\frac{|\sigma_{K,L}|}{d_{K,L}}(p_L^k - p_K^k) \quad \sigma = \sigma_{K,L} \in \mathcal{F}_h^{\text{int}}, \qquad (6.18b) $$

$$ U_{K,\sigma}^k := -\frac{|\sigma|}{d_{K,\sigma}}(0 - p_K^k) \quad \sigma \in \mathcal{F}_K \cap \mathcal{F}_h^{\text{ext}}, \qquad (6.18c) $$

$$ \xi_{K,\sigma}^{k-1} := \underline{\mathbf{K}}\left(\left|\frac{1}{2}(\xi_h^{k-1}|_K(\mathbf{x}_{K,\sigma}) + \xi_h^{k-1}|_L(\mathbf{x}_{K,\sigma}))\right|\right) \quad \sigma = \sigma_{K,L} \in \mathcal{F}_h^{\text{int}}, \qquad (6.18d) $$

$$ \xi_{K,\sigma}^{k-1} := \underline{\mathbf{K}}(|\xi_h^{k-1}|_K(\mathbf{x}_{K,\sigma})|) \quad \sigma \in \mathcal{F}_K \cap \mathcal{F}_h^{\text{ext}}, \qquad (6.18e) $$

$$ (\xi_h^{k-1}|_K)(\mathbf{x}) := \sum_{\sigma \in \mathcal{F}_K} U_{K,\sigma}^{k-1} \frac{1}{d|K|}(\mathbf{x} - \mathbf{a}_{K,\sigma}), \quad \mathbf{x} \in K, K \in T_h. \qquad (6.18f) $$

We see that (6.18) takes the form requested in (6.12).

Finally, (6.15) for an arbitrary iterative algebraic solver follows by the definition (6.14).

## 6.5 Face normal fluxes and potential point values

The key for the guaranteed a posteriori error estimate with inexpensive implementation and evaluation for the linear Darcy problem (5.1) in Theorem 5.18 were the face normal fluxes of Section 5.5 and the potential point values of Section 5.7. We will now identify the face normal fluxes and potential point values for the present nonlinear Darcy problem (6.1). In order to devise a guaranteed a posteriori error
estimate distinguishing the error components, we will also identify linearization error and algebraic error face fluxes.

Let a polytopal cell $K \in \mathcal{T}_H$ with its face $\sigma \in \mathcal{F}_K$ be fixed. For a given linearization step $k \ge 1$, algebraic step $i \ge 1$, and $j \ge 1$ additional algebraic iterations, we set

$$ (U_K^{k,i})_{\sigma} := (U(P^{k,i}))_{\sigma}, \qquad (6.19a) $$

$$ (U_{\text{lin},K}^{k,i})_{\sigma} := (U^{k-1}(P^{k,i}))_{\sigma} - (U(P^{k,i}))_{\sigma}, \qquad (6.19b) $$

$$ (U_{\text{alg},K}^{k,i})_{\sigma} := (U^{k-1}(P^{k,i+j}))_{\sigma} - (U^{k-1}(P^{k,i}))_{\sigma}. \qquad (6.19c) $$

Observe that $U^{k-1}(P^{k,i+j})$ are directly available as employed in (6.16) and that $U^{k-1}(P^{k,i})$ are directly available as employed in (6.15). In turn, $U(P^{k,i})$ are available upon plugging the current approximation $P^{k,i}$ into the nonlinear flux function of (6.11). We call $U_K^{k,i}$ discretization face fluxes since in the limit $k, i \to \infty$, they tend to $U(P)$ of (6.11). Similarly, we call $U_{\text{lin},K}^{k,i}$ linearization error face fluxes since they vanish with linearization and algebraic solver iterations $k, i \to \infty$. Finally, we call $U_{\text{alg},K}^{k,i}$ the algebraic error face fluxes since they vanish with algebraic solver iterations $i \to \infty$. Observe that

$$ (U_K^{k,i})_{\sigma} + (U_{\text{lin},K}^{k,i})_{\sigma} + (U_{\text{alg},K}^{k,i})_{\sigma} = (U^{k-1}(P^{k,i+j}))_{\sigma}, $$

so that from (6.16)

$$ \sum_{\sigma \in \mathcal{F}_K} \left( (U_K^{k,i})_{\sigma} + (U_{\text{lin},K}^{k,i})_{\sigma} + (U_{\text{alg},K}^{k,i})_{\sigma} \right) n_K \cdot n_{\sigma} = F_K - (R)_K^{k,i+j} \quad \forall K \in \mathcal{T}_H. \qquad (6.20) $$

The fluxes per simplex face are then obtained as in Definition 5.3.

The potential point values are obtained as in Definition 5.11 or 5.13 directly from the cell pressure heads $P^{k,i}$. This gives the vectors $Z_K^{k,i}$ and $Z_K^{k,i,\text{ext}}$ of potential point values in the vertices of the simplicial submesh $\mathcal{T}_K$ of each polytopal element $K \in \mathcal{T}_H$.

## 6.6 Fictitious flux and potential reconstructions

We now extend the reconstructions of Sections 5.6 and 5.8 to the present nonlinear setting.

We have to face the fact that the diffusion tensor $\underline{K}(|\nabla p|)$ is now unknown since $p$ is unknown. We thus modify Definition 5.5 to

$$ \underline{u}_h^{k,i}|_K := c_{\underline{K}}^{-1} C_{\underline{K}}^2 \arg \min_{v_h \in V_{h,N}, \nabla \cdot v_h = \text{constant}} \|v_h\|_K^2, \qquad (6.21) $$

where the data in the spaces $V_{h,N}^K$ are now given by $U_K^{k,i}$ from (6.19a) in place of $U_\sigma$ (we systematically use the convention (5.12) (without change of the notation) to define the fluxes of the simplicial submesh faces in case they subdivide the original polytopal faces). Note that we cannot in general impose here $\nabla \cdot (\underline{u}_h^{k,i})|_K = f|_K$. This discretization flux reconstruction is again not needed to evaluate the error $c_{\underline{K}}^{1/2} \|\underline{u} - \underline{u}_h^{k,i}\|_K$.) For appropriate scaling, we have replaced $\underline{K}^{-1}$ by $c_{\underline{K}}^{-1} C_{\underline{K}}^2 \underline{I}$. It is important to stress that just like in (5.16), problem (6.21) is linear and can be written in a matrix form as in (5.18), with the same matrices, just replacing $\underline{K}^{-1}$ by $c_{\underline{K}}^{-1} C_{\underline{K}}^2 \underline{I}$. This in particular leads to the definition of the MFE element matrix $\hat{A}_{\text{MFE},K}$ by (5.20) also in the present nonlinear setting, up to the replacement of $\underline{K}^{-1}$ by $c_{\underline{K}}^{-1} C_{\underline{K}}^2 \underline{I}$.

Similarly, we define the linearization error flux reconstruction $\underline{u}_{\text{lin},h}^{k,i}$ by the lifting (6.21), employing $U_{\text{lin},K}^{k,i}$ from (6.19b) in the definition of the space $V_{h,N}^K$, and the algebraic error flux reconstruction $\underline{u}_{\text{alg},h}^{k,i}$ by the lifting (6.21), employing $U_{\text{alg},K}^{k,i}$ from (6.19c) in the definition of the space $V_{h,N}^K$. The flux reconstructions $\underline{u}_{\text{lin},h}^{k,i}$ and $\underline{u}_{\text{alg},h}^{k,i}$ will only be used as a theoretical vehicle below. From (6.20), we observe by the Green theorem that

$$ (\nabla \cdot (\underline{u}_h^{k,i} + \underline{u}_{\text{lin},h}^{k,i} + \underline{u}_{\text{alg},h}^{k,i}))|_K = f|_K - |K|^{-1}(R)_K^{k,i+j} \quad \forall K \in \mathcal{T}_H. \qquad (6.22) $$

Finally, the potential reconstruction $\zeta_h^{k,i}$ is simply obtained from the element vectors $Z_K^{k,i}$ and $Z_K^{k,i,\text{ext}}$ assembled in Section 6.5 as in Definition 5.14, i.e., setting

$$ \zeta_h^{k,i}(a) := Z_a^{k,i} \quad \forall a \in V_h. $$
Congruently, for any $K \in \mathcal{T}_H$, we define the stiffness matrix $\hat{\mathbb{S}}_{\text{FE},K}$ as in (5.32), with $c_{\tilde{K}}^{-2} C_{\tilde{K}} \mathbf{I}$ in place of $\mathbf{K}$, whereas the mass matrix $\hat{\mathbb{M}}_{\text{FE},K}$ is defined exactly as in (5.33).

## 6.7 A guaranteed a posteriori error estimate distinguishing the error components with inexpensive implementation and evaluation

Let $C_F$ be the constant from the Friedrichs inequality

$$ \|v\| \le C_F h_{\Omega} \|\nabla v\| \quad \forall v \in H_0^1(\Omega); \qquad (6.23) $$

there holds $1/(\pi d) \le C_F \le 1$, as $C_F h_{\Omega}$ is the square root of the reciprocal of the smallest eigenvalue of the Laplace operator on $\Omega$ with homogenous Dirichlet boundary condition. The main result of this section is:

**Theorem 6.2** (A guaranteed a posteriori error estimate distinguishing the error components with inexpensive implementation and evaluation).

Let $f$ be constant on each $K \in \mathcal{T}_H$ and let $u$ be given by (6.6)-(6.7). Consider any polytopal discretization of the form (6.11), any iterative linearization (6.12) on step $k \ge 1$, and any iterative algebraic solver (6.15) on step $i \ge 1$. Consider $j \ge 1$ additional algebraic solver steps leading to (6.16). Let, for each polytopal element $K \in \mathcal{T}_H$, the element vectors of face normal fluxes $U_K^{k,i}$, $U_{\text{lin},K}^{k,i}$, and $U_{\text{alg},K}^{k,i}$ be given by (6.19). Let the element vectors of potential point values $Z_K^{k,i}$ and $Z_K^{k,i,\text{ext}}$ be constructed from the values $P^{k,i}$ following Section 6.5. Let finally the element matrices $\hat{\mathbb{A}}_{\text{MFE},K}$ and $\hat{\mathbb{S}}_{\text{FE},K}$ be defined by (5.20) and (5.32), where one respectively takes the multiples of the identity matrix $c_{\tilde{K}}^{-1} C_{\tilde{K}}^2 \mathbf{I}$ in place of $\mathbf{K}^{-1}$ and $c_{\tilde{K}}^{-2} C_{\tilde{K}} \mathbf{I}$ in place of $\mathbf{K}$. Let also $\hat{\mathbb{M}}_{\text{FE},K}$ be given by (5.33). Then there holds

$$ c_{\tilde{K}}^{1/2} \|u - u_h^{k,i}\| \le \eta_{\text{sp}}^{k,i} + \eta_{\text{lin}}^{k,i} + \eta_{\text{alg}}^{k,i} + \eta_{\text{rem}}^{k,i} \qquad (6.24) $$

with

$$ \eta_{\bullet}^{k,i} = \left\{ \sum_{K \in \mathcal{T}_H} \left( \eta_{\bullet,K}^{k,i} \right)^2 \right\}^{1/2}, \quad \bullet = \{\text{sp}, \text{lin}, \text{alg}, \text{rem}\}, $$

where the (spatial) discretization estimators are given by

$$ (\eta_{\text{sp},K}^{k,i})^2 := (U_K^{k,i})^t \hat{\mathbb{A}}_{\text{MFE},K} U_K^{k,i} + (Z_K^{k,i})^t \hat{\mathbb{S}}_{\text{FE},K} Z_K^{k,i} + 2c_{\tilde{K}}^{-1} C_{\tilde{K}} \left[ (U_K^{k,i,\text{ext}})^t Z_K^{k,i,\text{ext}} - F_K |K|^{-1} 1^t \hat{\mathbb{M}}_{\text{FE},K} Z_K^{k,i} \right], $$

the linearization estimators by

$$ (\eta_{\text{lin},K}^{k,i})^2 := (U_{\text{lin},K}^{k,i})^t \hat{\mathbb{A}}_{\text{MFE},K} U_{\text{lin},K}^{k,i}, $$

the algebraic estimators by,

$$ (\eta_{\text{alg},K}^{k,i})^2 := (U_{\text{alg},K}^{k,i})^t \hat{\mathbb{A}}_{\text{MFE},K} U_{\text{alg},K}^{k,i}, $$

and the algebraic remainder estimators by

$$ \eta_{\text{rem},K}^{k,i} := c_{\tilde{K}}^{-1/2} C_{\tilde{K}} C_F h_{\Omega} |K|^{-1/2} |(R_K^{k,i+j})|. $$

Here, the flux reconstruction $u_h^{k,i} \in \mathcal{R}\mathcal{T}_0(\mathcal{T}_h) \cap \mathcal{H}(\text{div}, \Omega)$ is obtained following (6.21).

**Remark 6.3** (Algebraic remainder estimators $\eta_{\text{rem},K}^{k,i}$, choice of $j$ additional linear solver iterations, and mass balance). From (6.20) or (6.22), at “convergence” of the iterative algebraic solver (or for a direct solver)(recall that we neglect rounding errors), as the algebraic residuals $(R_K^{k,i+j})$ vanish, we recover exact mass balance on each linearization step $k \ge 1$. More importantly, for a sufficiently high number $j$ of additional algebraic iterations, we “almost” have the mass balance on each linearization step $k \ge 1$. We typically choose $j$ adaptively, following [140, 154, 124], so that this mass balance misfit related to the algebraic residuals $(R_K^{k,i+j})$ is negligible: we choose $j$ such that $\eta_{\text{rem}}^{k,i}$ becomes negligible in comparison
with the other estimators in (6.24). The presence of the terms $(\mathbb{R})_K^{k,i+j}$ is the disadvantage of the j additional iterations. In turn, the advantage of the j of additional algebraic iterations is the simplicity of the formula (6.19c). A more involved choice can be made following [210]. Then no j additional iterations are needed, the term $|K|^{-1}(\mathbb{R})_K^{k,i+j}$ in (6.20) and (6.22) and the term $\eta_{\text{rem}}^{k,i}$ in (6.24) vanish, there is an exact mass balance on each linearization step k ≥ 1 and each linear algebraic step i ≥ 1, and the a posteriori bound (6.24) is typically slightly more precise. Remark finally that these two choices respectively correspond to nonzero and zero $\rho_h^{k,i}$ in the general discussion in Section 2.7.

*Proof.* The proof follows in spirit of previous works by Kim [163], Jiránek et al. [154], and Ern and Vohralík [124, 125]. Let $\zeta \in H_0^1(\Omega)$ be the solution of

$$ (\underline{\mathbf{K}}(|\nabla\zeta|)\nabla\zeta, \nabla v) = -(\mathbf{u}_h^{k,i}, \nabla v) \quad \forall v \in H_0^1(\Omega) \qquad (6.25) $$

and set $\zeta := -\underline{\mathbf{K}}(|\nabla\zeta|)\nabla\zeta$, similarly to (6.7). Remark that if there holds $\mathbf{u}_h^{k,i} \in \mathbf{H}(\text{div}, \Omega)$ with $\nabla \cdot \mathbf{u}_h^{k,i} = f$ (if (6.11) is satisfied exactly, no iterative linearization and exact algebraic solve), then, by the Green theorem, $-(\mathbf{u}_h^{k,i}, \nabla v) = (f, v)$, so that $\zeta$ coincides with the solution p of (6.6).

By the triangle inequality,

$$ c_{\underline{\mathbf{K}}}^{1/2} \|\mathbf{u} - \mathbf{u}_h^{k,i}\| \le c_{\underline{\mathbf{K}}}^{1/2} \|\mathbf{u} - \zeta\| + c_{\underline{\mathbf{K}}}^{1/2} \|\zeta - \mathbf{u}_h^{k,i}\|. \qquad (6.26) $$

For the first term above, the strong monotonicity implies

$$
\begin{align*}
c_{\underline{\mathbf{K}}} \|\mathbf{u} - \zeta\|^2 &\le (\underline{\mathbf{K}}(|\nabla p|)\nabla p - \underline{\mathbf{K}}(|\nabla \zeta|)\nabla \zeta, \nabla(p - \zeta)) && \text{by (6.3a) and (6.2)} \\
&= (\underline{\mathbf{K}}(|\nabla p|)\nabla p + \mathbf{u}_h^{k,i}, \nabla(p - \zeta)) && \text{by (6.25)} \\
&= (f, p - \zeta) + (\mathbf{u}_h^{k,i}, \nabla(p - \zeta)). && \text{by (6.6)}
\end{align*}
$$

Note that by the Lipschitz-continuity (6.3b) together with (6.2),

$$ \|\nabla(p - \zeta)\| \le C_{\underline{\mathbf{K}}} \|\mathbf{u} - \zeta\|. \qquad (6.27) $$

Thus,

$$
\begin{align}
c_{\underline{\mathbf{K}}}^{1/2} \|\mathbf{u} - \zeta\| &= \frac{(f, p - \zeta) + (\mathbf{u}_h^{k,i}, \nabla(p - \zeta))}{c_{\underline{\mathbf{K}}}^{1/2} \|\mathbf{u} - \zeta\|} \tag{6.28} \\
&\le c_{\underline{\mathbf{K}}}^{-1/2} C_{\underline{\mathbf{K}}} \max_{\varphi \in H_0^1(\Omega), \|\nabla \varphi\|=1} \{(f, \varphi) + (\mathbf{u}_h^{k,i}, \nabla \varphi)\}. \nonumber
\end{align}
$$

To estimate (6.28) by a computable quantity, we use (6.22). Fix now $\varphi \in H_0^1(\Omega)$ with $\|\nabla \varphi\| = 1$. Then adding and subtracting $(\mathbf{u}_{\text{lin},h}^{k,i} + \mathbf{u}_{\text{alg},h}^{k,i}, \nabla \varphi)$ and using the Green theorem,

$$
\begin{align*}
(f, \varphi) + (\mathbf{u}_h^{k,i}, \nabla \varphi) &= (f - \nabla \cdot (\mathbf{u}_h^{k,i} + \mathbf{u}_{\text{lin},h}^{k,i} + \mathbf{u}_{\text{alg},h}^{k,i}), \varphi) - (\mathbf{u}_{\text{lin},h}^{k,i} + \mathbf{u}_{\text{alg},h}^{k,i}, \nabla \varphi) \\
&= \sum_{K \in \mathcal{T}_H} \{|K|^{-1}((\mathbb{R})_K^{k,i+j}, \varphi)_K - (\mathbf{u}_{\text{lin},h}^{k,i} + \mathbf{u}_{\text{alg},h}^{k,i}, \nabla \varphi)_K\}.
\end{align*}
$$

The Cauchy–Schwarz inequality gives

$$
\begin{align*}
\sum_{K \in \mathcal{T}_H} - (\mathbf{u}_{\text{lin},h}^{k,i} + \mathbf{u}_{\text{alg},h}^{k,i}, \nabla \varphi)_K &\le \sum_{K \in \mathcal{T}_H} \|\mathbf{u}_{\text{lin},h}^{k,i} + \mathbf{u}_{\text{alg},h}^{k,i}\|_K \|\nabla \varphi\|_K \le \|\mathbf{u}_{\text{lin},h}^{k,i} + \mathbf{u}_{\text{alg},h}^{k,i}\|_K \underbrace{\|\nabla \varphi\|}_{=1} \\
&\le \|\mathbf{u}_{\text{lin},h}^{k,i}\| + \|\mathbf{u}_{\text{alg},h}^{k,i}\|.
\end{align*}
$$

Similarly, again the Cauchy–Schwarz inequality leads to

$$ \sum_{K \in \mathcal{T}_H} |K|^{-1} |((\mathbb{R})_K^{k,i+j}, \varphi)_K| \le \sum_{K \in \mathcal{T}_H} |K|^{-1/2} |(\mathbb{R})_K^{k,i+j}| \|\varphi\|_K \le \left\{ \sum_{K \in \mathcal{T}_H} |K|^{-1} |(\mathbb{R})_K^{k,i+j}|^2 \right\}^{1/2} \|\varphi\|_{\Omega}. \quad (6.29) $$

Summing up these developments with the Friedrichs inequality (6.23), infer from (6.28)

$$ c_{\underline{\mathbf{K}}}^{1/2} \|\mathbf{u} - \zeta\| \le c_{\underline{\mathbf{K}}}^{-1/2} C_{\underline{\mathbf{K}}} \left[ C_F h \Omega \left\{ \sum_{K \in \mathcal{T}_H} |K|^{-1} |(\mathbb{R})_K^{k,i+j}|^2 \right\}^{1/2} + \|\mathbf{u}_{\text{lin},h}^{k,i}\| + \|\mathbf{u}_{\text{alg},h}^{k,i}\| \right]. $$
Finally, to evaluate the norm $||u_{\text{lin},h}^{k,i}||$, we use

$$||u_{\text{lin},h}^{k,i}||^2 = \sum_{K \in \mathcal{T}_H} ||u_{\text{lin},h}^{k,i}||_K^2$$

and employ the constant matrix $c_{\tilde{K}}^{-1}C_{\tilde{K}}^2I$ in place of $\tilde{K}^{-1}$ in Lemma 5.8 to infer

$$c_{\tilde{K}}^{-1}C_{\tilde{K}}^2 ||u_{\text{lin},h}^{k,i}||_K^2 = (U_{\text{lin},K}^{k,i})^t \hat{\mathbb{A}}_{\text{MFE},K} U_{\text{lin},K}^{k,i} \quad (6.30)$$

for each polytopal element $K \in \mathcal{T}_H$. We proceed similarly for $||u_{\text{alg},h}^{k,i}||$, which gives

$$c_{\tilde{K}}^{-1}C_{\tilde{K}}^2 ||u_{\text{alg},h}^{k,i}||_K^2 = (U_{\text{alg},K}^{k,i})^t \hat{\mathbb{A}}_{\text{MFE},K} U_{\text{alg},K}^{k,i}. \quad (6.31)$$

We finally estimate the second term in (6.26). For $v \in H_0^1(\Omega)$ arbitrary, the strong monotonicity in particular leads to

$$ \begin{aligned} c_{\tilde{K}} ||\zeta - u_h^{k,i}||^2 &\le (\tilde{K}(|\nabla \zeta|)\nabla \zeta + u_h^{k,i}, \nabla \zeta + \tilde{K}(|u_h^{k,i}|)u_h^{k,i}) & \text{by (6.3a) and (6.2)} \\ &= (\tilde{K}(|\nabla \zeta|)\nabla \zeta + u_h^{k,i}, \nabla v + \tilde{K}(|u_h^{k,i}|)u_h^{k,i}) & \text{by (6.25)} \\ &\le ||\zeta - u_h^{k,i}|| ||\tilde{K}(|u_h^{k,i}|)u_h^{k,i} + \nabla v|| & \text{by Cauchy-Schwarz} \end{aligned} $$

so that,

$$c_{\tilde{K}}^{1/2} ||\zeta - u_h^{k,i}|| \le c_{\tilde{K}}^{-1/2} \min_{v \in H_0^1(\Omega)} ||\tilde{K}(|u_h^{k,i}|)u_h^{k,i} + \nabla v|| \quad (6.32a)$$

$$\le c_{\tilde{K}}^{-1/2} ||\tilde{K}(|u_h^{k,i}|)u_h^{k,i} + \nabla \zeta_h^{k,i}|| \quad (6.32b)$$

for an arbitrary $\zeta_h^{k,i} \in H_0^1(\Omega)$. We pick the potential reconstruction from Section 6.6. For the inexpensive implementation and evaluation, we proceed as in Section 5.9. Consequently, for each polytopal mesh element $K \in \mathcal{T}_H$, employing here also the eigenvalue hypothesis (6.4),

$$ \begin{aligned} & ||\tilde{K}(|u_h^{k,i}|)u_h^{k,i} + \nabla \zeta_h^{k,i}||_K^2 \\ &= ||[\tilde{K}(|u_h^{k,i}|)]^{1/2} ([\tilde{K}(|u_h^{k,i}|)]^{1/2}u_h^{k,i} + [\tilde{K}(|u_h^{k,i}|)]^{-1/2}\nabla \zeta_h^{k,i})||_K^2 \\ &\le C_{\tilde{K}} ||[\tilde{K}(|u_h^{k,i}|)]^{1/2}u_h^{k,i} + [\tilde{K}(|u_h^{k,i}|)]^{-1/2}\nabla \zeta_h^{k,i}||_K^2 & \text{by (6.4)} \\ &= C_{\tilde{K}} \left( (\tilde{K}(|u_h^{k,i}|)u_h^{k,i}, u_h^{k,i})_K + 2(u_h^{k,i}, \nabla \zeta_h^{k,i})_K + ([\tilde{K}(|u_h^{k,i}|)]^{-1}\nabla \zeta_h^{k,i}, \nabla \zeta_h^{k,i})_K \right) \\ &\le C_{\tilde{K}}^2 ||u_h^{k,i}||_K^2 + 2C_{\tilde{K}}(u_h^{k,i}, \nabla \zeta_h^{k,i})_K + c_{\tilde{K}}^{-1}C_{\tilde{K}} ||\nabla \zeta_h^{k,i}||_K^2 & \text{by (6.4)} \end{aligned} $$

For the first term above, we proceed as in (6.30). For the second one, the development (5.39) remains unchanged. Finally, for the third term, we use (5.34) where $\tilde{K}$ has been replaced by $c_{\tilde{K}}^{-2}C_{\tilde{K}}I$ to see that

$$c_{\tilde{K}}^{-2}C_{\tilde{K}} ||\nabla \zeta_h^{k,i}||_K^2 = (Z_K^{k,i})^t \hat{\mathbb{S}}_{\text{FE},K} Z_K^{k,i}.$$

The proof is finished by combining the two estimates on $c_{\tilde{K}}^{1/2}||u - \zeta||$ and $c_{\tilde{K}}^{1/2}||\zeta - u_h^{k,i}||$. $\square$

**Remark 6.4** (Inexpensive implementation and evaluation). Note that the estimate of Theorem 6.2 still takes the same simple matrix-vector multiplication form of Theorem 5.18, with in particular the same element matrices $\hat{\mathbb{A}}_{\text{MFE},K}$, $\hat{\mathbb{S}}_{\text{FE},K}$, $\hat{\mathbb{M}}_{\text{FE},K}$ as in the linear case. Similarly to Corollaries 5.19 and 5.20, a further simplification is to replace $\hat{\mathbb{A}}_{\text{MFE},K}$ by the element matrices $\hat{\mathbb{A}}_K$ of the given scheme whenever they are available, as per Assumption 5.2.

**Remark 6.5** (Distinction of error components). The distinction of the different error components is based on the face fluxes (6.19), directly available from any finite volume discretization (6.11)-(6.16). Beyond simplicity, a remarkable property is that all the error component estimators have the same physical units (that of energy flux error) and take the same form. This stands in contrast to the usual practice where the algebraic error is typically treated by the (relative) $L^2$-norm of the algebraic residual (here, from (6.33), we rather employ the energy-scaled $H^{-1}(\Omega)$-norm of the residual) and the linearization error estimator is often the $L^\infty$-norm of the difference of two consecutive iterates. This unified treatment of the different error components is a consequence of our holistic approach. In practice, the distinction of error components is to be used in an adaptive inexact Newton algorithm such as Algorithm 7.1 of the next section.