Regular Article
Comparison of finitevolume schemes for diffusion problems
Institute for Modelling Hydraulic and Environmental Systems, University of Stuttgart, Pfaffenwaldring 61, 70569 Stuttgart, Germany
^{*} Correspondiong author: martin.schneider@iws.unistuttgart.de
Received:
13
February
2018
Accepted:
18
September
2018
We present an abstract discretization framework and demonstrate that various cellcentered and hybrid finitevolume schemes fit into it. The different schemes considered in this work are then analyzed numerically for an elliptic model problem with respect to the properties consistency, coercivity, extremum principles, and sparsity. The test cases presented comprise of two and threedimensional setups, mildly and highly anisotropic tensors and grids of different complexities. The results show that all schemes show a similar convergence behavior, except for the twopoint flux approximation scheme, and seem to be coercive. Furthermore, they confirm that linear schemes, in contrast to nonlinear schemes, are in general neither positivitypreserving nor satisfy discrete minimum or maximum principles.
© M. Schneider et al., published by IFP Energies nouvelles, 2018
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Numerical simulations play a crucial role in a wide range of geotechnical engineering applications, thus, it is of great importance that the schemes used to produce these results are reliable and robust. Additionally, efficient schemes are desirable due to the often very large spatial scales involved and the uncertain nature of the subsurface, which might require statistical analyses with a large number of model runs to be performed. Appropriate schemes have to be chosen applicationdependent subject to the particular requirements. Besides consistency, efficiency, reliability, and robustness, local mass conservation is essential when performing subsurface flow simulations. This is why the most commonly used schemes for subsurface flow simulations are cellcentered finitevolume methods, such as cellcentered Galerkin methods [1] or MultiPoint Flux Approximation (MPFA) methods [2–7]; or hybrid, mixed, and mimetic (HMM) schemes, such as the Hybrid FiniteVolume (HFV) schemes [8, 9], the Mixed FiniteElement (MFE) [10, 11] or the Mimetic FiniteDifference (MFD) methods [12, 13]. All these schemes are closely related and they can be either presented within a finitevolume framework (which is done in this article) or within a finiteelement framework. Such aspects have been presented in [8, 14–16]. The advantage of HMM methods is the fact that they can be applied to highly complex grids, e.g. cornerpoint grids. However, this comes with the cost of additional unknowns. Recently, monotone or discrete extremumprinciplespreserving schemes have been developed in [17–22], and applied to highly complex porous media applications in [23–25].
In this work, we compare different finitevolume schemes for an elliptic model problem regarding convergence, consistency, and efficiency. The latter is measured via the sparsity of the resulting linear systems of equations and the number of nonlinear solver iterations in the case of nonlinear schemes, while consistency is estimated based on discrete extremum principles and the linearity preservation property. Definitions of these properties as well as the numerical schemes considered in this work, covering a number of cellcentered and hybrid schemes, are presented in Section 2. The convergence of the schemes is investigated in Section 3.1, while in Section 3.2 the linearity preservation property is studied. The satisfaction of discrete extremum principles is considered in Section 3.3, before in Section 3.4 all schemes are applied to a threedimensional benchmark case. Finally, in Section 4 we summarize the results obtained in this article.
2 Discretization schemes
Let $\mathrm{\Omega}\subset {\mathbb{R}}^{d},d\in {\mathbb{N}}^{\mathrm{*}}$, be an open bounded connected polygonal domain with boundary ∂Ω, and ddimensional measure Ω. In the following, we consider the elliptic problem$$\{\begin{array}{cc}\nabla \xb7\left(\mathbf{\Lambda}\nabla u\right)=f& \mathrm{in}\mathrm{\Omega},\\ u=0& \mathrm{on}\partial \mathrm{\Omega}.\end{array}$$(1)
Here, we assume that f ∈ L ^{2}(Ω) and Λ is a symmetric tensorvalued function such that (s.t.) the spectrum of Λ(x) is contained in [α _{0}, β _{0}], with 0 < α _{0} < β _{0} < +∞, for almost every (a.e.) $\mathbf{x}\in \stackrel{\u0305}{\mathrm{\Omega}}$. These assumptions guarantee the existence of a weak solution ū (i.e. a solution of the variational formulation of problem (1)).
Remark 1. Within this section, homogeneous Dirichlet boundary conditions are assumed for ease of presentation. Other types of boundary conditions have for example been discussed within the gradient discretization framework [26].
In the following an admissible discretization is defined together with the notations that are used within this article.
Definition 1 (Admissible discretization). An admissible discretization $\mathcal{D}$ is a triplet $\mathcal{D}=(\mathcal{T},\mathcal{E},\mathcal{P})$ , where

$\mathcal{T}$ are the cells s.t. $\stackrel{\u0305}{\Omega}={\cup}_{K\in \mathcal{T}}\stackrel{\u0305}{K}$ . For each cell $K\in \mathcal{T}$ , $\leftK\right$ > 0 denotes the cell volume and $\partial K\stackrel{\scriptscriptstyle\mathrm{def}}{=}\stackrel{\u0305}{K}\backslash K$ its boundary. ${h}_{\mathcal{D}}\stackrel{\scriptscriptstyle\mathrm{def}}{=}{\mathrm{sup}}_{K\in \mathcal{T}}\mathrm{diam}\left(K\right)$ defines the discretization length.

$\mathcal{E}$ are the faces such that each face σ is included in a hyperplane of ${\mathbb{R}}^{d}$ with (d − 1)dimensional measure σ > 0. For each cell $K\in \mathcal{T},\hspace{0.5em}{\mathcal{E}}_{K}$ is the subset of $\mathcal{E}$ such that $\partial K={\cup}_{\sigma {\in \mathcal{E}}_{K}}\sigma $ . ${\mathcal{T}}_{\sigma}\{K\in \mathcal{T}\sigma \in {\mathcal{E}}_{K}\}$ contains the cells that share the face σ; the sets of inner and boundary faces are denoted by ${\mathcal{E}}_{\mathrm{int}}$ and ${\mathcal{E}}_{\mathrm{ext}}$ , respectively.

$\mathcal{P}=\{{\mathbf{x}}_{K}{\}}_{K\in \mathcal{T}}$ are the cell centers (not required to be the barycenters) s.t. x _{ K } ∈ K and K is starshaped with respect to x _{ K } . For all $K\in \mathcal{T}$ , $\sigma \in {\mathcal{E}}_{K},{d}_{K,\sigma}$ is given as the Euclidean distance between x _{ K } and σ.
For a more detailed definition of an admissible discretization see [3, 4, 21]. Let $\mathcal{M}$ be any set, then the cardinality of this set is indicated with ${n}_{\mathcal{M}}$. With this definition, the number of cells, i.e. the cardinality of the set $\mathcal{T}$, is given by ${n}_{\mathcal{T}}$. The face evaluation points are denoted as x _{ σ }, $\sigma \in \mathcal{E}$ (not required to be the barycenters). With ${\mathbf{\Lambda}}_{K}\stackrel{\scriptscriptstyle\mathrm{def}}{=}{\langle \mathbf{\Lambda}\rangle}_{K}\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{1}{\leftK\right}{\int}_{K}\mathbf{\Lambda}\mathrm{}\left(\mathbf{x}\right)\mathrm{d}x$ we define the averaged tensor on cell K, where the integral is meant componentwise. In addition, it is assumed that ${\mathbf{\Lambda}}_{\stackrel{\u0305}{K}}\in \left[{C}^{2}\right(\stackrel{\u0305}{K}){]}^{d\times d}$ for all $K\in \mathcal{T}$. Furthermore, n _{ K,σ } denotes the unit vector that is normal to σ and outward to K.
2.1 Abstract discretization framework
Integration of equation (1) over the control volume $K\in \mathcal{T}$ yields$${\int}_{\partial K}^{}\left(\mathbf{\Lambda}\nabla u\right)\xb7{\mathbf{n}}_{K}\mathrm{d}x={\int}_{K}^{}f\mathrm{d}x.$$(2)
Using Definition 1, the integral on the left of equation (2) can be spitted into integrals over the faces σ, which results in$$\sum _{\sigma \in {\mathcal{E}}_{K}}\mathrm{}{\int}_{\sigma}\mathrm{}\left(\mathbf{\Lambda}\nabla u\right)\xb7{\mathbf{n}}_{K,\sigma}\mathrm{d}x={\int}_{K}\mathrm{}f\mathrm{d}x.$$(3)
To discretize equation (3) we need to define flux approximations F _{ K,σ }, for each cell K and face $\sigma \in {\mathcal{E}}_{K}$, such that$${F}_{K,\sigma}\approx {\int}_{\sigma}\mathrm{}\left(\mathbf{\Lambda}\nabla u\right)\xb7{\mathbf{n}}_{K,\sigma}\mathrm{d}x,$$(4)
i.e. F _{ K,σ } is an approximation of the exact flux.
For this purpose, we define the following discrete solution spaces$${X}_{\mathcal{T}}\stackrel{\scriptscriptstyle\mathrm{def}}{=}\{v=({v}_{K}{)}_{K\in \mathcal{T}},{v}_{K}\in \mathbb{R}\},$$(5) $${X}_{\mathcal{E}}\stackrel{\scriptscriptstyle\mathrm{def}}{=}\{v=({v}_{\sigma}{)}_{\sigma \in \mathcal{E}},{v}_{\sigma}\in \mathbb{R}\},$$(6) $${X}_{\mathcal{D}}\stackrel{\scriptscriptstyle\mathrm{def}}{=}{X}_{\mathcal{T}}\times {X}_{\mathcal{E}}.$$(7)
For any element $v\in {X}_{\mathcal{D}}$, the cell unknowns are denoted by ${v}_{\mathcal{T}}\in {X}_{\mathcal{T}}$ and the face unknowns are denoted by ${v}_{\mathcal{E}}\in {X}_{\mathcal{E}}$. The space that takes into account the homogeneous zero Dirichlet boundary conditions is accordingly given as$${X}_{\mathcal{D},0}\stackrel{\scriptscriptstyle\mathrm{def}}{=}\{v\in {X}_{\mathcal{D}}{v}_{\sigma}=0,\hspace{1em}\forall \sigma \in {\mathcal{E}}_{\mathrm{ext}}\}.$$(8)
Additionally, the space of piecewise constant functions on $\mathcal{T}$ is defined as$${H}_{\mathcal{T}}\stackrel{\scriptscriptstyle\mathrm{def}}{=}\{v\in {L}^{2}(\mathrm{\Omega})\forall K\in \mathcal{T},\forall \mathbf{x}\in K,{v}_{K}\left(\mathbf{x}\right)={v}_{K},{v}_{K}\in \mathbb{R}\}.$$
Therefore, for all $v\in {H}_{\mathcal{T}}$ and for all $K\in \mathcal{T},{v}_{K}$ will denote the constant value of v on K. Furthermore, the extraction operator is defined for all $v\in {X}_{\mathcal{D}}$ as$${\mathrm{\Pi}}_{\mathcal{T}}:{X}_{\mathcal{D}}\mathrm{}\mapsto {H}_{\mathcal{T}}\mathrm{s.t}.,\forall K\in \mathcal{T},{\left({\mathrm{\Pi}}_{\mathcal{T}}v\right)}_{K}={v}_{K}.$$(9)
The presented finitevolume schemes differ in the choice of the fluxes ${F}_{K,\sigma}$ and the choice of the solution space ${X}_{h}\subseteq {X}_{\mathcal{D},0}$. Therefore, for all $K\in \mathcal{T}$, and for all $\sigma \in {\mathcal{E}}_{K}$, let ${F}_{K,\sigma}:{X}_{h}\mathrm{}\mapsto \mathbb{}\mathbb{R}$ be a numerical flux function that approximates the flux flowing out of K through σ such that the finitevolume scheme reads: Find u ∈ X _{ h } s.t.$$\sum _{\sigma \in {\mathcal{E}}_{K}}\mathrm{}{F}_{K,\sigma}\left(u\right)={\int}_{K}\mathrm{}f\mathrm{d}x,\hspace{1em}\forall K\in \mathcal{T},$$(10) $${F}_{K,\sigma}\left(u\right)+{F}_{L,\sigma}\left(u\right)=0,\hspace{1em}\forall \sigma \in {\mathcal{E}}_{\mathrm{int}},{\mathcal{T}}_{\sigma}=\left\{K,L\right\}.$$(11)
Cellcentered schemes:
For cellcentered schemes the solution space is given as$${X}_{h}=\{v=({v}_{\mathcal{T}},{v}_{\mathcal{E}})\in {X}_{\mathcal{D},0}{v}_{\sigma}={I}_{\sigma}{v}_{\mathcal{T}},\hspace{1em}\forall \sigma \in {\mathcal{E}}_{\mathrm{int}}\},$$(12)whereby, for each face σ, ${I}_{\sigma}\in {\mathcal{L}(X}_{\mathcal{T}};\mathbb{R})$ is a trace reconstruction operator. Here, ${\mathcal{L}(X}_{\mathcal{T}};\mathbb{R})$ is the space of linear operators that map some element $u\in {X}_{\mathcal{T}}$ to a constant value associated with the corresponding face. This means that for cellcentered schemes, the face unknowns are eliminated using some trace reconstruction operators $I=\{{I}_{\sigma}{\}}_{\sigma \in \mathcal{E}}$.
Here, only locally massconservative cellcentered schemes are considered such that for all $u\in {X}_{h}$, $\sigma \in {\mathcal{E}}_{\mathrm{int}}$ with ${\mathcal{T}}_{\sigma}=\left\{K,L\right\}$ $${F}_{K,\sigma}\left(u\right)+{F}_{L,\sigma}\left(u\right)=0.$$(13)
This means that equation (11) is fulfilled by construction and holds not only for the discrete solution.
Hybrid Mixed Mimetic (HMM) schemes:
For HMM schemes (see [14] for a detailed description of such schemes) the solution space is given as$${X}_{h}={X}_{\mathcal{D},0}.$$(14)
This means that HMM schemes introduce additional face unknowns, in contrast to cellcentered schemes where these unknowns are eliminated.
For the locally massconservative cellcentered schemes with solution space (12) and for the HMM schemes with solution space (14), it can be shown that problem (10)–(11) is equivalent to the problem: Find $u\in {X}_{h}$ s.t.$${a}_{\mathcal{D}}\left(u,v\right)={\int}_{\mathrm{\Omega}}\mathrm{}f{\mathrm{\Pi}}_{\mathcal{T}}v\mathrm{d}x,\hspace{1em}\forall v\in {X}_{h},$$(15)with the form, for all (u, v) ∈ [X _{ h }]^{2},$${a}_{\mathcal{D}}\left(u,v\right)\sum _{K\in \mathcal{T}}\mathrm{}\sum _{\sigma \in {\mathcal{E}}_{K}}\mathrm{}{F}_{K,\sigma}\left(u\right)\left({v}_{\sigma}{v}_{K}\right).$$(16)
Therefore, each solution of problem (10) and (11) also solves problem (15). Starting from the discrete problem (15), with ${a}_{\mathcal{D}}$ defined by (16), equation (11) is fulfilled for the cellcentered schemes by construction, whereas for the HMM schemes it is obtained by taking ${v}_{\sigma}=1,{v}_{{\sigma}^{\text{'}}}=0$ for all σ′ ≠ σ, and v _{ K } = 0 for all K $.$
Equation (10) is obtained by taking, for each $K\in \mathcal{T},{v}_{K}=1$ and ${v}_{{K}^{\text{'}}}=0$ for all ${K}^{\text{'}}\ne K$ and by using the flux continuity (11). This shows that the formulation (10) and (11) is equivalent to (15).
Remark 2. If the flux functions are linear, then the form (16) is linear in both arguments. However, in Section 2.3 we will also consider nonlinear flux functions such that the linearity is only given for the second argument.
The advantage of the formulation (15) is that properties like coercivity can easily be defined for the form ${a}_{\mathcal{D}}$. Such properties are briefly defined in the next section.
2.2 Properties of discretization schemes
Besides simplicity, parallelizability, computational efficiency, flexibility, and code integrability, mathematical and physical properties of discretization schemes, e.g. consistency, coercivity or monotonicity, can be defined. Unfortunately, to the author’s knowledge, there is currently no scheme that satisfies all these properties. Therefore, appropriate schemes have to be chosen application dependent. In general, the convergence of schemes can be proven if the scheme is consistent and coercive, as it has been done in [4, 8, 12, 21]. In the following, we briefly describe such fundamental properties for cellcentered schemes (with discrete solution space (12)) and for HMM schemes (with discrete solution space (14)).
Consistency
There is no unique discretizationindependent definition of consistency. For example, finitevolume schemes are in general not consistent in the finitedifference context [27]. In general, a scheme is consistent if the truncation error between discrete and continuous operators goes to zero for ${h}_{\mathcal{D}}\to 0$. In the context of finitevolume schemes we say that a scheme is consistent if the numerical flux approximates the exact flux for regular functions, meaning that$${F}_{K,\sigma}\left({\phi}_{\mathcal{D}}\right)={\stackrel{\u0305}{F}}_{K,\sigma}\left(\phi \right)+\mathcal{O}\left(\left\sigma \right\mathrm{diam}\left(K\right)\right),\hspace{1em}\forall \phi \in \mathcal{D},$$(17)where $\mathcal{D}\subset {C}_{0}\left(\stackrel{\u0305}{\mathrm{\Omega}}\right)$ is a test function space which is assumed to be dense in ${H}_{0}^{1}\left(\mathrm{\Omega}\right)$, and ${F}_{K,\sigma},{\stackrel{\u0305}{F}}_{K,\sigma}$ are the discrete and exact flux functions, respectively. Furthermore, ${\phi}_{\mathcal{D}}=({\phi}_{\mathcal{T}},{\phi}_{\mathcal{E}})\in {X}_{D}$ is defined as $({\phi}_{\mathcal{T}}{)}_{K}=\phi ({\mathbf{x}}_{K})$ for all $K\in \mathcal{T}$, and $({\phi}_{\mathcal{E}}{)}_{\sigma}=\phi ({\mathbf{x}}_{\sigma})$ for all $\sigma \in \mathcal{E}$. More details can for example be found in [21].
Coercivity
A scheme is called uniformly coercive if the form ${a}_{\mathcal{D}}$, defined in (16), satisfies the following estimate:$${a}_{\mathcal{D}}\left(u,u\right)\ge \gamma {\left\rightu\left\right}_{{X}_{h}}^{2},\hspace{1em}\forall u\in {X}_{h},$$(18)with some appropriate discrete norm ${\left\right\xb7\left\right}_{{X}_{h}}$, and γ > 0 that is independent of u and ${h}_{\mathcal{D}}$. When talking about coercivity in the following, we always refer to the uniform coercivity property.
Minimum and maximum principles
It is desirable that the discrete solution satisfies properties of the exact solution. An important property from a physical point of view is the minimum and maximum principle. Schemes that satisfy these principles prevent oscillations of the discrete solutions, such that the discrete solution remains within physical bounds. Such extremum principles can be found in [28]. For the definition of discrete minimum and maximum principles, we refer to the overview provided by [29]. Furthermore, a scheme is said to be monotone if the discretization matrix $\mathbb{A}$ is monotone, which means that all entries of its inverse are nonnegative i.e. ${\mathbb{A}}^{1}\ge 0$. The monotonicity property is for example satisfied by Mmatrices. Mmatrices are monotone matrices with nonpositive offdiagonal entries (see, e.g. [30]).
Sparsity
For solving largescale systems, it is indispensable that the discretization matrices are sparse, which means that the number of entries (noe) in the matrix is significantly smaller than ${n}_{\mathcal{T}}^{2}$. The stencil denotes the noe of each row. Therefore, schemes resulting in small stencils are important for largescale problems.
2.3 Cellcentered schemes
Within this section, we present different cellcentered finitevolume schemes. Here, the discrete solution space is given as (12), and the face values are eliminated using some trace reconstruction operators $\{{I}_{\sigma}{\}}_{\sigma \in \mathcal{E}}$. As already mentioned, only locally massconservative schemes are considered here such that equation (11) holds by construction.
Family of weighted schemes
An established idea to obtain monotone or extremumprinciplespreserving schemes, as those developed in [17–20, 22, 23, 31, 32], is to combine, for each interior edge $\sigma {\mathrm{}\in \mathcal{E}}_{\mathrm{int}}$, with ${\mathcal{T}}_{\sigma}=\{K,L\}$, two consistent linear flux approximations ${\stackrel{\u0303}{F}}_{K,\sigma}\left(u\right)$ and ${\stackrel{\u0303}{F}}_{L,\sigma}\left(u\right)$ as a convex combination, using solutiondependent weights. Thus, the final flux F _{ K,σ }(u) is given as$$\begin{array}{cc}& {F}_{K,\sigma}\left(u\right)={\mu}_{K,\sigma}\left(u\right){\stackrel{\u0303}{F}}_{K,\sigma}\left(u\right){\mu}_{L,\sigma}\left(u\right){\stackrel{\u0303}{F}}_{L,\sigma}\left(u\right),\\ \mathrm{with}& {\mu}_{K,\sigma}\left(u\right)\ge 0,{\mu}_{L,\sigma}\left(u\right)\ge 0,\\ \mathrm{and}& {\mu}_{K,\sigma}\left(u\right)+{\mu}_{L,\sigma}\left(u\right)=1.\end{array}$$(19)
For any $K\in \mathcal{T}$ and $\sigma \in {\mathcal{E}}_{K}\cap {\mathcal{E}}_{\mathrm{int}}$, the linear flux ${\stackrel{\u0303}{F}}_{K,\sigma}\left(u\right)$ is built in order to ensure a strong consistency property, such that the total flux ${F}_{K,\sigma}$ satisfies a weak form of consistency, see [21] for detailed information. A general definition of these subfluxes is$${\stackrel{\u0303}{F}}_{K,\sigma}\left(u\right)\stackrel{\scriptscriptstyle\mathrm{def}}{=}\left\sigma \right\sum _{{\sigma}^{\text{'}}\in {S}_{K,\sigma}}^{}{\alpha}_{K,\sigma {\sigma}^{\text{'}}}\left({I}_{{\sigma}^{\text{'}}}u{u}_{K}\right),$$(20)with face stencil ${S}_{K,\sigma}$ and trace reconstruction operator ${I}_{{\sigma}^{\text{'}}}$.
The coefficients μ _{ K,σ } and μ _{ L,σ } are chosen to improve some of the previously mentioned properties, e.g. monotonicity. These coefficients may nonlinearly depend on the unknown u. Therefore, the form ${a}_{\mathcal{D}}$ is also nonlinear in its first argument. This can be prevented by introducing an additional argument for the fluxes and correspondingly for the form ${a}_{\mathcal{D}}$. The details are not presented here, for a detailed description see [21].
The final flux function (19) is locally massconservative, whereas the subfluxes $\stackrel{\u0303}{F}$ are in general not massconservative. In the following, different schemes that fit into this family are shortly introduced.
AvgMPFA scheme
The easiest choice of coefficients is μ _{ K,σ } = μ _{ L,σ } = 0.5, which results in a linear finitevolume scheme that is in the following denoted as AvgMPFA scheme.
NLTPFA schemes
To derive a Nonlinear TwoPoint Flux Approximation (NLTPFA), the different terms are reordered such that the flux is written as$${F}_{K,\sigma}\left(u\right)={t}_{L,\sigma}\left(u\right){u}_{L}{t}_{K,\sigma}\left(u\right){u}_{K}\underset{{R}_{K,\sigma}\left(u\right)}{\mathrm{}\underbrace{\left({\mu}_{L,\sigma}\left(u\right){\lambda}_{L,\sigma}\left(u\right){\mu}_{K,\sigma}\left(u\right){\lambda}_{K,\sigma}\left(u\right)\right)}}.$$(21)
The idea of the NLTPFA scheme is to choose the weights such that R _{ K,σ }(u) = 0. From a numerical point of view, it is sufficient that R _{ K,σ }(u) ≤ ϵ. This is given, under the assumption that λ _{ K,σ } λ _{ L,σ } ≥ 0, for the choice$$\begin{array}{l}{\mu}_{K,\sigma}\left(u\right)=\frac{\left{\lambda}_{L,\sigma}\left(u\right)\right\mathrm{}+\mathrm{}\u03f5}{\left{\lambda}_{K,\sigma}\left(u\right)\right\mathrm{}+\mathrm{}\left{\lambda}_{L,\sigma}\left(u\right)\right\mathrm{}+\mathrm{}2\u03f5},\\ {\mu}_{L,\sigma}\left(u\right)=\frac{\left{\lambda}_{K,\sigma}\left(u\right)\right\mathrm{}+\mathrm{}\u03f5}{\left{\lambda}_{K,\sigma}\left(u\right)\right\mathrm{}+\mathrm{}\left{\lambda}_{L,\sigma}\left(u\right)\right\mathrm{}+\mathrm{}2\u03f5}.\end{array}$$(22)
Further details and the discussion of the case λ _{ K,σ } λ _{ L,σ } < 0 can be found in [21], where the monotonicity of the NLTPFA scheme is also discussed. Here, ϵ is set to 10^{−8} but any other small number could also be chosen.
NLMPFA scheme
A Nonlinear MultiPoint Flux Approximation (NLMPFA) is derived by splitting the subfluxes into a twopoint part and a residual flux part$$\begin{array}{l}{\stackrel{\u0303}{F}}_{K,\sigma}\left(u\right)=\left\sigma \right{\beta}_{\sigma}\left({u}_{L}{u}_{K}\right)+{\stackrel{\u0303}{F}}_{K,\sigma}^{\mathrm{res}}\left(u\right),\\ {\stackrel{\u0303}{F}}_{L,\sigma}\left(u\right)=\left\sigma \right{\beta}_{\sigma}\left({u}_{K}{u}_{L}\right)+{\stackrel{\u0303}{F}}_{L,\sigma}^{\mathrm{res}}\left(u\right),\end{array}$$(23)with some appropriate ${\beta}_{\sigma}$ and ${\stackrel{\u0303}{F}}_{K,\sigma}^{\mathrm{res}},{\stackrel{\u0303}{F}}_{L,\sigma}^{\mathrm{res}}$. Here, the weights are chosen as$$\begin{array}{l}{\mu}_{K,\sigma}\left(u\right)=\frac{\left{\stackrel{\u0303}{F}}_{L,\sigma}^{\mathrm{res}}\left(u\right)\right\mathrm{}+\mathrm{}\u03f5}{\left{\stackrel{\u0303}{F}}_{K,\sigma}^{\mathrm{res}}\left(u\right)\right\mathrm{}+\mathrm{}\left{\stackrel{\u0303}{F}}_{L,\sigma}^{\mathrm{res}}\left(u\right)\right+\mathrm{}2\u03f5},\\ {\mu}_{L,\sigma}\left(u\right)=\frac{\left{\stackrel{\u0303}{F}}_{K,\sigma}^{\mathrm{res}}\left(u\right)\right\mathrm{}+\mathrm{}\epsilon}{\left{\stackrel{\u0303}{F}}_{K,\sigma}^{\mathrm{res}}\left(u\right)\right\mathrm{}+\mathrm{}\left{\stackrel{\u0303}{F}}_{L,\sigma}^{\mathrm{res}}\left(u\right)\right\mathrm{}+\mathrm{}2\u03f5}.\end{array}$$(24)
Again, ϵ is added due to numerical reasons. The final flux is then given by (19) with weights (24). It can be shown that this scheme satisfies discrete extremum principles.
Remark 3. Without going into detail, we shortly summarize the main differences between the NLTPFA and the NLMPFA scheme. For a detailed description we refer to [21]. Both schemes mainly differ in the choice of the weights μ _{ K,σ } and μ _{ L,σ } (the coefficients ${\alpha}_{K,\sigma {\sigma}^{\text{'}}}$ in the subfluxes (20) may also vary). The weights of the NLTPFA scheme are defined by ${\lambda}_{K,\sigma}$ and ${\lambda}_{L,\sigma}$ which can be written as a sum of solution values, i.e. λ _{ K,σ } = ∑ω _{ M } u _{ M }. Whereas, those of the NLMPFA scheme are defined by the residual fluxes ${\stackrel{\u0303}{F}}_{K,\sigma}^{\mathrm{res}}$ and ${\stackrel{\u0303}{F}}_{L,\sigma}^{\mathrm{res}}$ which can be written as a sum of solution value differences, i.e. ${\stackrel{\u0303}{F}}_{K,\sigma}^{\mathrm{res}}=\sum \mathrm{}{\omega}_{M}({u}_{M}{u}_{K})$.
Linear TPFA/MPFA schemes
Here, we shortly demonstrate that wellestablished schemes such as the TwoPoint Flux Approximation (TPFA) or the MultiPoint Flux Approximation (MPFA) also fit into this family of schemes. The difference here is that the subfluxes are constructed to satisfy flux continuity, i.e. ${\stackrel{\u0303}{F}}_{K,\sigma}+{\stackrel{\u0303}{F}}_{L,\sigma}=0$. Therefore, the final flux is given as ${F}_{K,\sigma}={\stackrel{\u0303}{F}}_{K,\sigma}$ (20), independent of the weights μ _{ K,σ }, μ _{ L,σ }.
For the TPFA scheme it holds that ${\mathrm{S}}_{K,\sigma}=\left\{\sigma \right\}$ and the coefficients are$${\alpha}_{K,\sigma \sigma}=\frac{{\mathbf{n}}_{K,\sigma}^{T}{\mathbf{\Lambda}}_{\mathbf{K}}\left({\mathbf{x}}_{\sigma}{\mathbf{x}}_{K}\right)}{{\left\left{\mathbf{x}}_{\sigma}{\mathbf{x}}_{K}\right\right}_{2}^{2}},{\alpha}_{L,\sigma \sigma}=\frac{{\mathbf{n}}_{L,\sigma}^{T}{\mathbf{\Lambda}}_{L}\left({\mathbf{x}}_{\sigma}{\mathbf{x}}_{L}\right)}{{\left\left{\mathbf{x}}_{\sigma}{\mathbf{x}}_{L}\right\right}_{2}^{2}}.$$(25)
The trace reconstruction operator is$${I}_{\sigma}u=\frac{{\alpha}_{K,\sigma \sigma}{u}_{K}+{\alpha}_{L,\sigma \sigma}{u}_{L}}{{\alpha}_{K,\sigma \sigma}+{\alpha}_{L,\sigma \sigma}},$$(26)such that the final flux reads$${F}_{K,\sigma}\left(u\right)=\left\sigma \right\frac{{\alpha}_{K,\sigma \sigma}{\alpha}_{L,\sigma \sigma}}{{\alpha}_{K,\sigma \sigma}+{\alpha}_{L,\sigma \sigma}}\left({u}_{L}{u}_{K}\right).$$(27)
This can be similarly shown for MPFA schemes but the details are not presented here.
2.4 Hybrid/mimetic schemes
Within this section, we present hybrid and mimetic schemes. Here, the discrete solution space is given as (14), such that additional face unknowns are introduced. Here, the fluxes are given, for all $u\in {X}_{h},K\in \mathcal{T},\sigma \in {\mathcal{E}}_{K}$, as$${F}_{K,\sigma}\left(u\right)\left\sigma \right\sum _{{\sigma}^{\mathrm{\prime}}\in {\mathcal{E}}_{K}}\mathrm{}{\alpha}_{K,\sigma {\sigma}^{\mathrm{\prime}}}\left({u}_{{\sigma}^{\mathrm{\prime}}}{u}_{K}\right),$$(28)with the face unknowns ${u}_{{\sigma}^{\mathrm{\prime}}}$, the cell unknowns u _{ K }, and coefficients ${\alpha}_{K,\sigma {\sigma}^{\mathrm{\prime}}}$. These coefficients are chosen such that the resulting scheme is coercive and consistent.
Defining the matrices,$${\mathbf{u}}_{{\mathcal{E}}_{K}}\stackrel{\scriptscriptstyle\mathrm{def}}{=}({u}_{\sigma}{)}_{\sigma \in {\mathcal{E}}_{K}}=\left(\begin{array}{c}{u}_{{\sigma}_{1}}\\ \vdots \\ {u}_{{\sigma}_{{n}_{{\mathcal{E}}_{K}}}}\end{array}\right),{\hspace{1em}\mathbb{C}}_{K}=\mathrm{diag}({\sigma}_{1},\cdots ,{\sigma}_{{n}_{{\mathcal{E}}_{K}}}\left\right),$$(29)we obtain the cell flux vector as$${\mathbf{F}}_{K,{\mathcal{E}}_{K}}\stackrel{\scriptscriptstyle\mathrm{def}}{=}\left(\begin{array}{c}{F}_{K,{\sigma}_{1}}\\ \vdots \\ {F}_{K,{\sigma}_{{n}_{{\mathcal{E}}_{K}}}}\end{array}\right)={\mathbb{C}}_{K}{\mathbb{W}}_{K}{\mathbb{C}}_{K}\left({\mathbf{u}}_{{\mathcal{E}}_{K}}{u}_{K}\mathbf{e}\right),$$(30)with ${\mathbb{W}}_{K}={\mathbb{A}}_{K}{\mathbb{C}}_{K}^{1}$, ${\mathbb{A}}_{K,\sigma {\sigma}^{\mathrm{\prime}}}={\alpha}_{K,\sigma {\sigma}^{\mathrm{\prime}}}$, and e is the vector with entries equal to one. Additionally, let ${\mathbb{N}}_{K}$ be the conormal matrix, and let ${\mathbb{R}}_{K}$ be the matrix that contains the scaled distance vectors:$${\mathbb{N}}_{K}=\left(\begin{array}{c}{\mathbf{n}}_{K,{\sigma}_{1}}^{T}{\mathbf{\Lambda}}_{K}\\ \vdots \\ {\mathbf{n}}_{K,{\sigma}_{{n}_{{\mathcal{E}}_{K}}}}^{T}{\mathbf{\Lambda}}_{K}\end{array}\right),\hspace{0.5em}{\mathbb{R}}_{K}=\left(\begin{array}{c}\left{\sigma}_{1}\right{({\mathbf{x}}_{{\sigma}_{1}}{\mathbf{x}}_{K})}^{T}\\ \vdots \\ \left{\sigma}_{{n}_{{\mathcal{E}}_{K}}}\right{({\mathbf{x}}_{{\sigma}_{{n}_{{\mathcal{E}}_{K}}}}{\mathbf{x}}_{K})}^{T}\end{array}\right).$$(31)
Then, the consistency condition is given by$${\mathbb{N}}_{K}={\mathbb{W}}_{K}{\mathbb{R}}_{K}.$$(32)
From this consistency condition, we can derive a general form of the coefficient matrix ${\mathbb{W}}_{K}$:$${\mathbb{W}}_{K}={\mathbb{N}}_{K}({\mathbb{N}}_{K}^{T}{\mathbb{R}}_{K}{)}^{1}{\mathbb{N}}_{K}^{T}+{\mathbb{S}}_{K},$$(33)with a stabilization matrix ${\mathbb{S}}_{K}$ such that ${\mathbb{S}}_{K}{\mathbb{R}}_{K}=\mathbb{O}$. Choosing x _{ σ } as the barycenter of face σ and using the geometrical identity$$\sum _{\sigma \in {\mathcal{E}}_{K}}^{}{\left\sigma \right\mathbf{n}}_{K,\sigma}({\mathbf{x}}_{\sigma}{\mathbf{x}}_{K}{)}^{T}=K\mathbb{I},\hspace{1em}\forall \mathit{K}\in \mathcal{T},$$(34)it can be seen that$${\mathbb{N}}_{K}^{T}{\mathbb{R}}_{K}={\leftK\right\mathbf{\Lambda}}_{K}.$$(35)
Therefore, the matrix ${\mathbb{W}}_{K}$ is symmetric if the stabilization matrix ${\mathbb{S}}_{K}$ is chosen symmetric.
MFD scheme
A simple choice of ${\mathbb{S}}_{K}$ is$${\mathbb{S}}_{K}={\nu}_{K}\left(\mathbb{I}{\mathbb{R}}_{K}({\mathbb{R}}_{K}^{T}{\mathbb{R}}_{K}{)}^{1}{\mathbb{R}}_{K}^{T}\right),$$(36)with ${\nu}_{K}=\frac{1}{2}\mathrm{trace}\left({\mathbb{N}}_{K}\right({\mathbb{N}}_{K}^{T}{\mathbb{R}}_{K}{)}^{1}{\mathbb{N}}_{K}^{T})$. This results in the final matrix$${\mathbb{W}}_{K}={\mathbb{N}}_{K}({\mathbb{N}}_{K}^{T}{\mathbb{R}}_{K}{)}^{1}{\mathbb{N}}_{K}^{T}+{\nu}_{K}\left(\mathbb{I}{\mathbb{R}}_{K}({\mathbb{R}}_{K}^{T}{\mathbb{R}}_{K}{)}^{1}{\mathbb{R}}_{K}^{T}\right).$$(37)
For this matrix, the consistency condition (32) is fulfilled by construction, whereas the coercivity can be proven by using some kind of stability condition, see [33].
HFV scheme
Another interesting scheme that fits into this framework, and also fits into the recently developed gradient discretization framework [15], is the Hybrid FiniteVolume (HFV) scheme [8]. The main idea of this scheme is the construction of a consistent discrete gradient, which at the same time defines the form ${a}_{\mathcal{D}}$, such that the consistency is naturally fulfilled. In the following we shortly introduce the scheme and present the main ideas.
Let us recall the weak formulation of problem (1): Find $\stackrel{\u0305}{u}\in {H}_{0}^{1}\left(\mathrm{\Omega}\right)$ such that$${\int}_{\mathrm{\Omega}}\mathrm{}\mathbf{\Lambda}\nabla \stackrel{\u0305}{u}\xb7\nabla v\mathrm{d}x={\int}_{\mathrm{\Omega}}\mathrm{}\mathrm{fv}\mathrm{d}x,\hspace{1em}\forall v\in {H}_{0}^{1}\left(\mathrm{\Omega}\right).$$
The idea of symmetric gradient discretization schemes is to replace the operator ∇ with a consistent approximation ${\widehat{\nabla}}_{\mathcal{D}}$ such that the form of the discrete problem (15) is given by$${a}_{\mathcal{D}}\left(u,v\right){\int}_{\mathrm{\Omega}}^{}\mathbf{\Lambda}{\widehat{\nabla}}_{\mathcal{D}}u\xb7{\widehat{\nabla}}_{\mathcal{D}}v\mathrm{d}x.$$(38)
Therefore, the scheme is defined by the discrete gradient reconstruction operator ${\widehat{\nabla}}_{\mathcal{D}}$. First, let us define a discrete gradient on each cell $K\in \mathcal{T}$ as$${\widehat{\nabla}}_{K}u=\frac{1}{\leftK\right}\sum _{\sigma \in {\mathcal{E}}_{K}}\mathrm{}\left\sigma \right({u}_{\sigma}{u}_{K}){\mathbf{n}}_{K,\sigma},$$(39)the consistency of this formula follows thanks to (34). However, to end up in a coercive form, we need an additional stabilization term, which is defined as$$\begin{array}{l}{\mathbf{S}}_{K,\sigma}u={{S}_{K,\sigma}\mathbf{n}}_{K,\sigma},\\ \mathrm{with}{S}_{K,\sigma}\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{\sqrt{d}}{{d}_{K,\sigma}}\left({u}_{\sigma}{u}_{K}{\widehat{\nabla}}_{K}u\xb7\left({\mathbf{x}}_{\sigma}{\mathbf{x}}_{K}\right)\right).\end{array}$$(40)
This results in the final discrete gradient: For all $K\in \mathcal{T},\sigma \in {\mathcal{E}}_{K}$ $${\widehat{\nabla}}_{\mathcal{D}}u\left(\mathbf{x}\right)\stackrel{\scriptscriptstyle\mathrm{def}}{=}{\widehat{\nabla}}_{K}u+{{S}_{K,\sigma}\mathbf{n}}_{K,\sigma},\mathrm{for}\mathrm{a.e}.\hspace{0.5em}\mathbf{x}\in {\Delta}_{K,\sigma},$$(41)where ${\mathbf{\Delta}}_{K,\sigma}$ is the convex hull that includes the face σ and the cell center x _{ K }. Inserting these discrete gradients into the form (38) and reordering of terms lead to the form that is defined in (16). Thus, the fluxes and the coefficient matrix ${\mathbb{W}}_{K}$ (30) can also be identified. Therefore, the hybrid finitevolume scheme also belongs to the family of hybrid mimetic schemes. The coefficient matrix, a detailed description, and the proof of consistency and coercivity can be found in [8].
Remark 4. For a more detailed summary of finitevolume schemes we refer to [34, 35].
3 Numerical experiments
In this chapter, the behavior of the different finitevolume schemes, that have been presented in the last section, is investigated. The AvgMPFA, NLTPFA, and NLMPFA schemes need the coefficients ${\alpha}_{K,\sigma {\sigma}^{\mathrm{\prime}}}$, that occur in the subfluxes (20). These coefficients are calculated using optimization strategies together with a PrimalDual Simplex Method provided by the opensource library GNU Linear Programming Kit ^{1} (GLPK). Detailed explanations can be found in [21, 24]. Furthermore, the harmonic averaging interpolator [36] is used as reconstruction operator I _{ σ }, which is needed to define the subfluxes (20).
The mimetic finitedifference scheme with fluxes (30) and local cell matrix (37) is denoted as MFD scheme. The hybrid finitevolume scheme defined by the discrete gradients (41) is denoted as HFV scheme. Finally, with TPFA we denote the scheme with fluxes (27) and with MPFAO the scheme that is described in [37]. The Box method [38, 39] is a vertexcentered finitevolume scheme that uses finiteelement basis functions on each cell to calculate fluxes over subcontrol volume faces.
All simulations are performed using the opensource simulator DuMu^{x} [40], which comes in the form of an additional DUNE module [41]. Newton’s method is used to solve the nonlinear systems of equations occurring for the nonlinear finitevolume schemes, where the iteration loop is stopped if the absolute residual is below 10^{−5} with respect to the Euclidean norm.
In this section, more general boundary conditions are considered, where the Dirichlet conditions are taken into account using the function $g\in {H}^{\frac{1}{2}}\left(\partial \mathrm{\Omega}\right)$ such that the weak solution $\stackrel{\u0305}{u}\in {H}^{1}\left(\mathrm{\Omega}\right)$ has to satisfy $T\stackrel{\u0305}{u}=g$ (where $T:{H}^{1}\left(\mathrm{\Omega}\right)\mapsto {L}^{2}\left(\partial \mathrm{\Omega}\right)$ denotes the trace operator).
A discrete seminorm on the space ${X}_{\mathcal{D}}$ is given by$$\left\rightu{}_{\mathcal{D}}\stackrel{\scriptscriptstyle\mathrm{def}}{=}{\left(\sum _{K\in \mathcal{T}}\mathrm{}\left(\sum _{\sigma \in {\mathcal{E}}_{K}}\mathrm{}\frac{\left\sigma \right}{{d}_{K,\sigma}}({v}_{\sigma}{v}_{K}{)}^{2}\right)\right)}^{1/2},$$(42)which is a norm on the space ${X}_{\mathcal{D},0}$. For the cellcentered schemes, we define the following discrete norm (that takes into account the Dirichlet data)$$\left\rightu{}_{\mathcal{T}}\stackrel{\scriptscriptstyle\mathrm{def}}{=}{\left(\sum _{K\in \mathcal{T}}\mathrm{}\left(\sum _{\sigma \in {\mathcal{E}}_{K}\cap {\mathcal{E}}_{\mathrm{int}}}\mathrm{}\sigma {d}_{K,\sigma}{\left(\frac{{v}_{L}{v}_{K}}{{d}_{K,\sigma}+{d}_{L,\sigma}}\right)}^{2}+\sum _{\sigma \in {\mathcal{E}}_{K}\cap {\mathcal{E}}_{\mathrm{ext}}}\mathrm{}\frac{\left\sigma \right}{{d}_{K,\sigma}}({\langle g\rangle}_{\sigma}{v}_{K}{)}^{2}\right)\right)}^{1/2},$$(43)whereas 〈g〉_{ σ } denotes the average of g on the face σ. Using the CauchySchwarz inequality, one can show that (see [8] for more details)$$\left\rightu{}_{\mathcal{T}}\le \leftu\right{}_{\mathcal{D}},\hspace{0.5em}\forall u\in {X}_{\mathcal{D},g},$$(44)with$${X}_{\mathcal{D},g}\stackrel{\scriptscriptstyle\mathrm{def}}{=}\{v\in {X}_{\mathcal{D}}{v}_{\sigma}={\langle g\rangle}_{\sigma},\hspace{0.5em}\forall \sigma \in {\mathcal{E}}_{\mathrm{ext}}\}.$$(45)
Owing to inequality (44), the discrete norm (43) will also be used for the hybrid and mimetic schemes, although it is not a discrete norm for these schemes since it does not depend on the face unknowns. For measuring the coercivity of the schemes, the following estimate is defined$${e}_{\mathcal{D}}\left(u\right)\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{{a}_{\mathcal{D}}\left(u,u\right)}{{\left\leftu\right\right}_{\mathcal{T}}^{2}}.$$(46)
In the next section the coercivity estimate is evaluated for the numerical solution. Please note that this is not sufficient to show the coercivity of the schemes, but it serves as a good indicator.
In the following test cases, the properties that have been mentioned in Section 2.2 are investigated numerically. Hereby, the numerical and exact solutions are denoted as u _{ n } (where n indicates the grid refinement level) and $\stackrel{\u0305}{u}$, respectively.
Remark 5. In this section, hexahedral grids are considered (or quadrilateral grids for the case of a twodimensional domain). Such grids are, in general, not admissible in the sense of Definition 1, because the faces are usually nonplanar. Therefore, the averaged normal vectors are calculated which are still denoted by n _{ K,σ }. For more details, we refer to [24, 37, 42]. In [42], the authors suggest to introduce, for strongly curved faces, additional degrees of freedom accounting for tangential velocities. These additional velocities are related to the planar face which is defined through the averaged normal vector. Such additional degrees of freedom are not introduced here, because we did not observe any convergence rate reduction for the considered examples.
3.1 Convergence behavior
Within this section, the convergence rates of the different schemes are compared for threedimensional test problems with smooth solutions. The convergence of the family of schemes (19) has been proven, under the assumption of coercivity, in [21]. Here, we demonstrate that the convergence rates are similar to wellestablished schemes. These examples are based on our previous work [24], but here we are using different norms and a more challenging tensor for the highly anisotropic test case. In addition, the convergence behavior of the MPFAO, MFD, and HFV are investigated and compared.
The convergence behavior is investigated for the meshes shown in Figure 1. The checkerboard and the random hexahedral meshes are taken from the FVCA6 benchmark [43], whereas the nonconvex grid and the curved grid are generated with the Matlab Reservoir Simulation Toolbox (MRST) [44]. For all grids, we use the dunealugrid module [45]. Except for the checkerboard mesh, all grids exhibit nonplanar faces. Therefore, we calculate the integrated normal vectors, see [24]. For the MPFAO scheme this is applied to the subfaces of the control volumes, see [37].
Mildly anisotropic test case
The first test case is similar to test 1 from the FVCA6 benchmarks [43]. Here, the exact solution$$\stackrel{\u0305}{u}\left(x,y,z\right)=\mathrm{sin}\left(\pi x\right)\mathrm{sin}\left(\pi \left(y+\frac{1}{2}\right)\right)\mathrm{sin}\left(\pi \left(z+\frac{1}{3}\right)\right)+1,$$(47)is prescribed on the unit cube Ω = [0, 1]^{3}. Dirichlet conditions are set on the boundary as $g=\stackrel{\u0305}{u}$ on ∂Ω. The permeability tensor is$$\mathbf{\Lambda}=\left(\begin{array}{lll}1& 0.5& 0\\ 0.5& 1& 0.5\\ 0& 0.5& 1\end{array}\right).$$(48)
The discrete L ^{2}errors and H ^{1}errors are shown in Figure 2. Please note that the results of the AvgMPFA scheme are not shown since they only differ slightly from the NLTPFA results. Furthermore, the MPFAO scheme cannot be applied to hanging nodes on threedimensional grids, which is why there are no results of the MPFAO scheme for the checkerboard mesh.
Fig. 2 Logarithm of L ^{2}error (left) and discrete H ^{1}error (right) for convergence test case one. From top to bottom: checkerboard, random, nonconvex, curved mesh. 
It can be seen that for increasing mesh refinement, the schemes converge with second order accuracy in the L ^{2}norm and at least first order accuracy in the H ^{1}norm. The convergence rates of all schemes are quite similar. We also observe that the MFD scheme produces the smallest errors. The errors of the NLMPFA scheme are the largest ones for most grids and refinement levels. Additionally, the behavior of the MPFAO and HFV is similar, especially for the random and nonconvex grids. In this example, the NLTPFA scheme produces smaller errors than the NLMPFA scheme. Considering Table 1, we observe that the coercivity estimates ${e}_{\mathcal{D}}$ are bounded from below which indicates that all schemes are coercive for this test case on all grids. The NLTPFA scheme converges within two Newton iterations, whereas the Newton convergence of the NLMPFA is in general worse.
Highly anisotropic test case
The next test problem investigates a highly anisotropic permeability tensor$$\mathbf{\Lambda}=\frac{1}{{x}^{2}+{y}^{2}}\left(\begin{array}{lll}\beta {x}^{2}+{y}^{2}& \left(\beta 1\right)\mathrm{xy}& 0\\ \left(\beta 1\right)\mathrm{xy}& {x}^{2}+\beta {y}^{2}& 0\\ 0& 0& \left({z}^{2}+\beta \right)\left({x}^{2}+{y}^{2}\right)\end{array}\right),$$(49)with β = 10^{−2} and the exact solution$$\stackrel{\u0305}{u}\left(x,y,z\right)=\mathrm{sin}\left(2\pi x\right)\mathrm{sin}\left(2\pi y\right)\mathrm{sin}\left(2\pi z\right)+1.$$(50)
Again, Dirichlet conditions are set on the whole boundary corresponding to the exact solution.
Figure 3 shows the errors for test case two, with solution (50). Again, the results of the AvgMPFA scheme are not shown since they only differ slightly from the NLTPFA results. As in the previous example, no results could be obtained for the MPFAO scheme on the checkerboard mesh due to the nonconformity of the grid.
Fig. 3 Logarithm of L ^{2}error (left) and discrete H ^{1}error (right) for convergence test case one. From top to bottom: checkerboard, random, nonconvex, curved mesh. 
In Figure 3, we can see that the NLTPFA scheme produces now similar errors than the HFV and MPFAO scheme, and results in the smallest H ^{1}errors on the random grid. Considering Table 2, we observe that the coercivity estimates ${e}_{\mathcal{D}}$ are bounded from below which indicates that all schemes are coercive also for this test case. Again, the NLTPFA scheme converges within two Newton iterations for all grids. However, the NLMPFA needs much more Newton iterations for the checkerboard mesh. The worse convergence behavior is probably caused by the fact that the weights μ _{ K,σ }, μ _{ L,σ } are constructed differently for the NLMPFA and NLTPFA schemes, as explained in Remark 3. Sign changes for the functions λ _{ K,σ } are rare, whereas sign changes occur more often for the residual fluxes ${\stackrel{\u0303}{F}}_{K,\sigma}^{\mathrm{res}}$, which are used to define the weights of the NLMPFA scheme. This explains why the Newton method frequently passes points of nondifferentiability of the NLMPFA weights. In addition, for the checkerboard mesh there are some cells and faces where some of the coefficients ${\alpha}_{K,\sigma {\sigma}^{\mathrm{\prime}}}$ in the subflux definition (20) are negative. These two arguments might be the reason why the NLMPFA scheme needs 623 Newton iterations on the finest refinement level.
In the last two test cases, it has been observed that the NLTPFA, NLMPFA, and AvgMPFA that use the conormal decomposition behave similar to wellestablished linear schemes such as the MPFAO, MFD, or HFV schemes. Especially for highly heterogeneous tensors, the convergence behavior of the NLTPFA scheme is similar to the behavior of the MPFAO or HFV scheme. In addition, the Newton converges much better for the NLTPFA scheme than for the NLMPFA scheme.
Remark 6. For the above test cases, the classical linear TPFA method does not converge. This is wellknown for nonKorthogonal grids, see for example [6]. Therefore, the TPFA scheme has not been considered so far.
3.2 Linearity preservation
Within this section, the linearity preservation property of the schemes is investigated, which is a good indicator for consistency of schemes. This example is based on our recent work [21]. The considered domain and the grid are shown in Figure 4 (right). The domain consists of two subdomains ${\mathrm{\Omega}}_{1}$ and ${\mathrm{\Omega}}_{2}$. The transition from ${\mathrm{\Omega}}_{1}$ to ${\mathrm{\Omega}}_{2}$ is located at $x=0.6$, and the tensors are chosen as$${\mathbf{\Lambda}}_{1}=\left(\begin{array}{lll}3& 1& 0\\ 1& 3& 0\\ 0& 0& 1\end{array}\right),\hspace{1em}{\mathbf{\Lambda}}_{2}=\left(\begin{array}{lll}10& 3& 0\\ 3& 10& 0\\ 0& 0& 1\end{array}\right).$$(51)
Fig. 4 Exact solution for linearity preservation test case as defined in (52) (left); Grid used for the spatial discretization (right). Dirichlet conditions are set at the domain boundary, equal to the exact solution $\stackrel{\u0305}{u}$ (modified after [21]). 
The exact solutions in the subdomains are$${\stackrel{\u0305}{u}}_{1}=14x+y+z,\hspace{1em}{\stackrel{\u0305}{u}}_{2}=4x+y+z+6.$$(52)
Figure 4 (left) depicts the exact solution. Please note that the exact solution and the corresponding flux function are globally continuous within the domain. It can also be seen that the grid is nonmatching at the transition of the subdomains. Such nonmatching grids often occur in faulted geological environments. The grid in Figure 4 is defined by means of the standard cornerpoint grid format and has been generated with the Matlab Reservoir Simulation Toolbox (MRST) [44]. To read in the grid, the opmgrid module from the Open Porous Media (OPM) initiative ^{2} is used, which supports the standard cornerpoint grid format, see [24, 46] for more information about cornerpoint grids.
Table 3 lists the discrete error norms, the number of entries in the Jacobian matrix (noe), and the number of Newton iterations needed for the simulation run. It can be seen that all schemes except the TPFA scheme reproduce the exact solution, because the errors are within the range of the nonlinear and linear solver tolerance, whereas the errors of the linear TPFA scheme are approximately five orders of magnitude higher. It is wellknown that the errors of the linear TPFA scheme are in $\mathcal{O}\left(1\right)$ for nonKorthogonal grids. However, the improved accuracy of the other schemes comes with the cost of a larger face flux stencil, which is the reason why the corresponding Jacobian matrices are denser than the one of the TPFA scheme. When using Picard’s method instead of Newton’s method, the number of entries is the same for the NLTPFA and TPFA scheme. It can also be seen that the MFD and HFV schemes result in a four times denser matrix than the NLTPFA, NLMPFA, or AvgMPFA schemes. This is due to the fact that the calculation of the coefficient matrix (37) for the MFD and HFV schemes, in general, include all face information, i.e. ${\mathcal{S}}_{K,\sigma}={\mathcal{E}}_{K}$. Additionally, flux conservation is weakly enforced using the additional equation (11) that is assembled into the global discretization matrix. The MPFAO or Box scheme cannot be easily applied to this test case because of the nonmatching interface between the subdomains.
Discrete error norms, number of entries in the Jacobian matrix (noe), and the number of Newton iterations (NIt) needed for the linearity preservation test case.
3.3 Discrete extremum principles
The following two examples investigate whether the schemes satisfy discrete extremum principles. As already mentioned before, in general, consistent linear schemes do not satisfy extremum principles, which motivates the construction of nonlinear schemes as introduced above. The setups for these test cases are taken from [21].
First test case
The first example investigates a test case without a source term. The domain and the grid are shown in Figure 5, with an inner and an outer boundary. The Dirichlet values u = 10^{5} and u = 0 are set at the inner and outer boundaries, respectively. Therefore, the solution is expected to be within these bounds. The tensor Λ is chosen homogeneous but anisotropic as$$\mathbf{\Lambda}=\mathbf{R}\left(\frac{\pi}{6}\right)\left(\begin{array}{ll}1000& 0\\ 0& 1\\ & \end{array}\right){\mathbf{R}}^{1}\left(\frac{\pi}{6}\right),\mathrm{with}\mathbf{R}\left(\alpha \right)=\left(\begin{array}{ll}\mathrm{cos}\alpha & \mathrm{sin}\alpha \\ \mathrm{sin}\alpha & \mathrm{cos}\alpha \\ & \end{array}\right).$$(53)
Fig. 5 Unstructured grid used for the first discrete extremum principle test case. 
Figure 6 shows the numerical solutions of the different schemes on a three times refined grid. The TPFA scheme produces no under or overshoots, and therefore satisfies extremum principles. This is wellknown since the resulting discretization matrix is an Mmatrix. However, the linear TPFA scheme is not consistent for this test case which explains why the anisotropy is not reflected by the TPFA solution. The maximum principle is violated by the MPFAO, MFD, and HFV schemes, with overshoots of 8% for the MPFAO scheme and up to 4% for the MFD and HFV schemes. The minimum principle is violated by all consistent linear schemes. The undershoots of the AvgMPFA scheme are above 4%, those of the Box scheme above 2%. The undershoots of the other schemes are below 1%. For this test case, both nonlinear schemes satisfy the discrete extremum principles. The small negative undershoots of the NLMPFA scheme are caused by Newton’s method. These undershoots can be prevented by using other nonlinear solvers such as Picard’s method or enhanced solvers [47]. Furthermore, we also point out that the NLTPFA scheme in general does not satisfy the minimum principle. A similar test case has been considered in [47, 48] with outer Dirichlet boundary conditions above zero. In that case undershoots can also be observed by the NLTPFA scheme.
Fig. 6 Solution of Box, AvgMPFA, NLTPFA, and NLMPFA schemes (first row); solution of MPFAO, TPFA, MFD, and HFV schemes (second row) for the first discrete extremum principle test case calculated on the grid that arises after three times refinement of the grid shown in Figure 5. The tensor is specified in (53). The range of the different solutions is depicted below the different drawings. 
Second test case
The previous example has demonstrated that the NLTPFA and NLMPFA schemes are positivitypreserving. In the following example, it is demonstrated that the NLMPFA scheme satisfies the maximum principle, in contrast to the NLTPFA scheme. This test case has been introduced in [49]. The boundary conditions are Neumann noflow on the whole boundary, i.e. (Λ∇u) · n = 0 on ∂Ω, and Ω = [0, 1]^{2} is discretized with a regular Cartesian grid, with 77 × 77 cells. The tensor Λ is chosen as (53) with α = 67.5°. We add the constraints u = 0 and u = 1 at the cells with centers $(\frac{7}{22},\frac{1}{2})$ and $(\frac{15}{22},\frac{1}{2})$, respectively. Therefore, the solution is bounded, such that $0\le \stackrel{\u0305}{u}\le 1$.
Figure 7 depicts the numerical solutions of the different schemes as surface plots. The TPFA obviously again satisfies the minimum and maximum principle due to the Mmatrix property but again the anisotropy is not correctly reproduced. The only consistent scheme that fulfills the maximum principle is the NLMPFA scheme, all the other schemes produce overshoots. The Box and AvgMPFA schemes produce the highest over and undershoots. Those of the Box scheme are above 26%. The overshoots of the NLTPFA scheme are higher than those of the MPFAO, MFD, and HFV schemes. But the NLTPFA scheme produces no undershoots since it is positivitypreserving. Furthermore, the undershoots of the MPFAO, MFD, and HFV are below 3%.
Fig. 7 Solution of Box, AvgMPFA, NLTPFA, and NLMPFA schemes (first row); solution of MPFAO, TPFA, MFD, and HFV schemes (second row) for the second discrete extremum principle test case. The solutions are depicted as a surface plot. The tensor $\mathbf{\Lambda}$ used for this test case is specified in (53) with α = 67.5°. On the whole boundary Neumann noflow conditions are set, whereby two wells are located within the domain where the values of u are fixed. The range of the different solutions is depicted below the different drawings. 
The above test cases exhibit how nonlinear schemes are capable to reproduce physical solutions, whereas linear schemes can produce negative values. When solving highly complex partial differential equations, where secondary variables nonlinearly depend on primary variables, such negative values can strongly influence the efficiency of the scheme, in terms of linear and nonlinear solver convergence. It should be mentioned again that the NLTPFA is only positivitypreserving, whereas the NLMPFA satisfies discrete extremum principles.
3.4 Benchmark: Northeast German Basin
The following example is a synthetic model inspired by the 3D Northeast German Basin model presented in [50]. This test case has recently been published in [21], where the results of the AvgMPFA, NLTPFA, NLMPFA, TPFA, and Box schemes have been presented. Here, we additionally present the results of the MPFAO, MFD, and HFV schemes. However, the model setup is analogous to [21]. The data and the approximate geometry of the basin are provided by IFPEN. The different facies of the Northeast German Basin are depicted in Figure 8, where the domain is reflected such that −z is oriented in depth direction.
Fig. 8 Different facies of the Northeast German Basin, where the domain is scaled by a factor of five in zdirection and the different layers are shifted horizontally for better visibility of the features. The domain lengths in coordinate directions are approximately 169 km (in the xdirection), 165 km (in the ydirection), and 10.87 km (in the zdirection). The domain has been reflected at the zplane such that −z corresponds to the depth axis. At the top and bottom boundaries, Dirichlet conditions are set to 281.15 K and 423.15 K, respectively, whereas Neumann noflow conditions are set elsewhere (modified after [21]). 
Salt diapirs within this model create highly conductive regions, as shown in Figure 9, leading to thermal anomalies. A robust discretization with respect to the grid is required for this type of structure, in order to evaluate the temperature field and to perform thermohaline simulations.
Fig. 9 Thermal conductivity, calculated with the law (54) (choosing Λ _{ w } = 0.6) and porosity distribution of the Northeast German Basin, whereby the domain is scaled by a factor of five in zdirection. The salt domes correspond to the highly conductive regions (modified after [21]). 
Here, the stationary heat equation is solved, where corresponds to the thermal conductivity [W/(m · K)] and u to the temperature T [K]. The thermal conductivity is computed with the following law [51]$$\mathbf{\Lambda}={\mathbf{\Lambda}}_{w}^{\varphi}{\mathbf{\Lambda}}_{s}^{1\varphi}\mathbf{I},$$(54)where I is the identity matrix, Λ _{ w } and Λ _{s} denote the water and rock conductivities, and $\varphi $ the porosity. The conductivity of water is set to Λ _{ w } = 0.6. The thermal properties of the different facies are listed in Table 4. The porosity distribution and the thermal conductivities are shown in Figure 9. At the top and bottom boundaries, Dirichlet conditions are set to 281.15 K and 423.15 K, respectively, whereas Neumann noflow conditions are set elsewhere. The grid consists of 864,435 cells.
Thermal properties of facies of the Northeast German Basin.
Figure 10 shows the numerical solutions of the MPFAO and TPFA schemes. Additionally, the absolute difference between the MPFAO and the NLTPFA, TPFA, HFV, and Box schemes is depicted in this figure. It is observed that the TPFA scheme differs the most from the other schemes. The largest differences occur at the salt domes, where it seems that the TPFA scheme overestimates the temperature values. The NLTPFA and HFV schemes are in good agreement with the MPFAO results. The largest difference between the NLTPFA and MPFAO scheme is below 7.88 K, and for the HFV scheme below 2.29 K. This can also be seen in Table 5, which lists the relative discrete errors $\frac{\left\right{u}_{1}{{u}_{2}\left\right}_{{L}^{2}}}{{T}_{\mathrm{ref}}\sqrt{\left\mathrm{\Omega}\right}}$ between the schemes, with the total domain volume Ω ≈ 1.75e+14 m^{3} and the reference temperature T _{ref} = 352.15 K. It can be observed that the AvgMPFA, NLMPFA, and NLTPFA schemes produce similar solutions. The same applies to the HFV and MFD schemes. Furthermore, the TPFA and the Box schemes differ the most from the other schemes, which is in agreement with the results shown in Figure 10. Again, the number of entries of the NLTPFA, NLMPFA, and AvgMPFA is approximately twice the number of the TPFA scheme. Moreover, the most dense matrices are those of the HFV and MFD schemes.
Fig. 10 Solution of MPFAO and TPFA schemes (first row). Absolute difference between the MPFAO and the NLTPFA, TPFA, HFV, and Box schemes (second and third row). The results are shown for a part of the domain. 
Discrete relative errors $\frac{\left\right{u}_{1}{{u}_{2}\left\right}_{{L}^{2}}}{{T}_{\mathrm{ref}}\sqrt{\mathrm{\Omega}}}$ between the different schemes. Number of entries (noe) and Newton iterations (NIt).
4 Conclusion
Within this article, we have presented an abstract discretization framework for elliptic equations. Furthermore, it has been demonstrated that various finitevolume schemes fit into this framework, with the difference that hybrid schemes introduce additional face unknowns, whereas cellcentered schemes eliminate these face unknowns by using trace reconstruction operators. Properties like consistency, coercivity, extremum principles, and sparsity have been defined in Section 2.2 for schemes that belong to this discretization framework. Consistency and coercivity are particularly essential to prove the convergence of such schemes. Hybrid schemes are coercive by construction, whereas cellcentered schemes are generally not coercive. In Section 2.3, besides the linear AvgMPFA scheme, a monotone Nonlinear TwoPoint Flux Approximation (NLTPFA) and an extremumprinciplespreserving MultiPoint Flux Approximation (NLMPFA) have been presented. In Section 2.4, we have briefly introduced the wellknown mimetic [13, 52] and hybrid finitevolume [8] schemes. Note that all schemes have been implemented into the opensource simulator DuMu^{x} [40], allowing for a comparison of the different methods within the same software framework. The schemes have been analyzed numerically in Section 3. For the two test cases considered, involving a mildly and highly anisotropic tensor, it has been demonstrated that all schemes, except for the TPFA, show a similar convergence behavior and that all schemes seem to be coercive. The consistency of the NLTPFA, NLMPFA, AvgMPFA, MFD, and HFV schemes, for piecewise linear solutions, has been investigated in Section 3.2 for a nonmatching grid. In addition, we have shown in Section 3.3 that linear schemes are in general neither positivitypreserving nor satisfy discrete minimum or maximum principles, in contrast to nonlinear schemes. In the last example the Northeast German Basin has been considered, showing that the NLTPFA and NLMPFA schemes result in similar solutions than other wellestablished consistent schemes such as the MPFAO, MFD, or HFV schemes.
Acknowledgments
As previously mentioned, some of the presented test cases are based on our recent publication [21] that has been produced in collaboration with Léo Agélas and Guillaume Enchéry from IFPEN. Therefore, the authors would like to thank both of them for their support and the many fruitful conversations. Furthermore, we thank the German Research Foundation (DFG) for funding this work within SFB 1313.
References
 Di Pietro D.A. (2012) Cell centered galerkin methods for diffusive problems, Math. Model. Numer. Anal. 46, 1, 111–144. [CrossRef] [EDP Sciences] [Google Scholar]
 Aavatsmark I., Barkve T., Bøe Ø., Mannseth T. (1996) Discretization on nonorthogonal, quadrilateral grids for inhomogeneous, anisotropic media, J. Comput. Phys. 127, 1, 2–14. [Google Scholar]
 Agélas L., Di Pietro D., Droniou J. (2010) The G method for heterogeneous anisotropic diffusion on general meshes, M2AN Math. Model. Numer. Anal. 44, 4, 597–625. [Google Scholar]
 Agélas L., Guichard C., Masson R. (2010) Convergence of finite volume MPFA O type schemes for heterogeneous anisotropic diffusion problems on general meshes, Int. J. Finite Vol. 7, 2, 1–33. [Google Scholar]
 Agélas L., Masson R. (2008) Convergence of the finite volume MPFA O scheme for heterogeneous anisotropic diffusion problems on general meshes, C. R. Math. 346, 17, 1007–1012. [CrossRef] [Google Scholar]
 Edwards M., Rogers C. (1998) Finite volume discretization with imposed flux continuity for the general tensor pressure equation, Comput. Geosci. 2, 4, 259–290. [Google Scholar]
 Wolff M., Cao Y., Flemisch B., Helmig R., Wohlmuth B. (2013) Multipoint flux approximation Lmethod in 3d: numerical convergence and application to twophase flow through porous media, Radon Ser. Comput. Appl. Math., De Gruyter 12, 39–80. [Google Scholar]
 Eymard R., Gallouët T., Herbin R. (2010) Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes. SUSHI: a scheme using stabilization and hybrid interfaces, IAM J. Num. Anal. 30, 4, 1009–1043. [Google Scholar]
 Eymard R., Herbin R., Guichard C., Masson R. (2012) Vertexcentred discretization of multiphase compositional Darcy flows on general meshes, Comput. Geosci. 16, 4, 987–1005. [Google Scholar]
 Arnold D., Brezzi F. (1985) Mixed and nonconforming finite element methods: implementation, postprocessing and error estimates, Math. Model. Numer. Anal. 19, 1, 7–32. [Google Scholar]
 Raviart P.A., Thomas J.M. (1977) A mixed finite element method for 2nd order elliptic problems, in: I. Galligani, E. Magenes (eds), Mathematical aspects of finite element methods, Springer, Berlin, Heidelberg, pp. 292–315. [Google Scholar]
 Brezzi F., Lipnikov K., Shashkov M. (2005) Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes, SIAM J. Numer. Anal. 43, 5, 1872–1896. [Google Scholar]
 Brezzi F., Lipnikov K., Simoncini V. (2005) A family of mimetic finite difference methods on polygonal and polyhedral meshes, Math. Models Methods Appl. Sci. 15, 10, 1533–1551. [Google Scholar]
 Droniou J., Eymard R., Gallouët T., Herbin R. (2010) A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods, Math. Models Methods Appl. Sci. 20, 2, 265–295. [Google Scholar]
 Droniou J., Eymard R., Herbin R. (2016) Gradient schemes: generic tools for the numerical analysis of diffusion equations, ESAIM Math. Model. Numer. Anal. 50, 3, 749–781. [Google Scholar]
 Vohralík M., Wohlmuth B.I. (2013) Mixed finite element methods: implementation with one unknown per element, local flux expressions, positivity, polygonal meshes, and relations to other methods, Math. Models Methods Appl. Sci. 23, 5, 803–838. [Google Scholar]
 Danilov A., Vassilevski Y. (2009) A monotone nonlinear finite volume method for diffusion equations on conformal polyhedral meshes, Russ. J. Numer. Anal. Math. Modelling 24, 3, 207–227. [CrossRef] [Google Scholar]
 Lipnikov K., Svyatskiy D., Vassilevski Y. (2009) Interpolationfree monotone finite volume method for diffusion equations on polygonal meshes, J. Comput. Phys. 228, 3, 703–716. [Google Scholar]
 Lipnikov K., Svyatskiy D., Vassilevski Y. (2012) Minimal stencil finite volume scheme with the discrete maximum principle, Russ. J. Numer. Anal. Math. Modelling 27, 4, 369–385. [CrossRef] [Google Scholar]
 Potier C.L. (2005) Schéma volumes finis monotone pour des opérateurs de diffusion fortement anisotropes sur des maillages de triangles non structurés, C.R. Math. 341, 12, 787–792. [CrossRef] [MathSciNet] [Google Scholar]
 Schneider M., Agélas L., Enchéry G., Flemisch B. (2017) Convergence of nonlinear finite volume schemes for heterogeneous anisotropic diffusion on general meshes, J. Comput. Phys. 351, Supplement C, 80–107. [Google Scholar]
 Yuan G., Sheng Z. (2008) Monotone finite volume schemes for diffusion equations on polygonal meshes, J. Comput. Phys. 227, 12, 6288–6312. [Google Scholar]
 Schneider M., Flemisch B., Helmig R. (2017) Monotone nonlinear finitevolume method for nonisothermal twophase twocomponent flow in porous media, Int. J. Numer. Methods Fluids 84, 6, 352–381. [Google Scholar]
 Schneider M., Flemisch B., Helmig R., Terekhov K., Tchelepi H. (Apr 2018) Monotone nonlinear finitevolume method for challenging grids, Comput. Geosci. 22, 2, 565–586. [Google Scholar]
 Schneider M., Gläser D., Flemisch B., Helmig R. (2017) Nonlinear finitevolume scheme for complex flow processes on cornerpoint grids, in: C. Cancès, P. Omnes (eds), Finite volumes for complex applications VIII – Hyperbolic elliptic and parabolic problems, Springer International Publishing, pp. 417–425. [CrossRef] [Google Scholar]
 Droniou J., Eymard R., Gallouët T., Guichard C., Herbin R. (2018) The gradient discretisation method, HAL, https://hal.archivesouvertes.fr/hal01382358. [Google Scholar]
 Eymard R., Gallouët T., Herbin R. (2000) Finite volume methods, Handbook Numer. Anal. 7, 713–1018. [Google Scholar]
 Evans L. (1998) Partial differential equations, Graduate Studies in Mathematics. American Mathematical Society. [Google Scholar]
 Kuzmin D. (2010) A guide to numerical methods for transport equations, Universität Nürnberg, http://www.mathematik.unidortmund.de/~kuzmin/Transport.pdf. [Google Scholar]
 Berman A., Plemmons R.J. (1994) Nonnegative matrices in the mathematical sciences, Vol. 9, SIAM, Philadelphia, PA. [CrossRef] [Google Scholar]
 Droniou J., Potier C.L. (2011) Construction and convergence study of schemes preserving the elliptic local maximum principle, SIAM J. Numer. Anal. 49, 2, 459–490. [Google Scholar]
 Potier C.L. (2009) A nonlinear finite volume scheme satisfying maximum and minimum principles for diffusion operators, Int. J. Finite Vol. 6, 2, 1–20. [Google Scholar]
 Lipnikov K., Manzini G., Shashkov M. (2014) Mimetic finite difference method, J. Comput. Phys. 257, 1163–1227. [Google Scholar]
 Di Pietro D.A., Vohralík M. (2014) A review of recent advances in discretization methods, a posteriori error analysis, and adaptive algorithms for numerical modeling in geosciences, Oil Gas Sci. Technol.  Rev. IFP Energies nouvelles 69, 4, 701–729. [CrossRef] [Google Scholar]
 Droniou J. (2014) Finite volume schemes for diffusion equations: introduction to and review of modern methods, Math. Models Methods Appl. Sci. 24, 8, 1575–1619. [Google Scholar]
 Agélas L., Eymard R., Herbin R. (2009) A ninepoint finite volume scheme for the simulation of diffusion in heterogeneous media, C. R. Math. 347, 11, 673–676. [CrossRef] [MathSciNet] [Google Scholar]
 Aavatsmark I. (2002) An introduction to multipoint flux approximations for quadrilateral grids, Comput. Geosci. 6, 3–4, 405–432. [Google Scholar]
 Hackbusch W. (1989) On first and second order box schemes, Computing 41, 4, 277–296. [CrossRef] [MathSciNet] [Google Scholar]
 Helmig R. (1997) Multiphase flow and transport processes in the subsurface: a contribution to the modeling of hydrosystems, SpringerVerlag, Berlin. [Google Scholar]
 Hommel J., Ackermann S., Beck M., Becker B., Class H., Fetzer T., Flemisch B., Gläser D., Grüninger C., Heck K., Kissinger A., Koch T., Schneider M., Seitz G., Weishaupt K. (2016) DuMuX 2.10.0, https://doi.org/10.5281/zenodo.159007. [Google Scholar]
 Blatt M., Burchardt A., Dedner A., Engwer C., Fahlke J., Flemisch B., Gersbacher C., Gräser C., Gruber F., Grüninger C., Kempf D., Klöfkorn R., Malkmus T., Müthing S., Nolte M., Piatkowski M., Sander O. (2016) The distributed and unified numerics environment, version 2.4, Arc. Numer. Softw. 4, 100, 13–29. [Google Scholar]
 Brezzi F., Lipnikov K., Shashkov M. (2006) Convergence of mimetic finite difference method for diffusion problems on polyhedral meshes with curved faces, Math. Models Methods Appl. Sci. 16, 2, 275–297. [Google Scholar]
 Fort J., Fürst J., Halama J., Herbin R., Hubert F. (2011) Finite Volumes for Complex Applications VI: Problems and Perspectives, Springer, Heidelberg. [CrossRef] [Google Scholar]
 Krogstad S., Lie K., Møyner O., Nilsen H.M., Raynaud X., Skaflestad B. (2015) MRSTAD – an opensource framework for rapid prototyping and evaluation of reservoir simulation problems. SPE Reservoir Simulation Symposium, Houston, TX, February 23–25, Society of Petroleum Engineers. [Google Scholar]
 Alkämper M., Dedner A., Klöfkorn R., Nolte M. (2016) The DUNEALUGrid module, Arc. Numer. Softw. 4, 1, 1–28. [Google Scholar]
 Aarnes J.E., Krogstad S., Lie K.A. (2008) Multiscale mixed/mimetic methods on cornerpoint grids, Comput. Geosci. 12, 3, 297–315. [Google Scholar]
 Terekhov K., Mallison B., Tchelepi H. (2017) Cellcentered nonlinear finitevolume methods for the heterogeneous anisotropic diffusion problem, J. Comput. Phys. 330, 245–267. [Google Scholar]
 Kapyrin I., Nikitin K., Terekhov K., Vassilevski Y. (2014) Nonlinear monotone FV schemes for radionuclide geomigration and multiphase flow models, in: Finite volumes for complex applications VIIElliptic, parabolic and hyperbolic problems, Springer, pp. 655–663. [Google Scholar]
 Aavatsmark I., Eigestad G., Mallison B., Nordbotten J. (2008) A compact multipoint flux approximation method with improved robustness, Numer. Methods Partial Differ. Equ. 24, 5, 1329–1360. [Google Scholar]
 Scheck M., Bayer U. (1999) Evolution of the Northeast German Basin – inferences from a 3D structural model and subsidence analysis, Tectonophysics 3, 3, 145–169. [Google Scholar]
 Woodside W., Messmer J. (1961) Thermal conductivity of porous media. I. Unconsolidated sands, J. Appl. Phys. 32, 9, 1688–1699. [Google Scholar]
 Lipnikov K., Manzini G., Svyatskiy D. (2011) Analysis of the monotonicity conditions in the mimetic finite difference method for elliptic problems, J. Comput. Phys. 230, 7, 2620–2642. [Google Scholar]
All Tables
Discrete error norms, number of entries in the Jacobian matrix (noe), and the number of Newton iterations (NIt) needed for the linearity preservation test case.
Discrete relative errors $\frac{\left\right{u}_{1}{{u}_{2}\left\right}_{{L}^{2}}}{{T}_{\mathrm{ref}}\sqrt{\mathrm{\Omega}}}$ between the different schemes. Number of entries (noe) and Newton iterations (NIt).
All Figures
Fig. 1 Hexahedral meshes used for the convergence test cases (modified after [24]). 

In the text 
Fig. 2 Logarithm of L ^{2}error (left) and discrete H ^{1}error (right) for convergence test case one. From top to bottom: checkerboard, random, nonconvex, curved mesh. 

In the text 
Fig. 3 Logarithm of L ^{2}error (left) and discrete H ^{1}error (right) for convergence test case one. From top to bottom: checkerboard, random, nonconvex, curved mesh. 

In the text 
Fig. 4 Exact solution for linearity preservation test case as defined in (52) (left); Grid used for the spatial discretization (right). Dirichlet conditions are set at the domain boundary, equal to the exact solution $\stackrel{\u0305}{u}$ (modified after [21]). 

In the text 
Fig. 5 Unstructured grid used for the first discrete extremum principle test case. 

In the text 
Fig. 6 Solution of Box, AvgMPFA, NLTPFA, and NLMPFA schemes (first row); solution of MPFAO, TPFA, MFD, and HFV schemes (second row) for the first discrete extremum principle test case calculated on the grid that arises after three times refinement of the grid shown in Figure 5. The tensor is specified in (53). The range of the different solutions is depicted below the different drawings. 

In the text 
Fig. 7 Solution of Box, AvgMPFA, NLTPFA, and NLMPFA schemes (first row); solution of MPFAO, TPFA, MFD, and HFV schemes (second row) for the second discrete extremum principle test case. The solutions are depicted as a surface plot. The tensor $\mathbf{\Lambda}$ used for this test case is specified in (53) with α = 67.5°. On the whole boundary Neumann noflow conditions are set, whereby two wells are located within the domain where the values of u are fixed. The range of the different solutions is depicted below the different drawings. 

In the text 
Fig. 8 Different facies of the Northeast German Basin, where the domain is scaled by a factor of five in zdirection and the different layers are shifted horizontally for better visibility of the features. The domain lengths in coordinate directions are approximately 169 km (in the xdirection), 165 km (in the ydirection), and 10.87 km (in the zdirection). The domain has been reflected at the zplane such that −z corresponds to the depth axis. At the top and bottom boundaries, Dirichlet conditions are set to 281.15 K and 423.15 K, respectively, whereas Neumann noflow conditions are set elsewhere (modified after [21]). 

In the text 
Fig. 9 Thermal conductivity, calculated with the law (54) (choosing Λ _{ w } = 0.6) and porosity distribution of the Northeast German Basin, whereby the domain is scaled by a factor of five in zdirection. The salt domes correspond to the highly conductive regions (modified after [21]). 

In the text 
Fig. 10 Solution of MPFAO and TPFA schemes (first row). Absolute difference between the MPFAO and the NLTPFA, TPFA, HFV, and Box schemes (second and third row). The results are shown for a part of the domain. 

In the text 