© M. Schneider et al., published by IFP Energies nouvelles, 2018

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 application-dependent 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 cell-centered finite-volume methods, such as cell-centered Galerkin methods [1] or Multi-Point Flux Approximation (MPFA) methods [2–7]; or hybrid, mixed, and mimetic (HMM) schemes, such as the Hybrid Finite-Volume (HFV) schemes [8, 9], the Mixed Finite-Element (MFE) [10, 11] or the Mimetic Finite-Difference (MFD) methods [12, 13]. All these schemes are closely related and they can be either presented within a finite-volume framework (which is done in this article) or within a finite-element 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. corner-point grids. However, this comes with the cost of additional unknowns. Recently, monotone or discrete extremum-principles-preserving schemes have been developed in [17–22], and applied to highly complex porous media applications in [23–25].

In this work, we compare different finite-volume 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 cell-centered 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 three-dimensional benchmark case. Finally, in Section 4 we summarize the results obtained in this article.

2 Discretization schemes

Let , be an open bounded connected polygonal domain with boundary ∂Ω, and d-dimensional measure |Ω|. In the following, we consider the elliptic problem(1)

Here, we assume that f ∈ L ²(Ω) and Λ is a symmetric tensor-valued function such that (s.t.) the spectrum of Λ(x) is contained in [α ₀, β ₀], with 0 < α ₀ < β ₀ < +∞, for almost every (a.e.) . 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 is a triplet , where

1. are the cells s.t. . For each cell ,  > 0 denotes the cell volume and its boundary. defines the discretization length.

2. are the faces such that each face σ is included in a hyperplane of with (d − 1)-dimensional measure |σ| > 0. For each cell is the subset of such that . contains the cells that share the face σ; the sets of inner and boundary faces are denoted by and , respectively.

3. are the cell centers (not required to be the barycenters) s.t. x _K  ∈ K and K is star-shaped with respect to x _K . For all , is given as the Euclidean distance between x _K and σ.

For a more detailed definition of an admissible discretization see [3, 4, 21]. Let be any set, then the cardinality of this set is indicated with . With this definition, the number of cells, i.e. the cardinality of the set , is given by . The face evaluation points are denoted as x _σ , (not required to be the barycenters). With we define the averaged tensor on cell K, where the integral is meant component-wise. In addition, it is assumed that for all . 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 yields(2)

Using Definition 1, the integral on the left of equation (2) can be spitted into integrals over the faces σ, which results in(3)

To discretize equation (3) we need to define flux approximations F _K,σ , for each cell K and face , such that(4)

i.e. F _K,σ is an approximation of the exact flux.

For this purpose, we define the following discrete solution spaces(5) (6) (7)

For any element , the cell unknowns are denoted by and the face unknowns are denoted by . The space that takes into account the homogeneous zero Dirichlet boundary conditions is accordingly given as(8)

Additionally, the space of piecewise constant functions on is defined as

Therefore, for all and for all will denote the constant value of v on K. Furthermore, the extraction operator is defined for all as(9)

The presented finite-volume schemes differ in the choice of the fluxes and the choice of the solution space . Therefore, for all , and for all , let be a numerical flux function that approximates the flux flowing out of K through σ such that the finite-volume scheme reads: Find u ∈ X _h s.t.(10) (11)

Cell-centered schemes:

For cell-centered schemes the solution space is given as(12)whereby, for each face σ, is a trace reconstruction operator. Here, is the space of linear operators that map some element to a constant value associated with the corresponding face. This means that for cell-centered schemes, the face unknowns are eliminated using some trace reconstruction operators .

Here, only locally mass-conservative cell-centered schemes are considered such that for all , with (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(14)

This means that HMM schemes introduce additional face unknowns, in contrast to cell-centered schemes where these unknowns are eliminated.

For the locally mass-conservative cell-centered 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 s.t.(15)with the form, for all (u, v) ∈ [X _h ]²,(16)

Therefore, each solution of problem (10) and (11) also solves problem (15). Starting from the discrete problem (15), with defined by (16), equation (11) is fulfilled for the cell-centered schemes by construction, whereas for the HMM schemes it is obtained by taking for all σ′ ≠ σ, and v _K  = 0 for all K

Equation (10) is obtained by taking, for each and for all 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 . 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 cell-centered schemes (with discrete solution space (12)) and for HMM schemes (with discrete solution space (14)).

Consistency

There is no unique discretization-independent definition of consistency. For example, finite-volume schemes are in general not consistent in the finite-difference context [27]. In general, a scheme is consistent if the truncation error between discrete and continuous operators goes to zero for . In the context of finite-volume schemes we say that a scheme is consistent if the numerical flux approximates the exact flux for regular functions, meaning that(17)where is a test function space which is assumed to be dense in , and are the discrete and exact flux functions, respectively. Furthermore, is defined as for all , and for all . More details can for example be found in [21].

Coercivity

A scheme is called uniformly coercive if the form , defined in (16), satisfies the following estimate:(18)with some appropriate discrete norm , and γ > 0 that is independent of u and . 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 is monotone, which means that all entries of its inverse are non-negative i.e. . The monotonicity property is for example satisfied by M-matrices. M-matrices are monotone matrices with non-positive off-diagonal entries (see, e.g. [30]).

Sparsity

For solving large-scale 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 . The stencil denotes the noe of each row. Therefore, schemes resulting in small stencils are important for large-scale problems.

2.3 Cell-centered schemes

Within this section, we present different cell-centered finite-volume schemes. Here, the discrete solution space is given as (12), and the face values are eliminated using some trace reconstruction operators . As already mentioned, only locally mass-conservative schemes are considered here such that equation (11) holds by construction.

Family of weighted schemes

An established idea to obtain monotone or extremum-principles-preserving schemes, as those developed in [17–20, 22, 23, 31, 32], is to combine, for each interior edge , with , two consistent linear flux approximations and as a convex combination, using solution-dependent weights. Thus, the final flux F _K,σ (u) is given as(19)

For any and , the linear flux is built in order to ensure a strong consistency property, such that the total flux satisfies a weak form of consistency, see [21] for detailed information. A general definition of these sub-fluxes is(20)with face stencil and trace reconstruction operator .

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 is also nonlinear in its first argument. This can be prevented by introducing an additional argument for the fluxes and correspondingly for the form . The details are not presented here, for a detailed description see [21].

The final flux function (19) is locally mass-conservative, whereas the sub-fluxes are in general not mass-conservative. 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 finite-volume scheme that is in the following denoted as AvgMPFA scheme.

NLTPFA schemes

To derive a Nonlinear Two-Point Flux Approximation (NLTPFA), the different terms are reordered such that the flux is written as(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(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⁻⁸ but any other small number could also be chosen.

NLMPFA scheme

A Nonlinear Multi-Point Flux Approximation (NLMPFA) is derived by splitting the sub-fluxes into a two-point part and a residual flux part(23)with some appropriate and . Here, the weights are chosen as(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 in the sub-fluxes (20) may also vary). The weights of the NLTPFA scheme are defined by and 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 and which can be written as a sum of solution value differences, i.e. .

Linear TPFA/MPFA schemes

Here, we shortly demonstrate that well-established schemes such as the Two-Point Flux Approximation (TPFA) or the Multi-Point Flux Approximation (MPFA) also fit into this family of schemes. The difference here is that the sub-fluxes are constructed to satisfy flux continuity, i.e. . Therefore, the final flux is given as (20), independent of the weights μ _K,σ , μ _L,σ .

For the TPFA scheme it holds that and the coefficients are(25)

The trace reconstruction operator is(26)such that the final flux reads(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 , as(28)with the face unknowns , the cell unknowns u _K , and coefficients . These coefficients are chosen such that the resulting scheme is coercive and consistent.

Defining the matrices,(29)we obtain the cell flux vector as(30)with , , and e is the vector with entries equal to one. Additionally, let be the conormal matrix, and let be the matrix that contains the scaled distance vectors:(31)

Then, the consistency condition is given by(32)

From this consistency condition, we can derive a general form of the coefficient matrix :(33)with a stabilization matrix such that . Choosing x _σ as the barycenter of face σ and using the geometrical identity(34)it can be seen that(35)

Therefore, the matrix is symmetric if the stabilization matrix is chosen symmetric.

MFD scheme

A simple choice of is(36)with . This results in the final matrix(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 Finite-Volume (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 , 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 such that

The idea of symmetric gradient discretization schemes is to replace the operator ∇ with a consistent approximation such that the form of the discrete problem (15) is given by(38)

Therefore, the scheme is defined by the discrete gradient reconstruction operator . First, let us define a discrete gradient on each cell as(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(40)

This results in the final discrete gradient: For all (41)where 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 (30) can also be identified. Therefore, the hybrid finite-volume 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 finite-volume schemes we refer to [34, 35].

3 Numerical experiments

In this chapter, the behavior of the different finite-volume schemes, that have been presented in the last section, is investigated. The AvgMPFA, NLTPFA, and NLMPFA schemes need the coefficients , that occur in the sub-fluxes (20). These coefficients are calculated using optimization strategies together with a Primal-Dual Simplex Method provided by the open-source library GNU Linear Programming Kit ¹ (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 sub-fluxes (20).

The mimetic finite-difference scheme with fluxes (30) and local cell matrix (37) is denoted as MFD scheme. The hybrid finite-volume scheme defined by the discrete gradients (41) is denoted as HFV scheme. Finally, with TPFA we denote the scheme with fluxes (27) and with MPFA-O the scheme that is described in [37]. The Box method [38, 39] is a vertex-centered finite-volume scheme that uses finite-element basis functions on each cell to calculate fluxes over sub-control volume faces.

All simulations are performed using the open-source 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 finite-volume schemes, where the iteration loop is stopped if the absolute residual is below 10⁻⁵ 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 such that the weak solution has to satisfy (where denotes the trace operator).

A discrete seminorm on the space is given by(42)which is a norm on the space . For the cell-centered schemes, we define the following discrete norm (that takes into account the Dirichlet data)(43)whereas 〈g〉 _σ denotes the average of g on the face σ. Using the Cauchy-Schwarz inequality, one can show that (see [8] for more details)(44)with(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(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 , respectively.

Remark 5. In this section, hexahedral grids are considered (or quadrilateral grids for the case of a two-dimensional domain). Such grids are, in general, not admissible in the sense of Definition 1, because the faces are usually non-planar. 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 three-dimensional 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 well-established 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 MPFA-O, 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 non-convex grid and the curved grid are generated with the Matlab Reservoir Simulation Toolbox (MRST) [44]. For all grids, we use the dune-alugrid module [45]. Except for the checkerboard mesh, all grids exhibit non-planar faces. Therefore, we calculate the integrated normal vectors, see [24]. For the MPFA-O scheme this is applied to the sub-faces 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(47)is prescribed on the unit cube Ω = [0, 1]³. Dirichlet conditions are set on the boundary as on ∂Ω. The permeability tensor is(48)

The discrete L ²-errors and H ¹-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 MPFA-O scheme cannot be applied to hanging nodes on three-dimensional grids, which is why there are no results of the MPFA-O scheme for the checkerboard mesh.

It can be seen that for increasing mesh refinement, the schemes converge with second order accuracy in the L ²-norm and at least first order accuracy in the H ¹-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 MPFA-O and HFV is similar, especially for the random and non-convex grids. In this example, the NLTPFA scheme produces smaller errors than the NLMPFA scheme. Considering Table 1, we observe that the coercivity estimates 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(49)with β = 10⁻² and the exact solution(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 MPFA-O scheme on the checkerboard mesh due to the non-conformity of the grid.

In Figure 3, we can see that the NLTPFA scheme produces now similar errors than the HFV and MPFA-O scheme, and results in the smallest H ¹-errors on the random grid. Considering Table 2, we observe that the coercivity estimates 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 , which are used to define the weights of the NLMPFA scheme. This explains why the Newton method frequently passes points of non-differentiability of the NLMPFA weights. In addition, for the checkerboard mesh there are some cells and faces where some of the coefficients in the sub-flux 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 well-established linear schemes such as the MPFA-O, MFD, or HFV schemes. Especially for highly heterogeneous tensors, the convergence behavior of the NLTPFA scheme is similar to the behavior of the MPFA-O 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 well-known for non-K-orthogonal 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 sub-domains and . The transition from to is located at , and the tensors are chosen as(51)

The exact solutions in the sub-domains are(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 non-matching at the transition of the sub-domains. Such non-matching grids often occur in faulted geological environments. The grid in Figure 4 is defined by means of the standard corner-point grid format and has been generated with the Matlab Reservoir Simulation Toolbox (MRST) [44]. To read in the grid, the opm-grid module from the Open Porous Media (OPM) initiative ² is used, which supports the standard corner-point grid format, see [24, 46] for more information about corner-point 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 well-known that the errors of the linear TPFA scheme are in for non-K-orthogonal 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. . Additionally, flux conservation is weakly enforced using the additional equation (11) that is assembled into the global discretization matrix. The MPFA-O or Box scheme cannot be easily applied to this test case because of the non-matching interface between the sub-domains.

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⁵ 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(53)

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 well-known since the resulting discretization matrix is an M-matrix. 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 MPFA-O, MFD, and HFV schemes, with overshoots of 8% for the MPFA-O 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.

Second test case

The previous example has demonstrated that the NLTPFA and NLMPFA schemes are positivity-preserving. 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 no-flow on the whole boundary, i.e. (Λ∇u) · n = 0 on ∂Ω, and Ω = [0, 1]² 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 and , respectively. Therefore, the solution is bounded, such that .

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 M-matrix 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 MPFA-O, MFD, and HFV schemes. But the NLTPFA scheme produces no undershoots since it is positivity-preserving. Furthermore, the undershoots of the MPFA-O, MFD, and HFV are below 3%.

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 non-linearly 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 positivity-preserving, 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 MPFA-O, 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.

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.

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](54)where I is the identity matrix, Λ _w and Λ _s denote the water and rock conductivities, and 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 no-flow conditions are set elsewhere. The grid consists of 864,435 cells.

Figure 10 shows the numerical solutions of the MPFA-O and TPFA schemes. Additionally, the absolute difference between the MPFA-O 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 MPFA-O results. The largest difference between the NLTPFA and MPFA-O 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 between the schemes, with the total domain volume |Ω| ≈ 1.75e+14 m³ 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.

4 Conclusion

Within this article, we have presented an abstract discretization framework for elliptic equations. Furthermore, it has been demonstrated that various finite-volume schemes fit into this framework, with the difference that hybrid schemes introduce additional face unknowns, whereas cell-centered 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 cell-centered schemes are generally not coercive. In Section 2.3, besides the linear AvgMPFA scheme, a monotone Nonlinear Two-Point Flux Approximation (NLTPFA) and an extremum-principles-preserving Multi-Point Flux Approximation (NLMPFA) have been presented. In Section 2.4, we have briefly introduced the well-known mimetic [13, 52] and hybrid finite-volume [8] schemes. Note that all schemes have been implemented into the open-source 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 non-matching grid. In addition, we have shown in Section 3.3 that linear schemes are in general neither positivity-preserving 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 well-established consistent schemes such as the MPFA-O, 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

All Tables

All Figures