Abstract

A Discrete Duality Finite Volume (DDFV) method to solve on unstructured meshes the flow problems in anisotropic nonhomogeneous porous media with full Neumann boundary conditions is proposed in the present work. We start with the derivation of the discrete problem. A result of existence and uniqueness of a solution for that problem is given thanks to the properties of its associated matrix combined with adequate assumptions on data. Their theoretical properties, namely, stability and error estimates (in discrete energy norms and -norm), are investigated. Numerical test is provided.

1. Introduction and the Model Problem

Efficient schemes are required for addressing flow problems in geologically complex media. The most important criteria of efficiency are mass conservation in grid blocks, accurate approximation of Darcy velocity, capability for dealing with anisotropic flow on unstructured grids and diverse heterogeneities (relevant to absolute permeability, porosity, etc.), and robust and easy implementation. Schemes well known in the literature for meeting many of the previous criteria are the following: Mixed Finite Element methods (see, e.g., [1, 2]), Control-Volume Finite Element methods (see, e.g., [3, 4]), Mimetic Finite Difference methods (see, e.g., [5, 6] and references therein), Cell-Centered Finite Volume methods (see, e.g., [711] and certain references therein; see also [1214]), Multipoint Flux Approximation (see, e.g., [1518] and some contributions to convergence analysis of MPFA O-scheme like [19]) and Discrete Duality Finite Volume methods (DDFV methods for short).

The DDFV methods come in two formulations. The first formulation is based on interface flux computations for primary and dual meshes, accounting with the interface flux continuity (see, e.g., [20, 21]) and the second formulation of DDFV is based on pressure gradient reconstructions over a diamond grid (see [2224]). Note that this second formulation attracted the attention of some mathematicians as Andreianov, Boyer, and Hubert who have greatly contributed to its mathematical development. Indeed, key ideas involved in the pressure gradient reconstruction approach have been generalized by these authors to nonlinear operators of Leray-Lions type. Motivated by the possibility of increasing the order of convergence of the pressure gradient reconstruction method for nonlinear operators, Boyer and Hubert have proposed in [25] the so-called modified DDFV.

Beyond flow problems, we find the applications of DDFV methods in many areas: the numerical modeling of the surface erosion occurring at a fluid/soil interface undergoing a flow process in [26], the discretisation of partial differential equation appearing in image processing in [27], the assessment of nuclear waste repository safety in the context of simulating flow, transport in porous media in [28], and so on.

For presenting our analysis of DDFV method, let us consider the 2D diffusion problem consisting of finding a function which satisfies the following partial differential equation associated with nonhomogeneous Neumann boundary conditions:where is a given open polygonal domain, is its boundary, and are two given functions (defined, resp., in and , at least almost everywhere in Lebesgue measure sense), and is the unit normal vector to outward to . The permeability , with , may be a full matrix depending on solely spatial variables and obeying the following conditions.(i)The primary mesh is such that the discontinuity of lies on mesh interfaces.(ii)Symmetry is(iii)Uniform ellipticity and boundedness arewhere denotes the Euclidian norm in and where are the components of satisfyingwhere denotes the subspace of made of piecewise constant functions defined in .

One also assume that and satisfy the following conditions.(i)Compatibility condition is(ii)Null average condition is(iii)Let us emphasize that the novelty of this work is that Neumann boundary conditions are imposed on the whole boundary, which causes additional difficulties in the analysis.

2. A Finite Volume Formulation of the Model Problem

We briefly present finite volume formulation of the model problem (1)-(2) on unstructured meshes into the same spirit as [20, 25, 2931]. We assume thatthe diffusion matrix is a piecewise constant function over .This assumption is currently used at least in industrial problems (e.g., reservoir simulation problems). Indeed, a subsurface area is made up of a collection of various geologic formations that may be characterized at intermediate scales by averaged full permeability tensors over grid blocks of the primary grid: for more details on this topic, see [32, 33].

2.1. Formulation of the Discrete Problem

First of all, notice that the method under consideration is analyzed in this work for general polygonal domains covered with unstructured matching primary meshes made up of arbitrary convex polygons (see Figure 1). Let us consider some definitions needed in what follows.

Definition 1. A mesh defined on is compatible with the discontinuities of the permeability tensor if these discontinuities are located along the mesh interfaces.

Definition 2. One defines an edge-point as any point located over an edge and different from the extremities of that edge.

Definition 3. Two edge-points and are named “neighboring edge-points” if they share the same vertex in the sense that and belong to two different edges that intersect in .

Let us recall our main objectives in this work:(i)Compute at the cell-points (to be defined later) and at the interior vertices from the mesh the values of the unknown function as a solution (expected unique) of a linear system.(ii)Analyze the stability and the convergence of this solution in some discrete energy norm similar to the one introduced in [31].

In the context of unstructured primary meshes (including square primary meshes with cell-points chosen different from cell-centers), the definition of a discrete energy norm requires that the cell-points lie inside certain perimeters. In this connection, the main steps for defining the cell-points are as follows:(a)Choose arbitrarily a unique point (different from a vertex) on each edge of the mesh . This process generates a finite family of edge-points.(b)Join every pair of neighboring edge-points by a dotted straight line.By this way, we generate an auxiliary mesh denoted .(c)Fix arbitrarily a unique point inside each intersection of a primary cell and an auxiliary cell. These points define a finite family of cell-points.Figure 2 illustrates the location of the edge-points and cell-points within both primary and auxiliary meshes.

Remark 4. Note that, in the 3D framework, for a given primary mesh the corresponding auxiliary mesh is generated very easily as follows. Any primary cell involves a finite number of faces and to each face is assigned one and only one face point lying necessarily on the boundary of . Therefore, one could associate with these face points the smallest convex polygon containing all of them. This is the process leading to generation of the auxiliary mesh associated with the primary mesh under consideration.

We should emphasize that we need to take cell-points inside auxiliary cells (see Figure 2) in view of achieving the following goal which is to define a vector space for the piecewise constant solution (named “weak approximate” solution in the sequel) and to equip it with a discrete energy norm.

Remark 5. Note that there exists a DDFV approach based on pressure gradient reconstructions that can address nonlinear elliptic problems: see [25, 34] for more details. Nevertheless in these works, the choice of edge-points depends strongly on that of cell-points as each edge-point coincides (by Definition 2) with the intersection of a primary edge and the straight line joining two cell-points located in both sides of that edge.

From the boundary-value problem theory (see, e.g., [35]), system (1)-(2) possesses a unique solution in under assumptions (3)–(7). We have assumed that the diffusion coefficient is a piecewise constant function over . The discontinuities of naturally divide into a finite number of convex subdomains . We now make the additional assumption that the restriction over of the exact solution to system (1)-(2), denoted by , satisfies

Let us now focus on a finite volume formulation of problem (1)-(2) in terms of a linear system which should be derived from the elimination of auxiliary unknowns, namely, edge-point pressures, in flux balance equations over primary cells and also dual cells (to be introduced later). This linear system should involve the real numbers and as discrete unknowns which are expected to be reasonable approximations of (cell-point pressures) and (interior vertex pressures), respectively, where and and where stands for dual mesh. We now give a description of the procedure leading to the linear discrete system.

Let be a primary cell, where is the corresponding cell-point. We integrate the two sides of the mass balance equation (1) in . Applying Ostrogradsky’s theorem to the integral in the left-hand side leads to computing the flux across the boundary of . Thanks to a suitable quadrature formula this computation yields a relation involving edge-point pressures. Due to the flux continuity over the mesh interfaces, the edge-point pressures can be eliminated from the previous relation.

As an illustration of this technique for computing the fluxes across primary cell boundaries, we consider the edge associated with the primary cell (see Figure 3).

Let be the diffusion tensor of the primary cell Denoting as the unit normal vector to exterior to , the flux expression over the edge viewed as part of the boundary of is given bywhere is the truncation error and where we introduce the following necessary ingredients:where (resp., ) is the angle defined by the vectors and (resp., and ) and where (resp., ) denotes the unit normal vector to (resp., to ) exterior to the half-plane from containing the point and bordered by the straight line (resp., ). Notice that , . Moreover if the primary mesh () is regular in the sense of Definition 6, we have and therefore , where is a certain real number not depending on .

In what follows, we will need the following notations. We denote by the set of all edge-points (from the primary mesh of course), the subset of made up of internal edge-points, that is, edge-points lying on , the subset of made up of edge-points lying on edges included in the boundary of , and (for ) the subset of made up of edge-points lying on the boundary of the primary cell .

We now introduce the notion of regular primary mesh that should play a central role in the sequel.

Definition 6. The set defines a regular mesh system if the following conditions are fulfilled.
There exist and , both of them mesh independent, such thatwhere stands for the set of primary cells and where is the distance between a cell-point or a vertex from and an edge-point from .

Proposition 7. Assume that (i) the primary mesh is compatible with the discontinuities of the permeability tensor (see Definition 1), (ii) the set is a regular mesh system, and (iii) relation (8) is honored.
Then, there exists a strictly positive number mesh independent such that

Using the previous notations and thanks to the consistency of the flux approximation across cell edges (see Proposition 7), one reasonably can approximate the flux balance within any primary cell as follows:where , , and are such that and where . It is clear that the number of discrete unknowns and is greater than the number of equations in system (14). It is then natural to complete this system with discrete equations obtained from mass flux balance over dual cells; see Figure 4 for the definition of the dual mesh.

In what follows, we will need to use the notion of pseudoedge associated with dual cells. Let us define now this notion that is illustrated in Figure 5.

Definition 8. Let and be two cell-points from the primary mesh (i.e., ) such that the corresponding primary cells and are adjacent, and consider (recall that , for , is the set of edge-points lying in the boundary of the primary cell ). The line defines a pseudoedge denoted by and “centered” on , with extremities and .

We will say that a pseudoedge is associated with a dual cell if it is lying in the boundary of .

Remark 9. The boundary of each dual cell is a union of a finite number of pseudoedges (see Figure 4).

Let us now look for discrete flux balance equations over dual cells. Integrating the two sides of the balance equation (1) in a dual cell , applying Ostrogradsky’s theorem for the left-hand side, and exploiting Remark 9 lead towhere stands for the outward unit normal vector to the boundary of and where is a pseudoedge associated with the dual cell . Recall that is the set of edge-points lying in the boundary of the dual cell .

Let us look for a flux approximation across the pseudoedge viewed as part of the boundary of . Denoting by the exact flux over , it can be expressed by the relation

By using the same process, the computation of the flux across iswhere (recall that denotes the angle defined by the vectors and ).

Similarly, the computation of the flux across leads towhere We should now formulate a global estimate for the truncation errors associated with the flux approximation over the pseudoedge.

Proposition 10. Under the same assumptions as those of Proposition 7, there exists a positive number mesh independent such that

The above inequality shows the consistency of the flux approximation across the pseudoedge . We summarize what precedes aswhere we have setWe deduce from (22) that an approximate flux balance equation over any dual cell can be formulated as follows:where , , and are such that and and where , a boundary-vertex are such that and .

Systems (14) and (24) naturally suggest defining a finite volume formulation of the diffusion problem (1)-(2) as follows.

Find and such that

Remark 11. It is useful to note that

2.2. Existence of Solutions to System (25)-(26) Conditions for Uniqueness

(i) Matrix Properties of the DDFV Problem (25)-(26). It is easily seen that the symmetry of the matrix associated with the linear system (25)-(26) essentially follows from the symmetry of the diffusion coefficient (see assumption (3)). We assume that all the cell-points and all the vertices (with respect to the primary mesh) are numbered. Therefore one can identify and with two disjoint subsets of the set of positive integers. To fix ideas, let us setThen denotes the total number of cell-points and stands for the total number of vertices. On the other hand, define the subvectors and byand set wherewhere is a subvector with components defined by the right-hand side of (25) and is a subvector with components defined by the right-hand side of (26).

Remark 12. The matrix satisfies the following properties:(1);(2);(3).

Let us introduce two vectors of named and and defined by

Proposition 13 (characterization of Kernel space of ). (i) The matrix is singular.(ii) Moreover, let denote the subset of defined as follows: (named Kernel space of in the sequel); then we have

Sketch for the Proof(i)The singularity of is an immediate consequence of Remark 12.(ii)Define the space byand endow it with the seminorm defined asThen, find the above identification of the Kernel space of using the following Lemma.

Lemma 14. The matrix is positive; that is,

Proof. Let , where and .
It follows from what precedes thatDefinewhere are two adjacent primary cells sharing the edge as interface containing the edge-point . Then, relation (37) becomeswhere denotes the subset of made of primary cells adjacent to the domain boundary . Let us prove that the homogenized symmetric matrix is positive definite; that is,Settingit is easy to check thatwhere we have setwhich are strictly positive numbers.
Since the diffusion matrix is symmetric and positive definite (see assumptions (3)-(4)), Cauchy-Schwarz’s inequality for the inner product associated with ensures thatas either and or and are not collinear. Therefore, and thus is a symmetric and positive definite matrix.
It follows from what precedes that the matrix possesses strictly positive eigenvalues. Let be its least eigenvalue. So we haveRemarking that there exist two real numbers and strictly positive depending exclusively on the geological structure of the domain such thatthe following result is easily seen.
Lemma  15.Moreover Thanks to inequality (45) and to Lemma the proof of Lemma 14 ends.

Lemma 16. The matrix satisfies the relationwhere is a strictly positive number mesh independent and where is a seminorm defined on by relation (35).

Proof. It is conducted with the same arguments as those developed in the previous proof, except that one should go much farther by proving that there exists mesh independent such thatThe eigenvalues of the symmetric positive definite matrix satisfy the so-called characteristic equation associated with ; that is,The least eigenvalue of denoted bywhere , is a strictly positive number. One can easily deduce thatwhere is a strictly positive number. We should bound the quantities , , and by mesh independent strictly positive numbers. Let us start first with We consider a change of coordinates by moving from the initial Cartesian coordinates to a local one, namely, , where is the edge-point on the interface between the cells and and where is a vector orthogonal to and oriented such that the basis change matrix is a rotation. Denoting the permeability tensor of the cell by in the initial Cartesian coordinates and by in the local coordinates we haveSimilarly, we get for the cell Then it is easy to check thatthat is,where denotes the determinant. Similarly, we have for the cell It follows from what precedes thatwhere and are given, respectively, by relations (43).
On one hand, we can deduce from (4) and Definition 6 thatOn the other hand, we remark thatwhere the set (introduced in (8)) depends exclusively on the lithologic structure of the medium . Then we deduce from (59), (60), and (61) thatRemarking thatand exploiting again (4) and Definition 6 lead to the following inequality:Thanks again to (4) and Definition 6 one can easily check thatLemma 16 follows from the combination of (53), (62), (64), and (65).

Proposition 17 (discrete compatibility condition). The right-hand side of the discrete system (25)-(26) satisfies the following discrete compatibility condition:

Proof. First of all, note that the double summation in the previous proposition is equal to zero thanks to Remark 11. Let us consider two vectors of called and and defined byAccording to the matrix form of the discrete system (see relation (31)), we haveIt follows from Remark 12 thatThereforeSimilarly we haveThis ends the proof of Proposition 17.

An immediate consequence of the previous result is the following existence result for a solution to the discrete problem (25)-(26).

Proposition 18 (existence result). (1) The linear system (25)-(26) possesses an infinite number of solutions.
(2) More precisely, if is a solution to the discrete system (25)-(26), then is the set of all solutions for this system.

The question to know how to get the physical solution of the model problem has a natural answer; that is, one should follow the same way as the continuous problem analysis. Since the dimension of is , required are two (linearly independent) discrete versions of the null average condition (7). This naturally leads to the following result.

Proposition 19 (uniqueness result). Under the following (null average) conditions,the discrete problem (25)-(26) possesses a unique solution.

Proof. Let us consider the spaceOne knows from Lemma 14 thatThen it follows that, for all in the subspace of made up of such that conditions (72) are fulfilled, one has .

3. Stability and Error Estimates

We deal here with a theoretical analysis of the solution for the discrete system (25), (26), and (72). Recall that the existence and uniqueness of that solution (under explicit conditions) is proven in the previous section. We assume in what follows that the primary mesh is regular in the sense of Definition 6.

3.1. Preliminaries and a Stability Result

Let us consider the auxiliary mesh introduced in the previous section (see Figure 4). Note that each mesh cell of contains either one cell-point or one corner point and only one which should be lying inside or on the boundary of . In the sequel, a node is a cell-point or a corner point with respect to the primary mesh.

Definition 20. An auxiliary mesh cell is degenerate if the corresponding node is lying on the boundary of

In the sequel, we will say sometimes “auxiliary cell” instead of “auxiliary mesh cell.” We denote by ) the subspace of made of functions which satisfies conditions (72). This space is obviously nonempty as there is the null function. Moreover, the solution of the discrete system (25), (26), and (72) could clearly be identified with one (and only one) function from ). In the sequel we denote by such a function named “cellwise constant (approximate) solution” of the diffusion problem (1)-(2). Let us endow with the following discrete seminorm.

For all , definewhere “s.t.” means “such that.” As this seminorm is actually a discrete energy norm for the space ), we denote it by as far as functions from ) are concerned. A norm on the space is defined by the mapping We focus here on the proof of the stability of the solution to the system of (25), (26), and (72). For this purpose, we need to introduce as in [31] a preliminary result, namely, a discrete version of Poincaré inequality which is based upon the following ingredients. Consider the following spacesand the linear operators and defined as follows: with and , where and are, respectively, subspaces of and made of functions satisfying conditions (72)-(i) and (72)-(ii), respectively.

Let us equip the function spaces and with the following discrete energy norms:These mappings are only seminorms, respectively, for and that one can transform into norms for these spaces as follows:where we have set and . The following results are key ingredients for proving the stability of the discrete solution in the sense of the discrete energy norm defined previously on the space .

Lemma 21 (discrete versions of Poincaré inequality). We have the following inequalities:where represents diverse strictly positive numbers mesh independent.

Note that the right-hand side of inequalities (81) and (82) is, respectively, norms for and . Since is bounded, the previous lemma permits seeing that these norms are equivalent to standard ones, namely, (80). For proving the preceding lemma one can use the same arguments as in [36].

Lemma 22 (a key-result). Let be a regular mesh system (defined on) in the sense of Definition 6 and denote by the corresponding dual mesh. For , denote by the value of in the control volume . Let be a discrete trace function defined a.e. (for the ()-Lebesgue measure) by on , for all , where denotes the set of mesh elements adjacent to the domain boundary . ThenSimilarly, we havewhere stands for diverse positive numbers mesh independent.

Proof. See, for instance, [36].

Remark 23. It is more than useful to note that

Let us give now one of the main results of this section.

Proposition 24 (stability result). The cellwise constant approximate solution of the diffusion problem (1)-(2) satisfies the following inequality:where is a strictly positive real number not depending on the spatial discretization.

Proof. Multiplying (25) by and (26) by and summing over and , respectively, lead toRecall thatLet us set We know from Lemma 16 that there exists a strictly positive number mesh independent such thatOn the other hand, according to Remark 11 we haveCauchy-Schwarz’s inequality leads to where and are, respectively, discrete trace operators associated with grids and . Thanks to relation (4), Definition 6, and Cauchy-Schwarz inequality we getAccounting with the fact that satisfies null average condition on the grids and , Lemma 21 ensures thatwhereas Lemma 22 guarantees thatA double application of Cauchy-Schwarz’s inequality leads toCombining relations (93), (94), (95), and (96) on the one hand and using the inequality of Cauchy-Schwarz and Remark 23 on the other hand yieldRecall that represents diverse strictly positive numbers mesh independent. The stability result follows from relations (90) and (97). The proof of Proposition 24 ends.

3.2. Error Estimates

Following the original technique exposed in [28], we investigate in this subsection the error estimates for the finite volume approximate solution to the model problem. Recall that denotes the exact solution to the model problem while is the cellwise constant approximate solution obtained from a DDFV formulation of this model problem. View also as a representation in terms of function of the vector , where and . Recall that the vector is defined above as the unique solution of the DDFV discrete system, satisfying null average conditions on grids and (see (25)-(26) and (72)). Let us setNote that the components of the so-called error-vector can be viewed as values in auxiliary cells of the so-called error-function defined a.e. in . More precisely let us recall that any auxiliary cell is associated with a unique cell-point or a unique vertex. Let be the error corresponding to the auxiliary cell ; that is,Since the auxiliary mesh define a partition of , one can setwhere is the characteristic function of any cell .

Our first purpose is to show that the error-vector is a solution to a square linear system of the same type as the discrete system (25)-(26) in the sense that both of them are associated with the same matrix defined above. Note that unfortunately is not in as it does meet the null average conditions (72) on grids and . So the error estimates, which are our main purpose in this section, should be investigated in the sense of the norm defined by (76) on the space . In this connection we first should prove thatand in the next step we should prove that

(i) Let us start with proving that the components of the error-vector are solution to a square linear system of the same type as the discrete system (25)-(26).

Accounting with truncation errors, one can see from the DDFV flux computations that the exact flux across the edge viewed as part of the boundary of the primary cell is given bywhere the truncation error is given byand meets the following inequality (see Proposition 7):On the other hand, the exact fluxes across the pseudoedge for and across the edge for , both of them acting as part of the boundary of the dual cell , are given by the following relations (accounting with truncation errors):where the truncation errors and meet the following inequalities (according to Propositions 7 and 10):The system of equations satisfied by the exact nodal potentials and reads aswhere is defined in (23). Recall that is the set of neighboring cells of a given (primary or dual) cell.

Due to the conservation of our flux approximation schemes, the truncation errors naturally obey to the following relation:where and are two cell-points or corner points associated with primary or dual adjacent cells. Note that, in the previous equality, the truncation errors are written in a very formal way. This is with a view to involving both primary and dual adjacent cells. Let us introduce now the set of diamond cells called in what follows the diamond mesh and denoted as (the concept of diamond cell has been introduced for a different usage in earlier works of some authors: see, e.g., [22, 25, 34]). Each diamond cell is associated with one edge-point and vice versa (see Figure 6). A diamond cell is declared degenerate if the corresponding edge-point is lying on the boundary of . An example of a degenerate diamond cell is provided in Figure 6.

The following assumption plays a key-role in our proof of the convergence of the piecewise constant (approximate) solution .where is the Lebesgue measure in any spatial dimension.

Recall that we have set , for , and , for . An adequate linear combination of equations from system (108)-(109) with those from system (25)-(26) shows that the quantities and satisfy the following equations:

(ii) Let us now prove that , where is a mesh independent positive number.

Multiplying (112) by and (113) by and summing over and yield (thanks to Lemma 16)where is some strictly positive real number mesh independent.

Thanks to (109)-(110) (consequence of the conservation of the flux approximation scheme) we haveThereforewhere we have setIt is important to mention that, due to (105) and (107), we haveDefine Lebesgue measure of the diamond cell defined by the points and associated with and Lebesgue measure of the diamond cell defined by the points and associated with .

Therefore, thanks to Cauchy-Schwarz’s inequality, it follows from (116) thatOne concludes with the help of assumption (111) and relation (118) thatwhere is a strictly positive constant that is mesh independent.

(iii) Let us now prove that .

We start with settingRecall is the auxiliary grid and is a cell-point while is a vertex (with respect to the primary grid introduced in a preceding section).

DefineNote that and are a partition of in the sense thatTwo main steps are necessary in our technique to get the estimates of ; before exposing these steps, we develop some preliminaries.where and are operators from the function space into the function spaces and , respectively (see (78) for the definition of and ).

Let and be two real numbers and and two real functions defined almost everywhere in such thatDefineIt is then clear that the vectoris identifiable with a function denoted again by which lies in the space and thus in the space

Proposition 25. The function meets the following trivial properties:where stands for diverse positive numbers that are mesh independent.

Let us look for the estimates of :

(i)  Since the primary cells (with ) are convex and the function the Taylor-Lagrange expansion applies and gives rise to what follows: Since meets the null average condition (7) and aux honors condition (125)-(ii), integrating the two sides of (134) in a primary cell and summing on yieldwhen accounting with the fact that a.e. in , where represents diverse positive numbers that are mesh independent.

(ii)  Thanks again to Taylor-Lagrange it is easily seen that

(iii)  At last, we haveSummarizing what precedes we have proven the following.

Lemma 26. One has

Remarking that any dual cell (with ) is actually a finite union of convex homogeneous polygons sharing as a common vertex, the arguments that have led to the preceding lemma apply and give rise to the following.

Lemma 27. One has

Let us summarize the previous developments in the following proposition in which denotes the set of nodes (i.e., cell-points and vertices) with respect to the mesh system defined on .

Proposition 28 (error estimate result). Assume that the primary mesh is regular in the sense of Definition 6 and that the discontinuities in of the piecewise constant permeability tensor generate a finite number of convex subdomains over which the exact solution of (1)-(2) meets the following property:Under conditions (3), (4), (6), (11), and (111), the error-function associated with the error-vector from with components satisfies the following estimate:where represents mesh independent real number.

4. Numerical Test

The triangular coarse and fine meshes (see Figure 7) are also from FVCA5 Benchmark (see, e.g., [37]). Let us consider a diffusion problem formulated as (1)-(2), where the permeability tensor is defined by the following relation:The exact solution is taken to be . Furthermore, we shall ensure the uniqueness of the solution to (1)-(2) by enforcing the conditionNotations(i) : number of unknowns.(ii) : number of nonzero terms in the matrix.(iii) : the discrete flux balance; that is, , where , , , and are, respectively, the outward numerical fluxes through the boundaries , , , and (e.g., flux0 is an approximation of ), and , where denotes some point of the control volume ; note that the residual is a measure of the global conservation of the scheme.(iv) : value of the minimum of the approximate solution.(v) : value of the maximum of the approximate solution.(vi) , : approximations of the energy following the two expressions:Let us denote by the exact solution, by the mesh, and by the piecewise constant approximate solution.(vii) : relative discrete -norm of the error; that is, for instance,(viii) : discrete -norm of the error on the gradient; that is, for instance,(ix) : for (x) : for , same formula as above with ergrad instead of erl2.(xi) , , , : relative error between , , , and the corresponding flux of the exact solution:(xii) : order of convergence of the method for -norm of the solution as defined by erl2 with respect to the mesh size:where is the maximum of the diameter of the control volume.(xiii) : order of convergence of the method for the norm defined by Remark 23.(xiv) ocvgradl2: order of convergence of the method in -norm of the gradient as defined by ergrad with respect to the mesh size, same formula as above with ergrad instead of erl2.Comments about the Numerical Test Results. The results of numerical computations of pressure and its gradient show a convergence of order two in -norm and a convergence of order 1.4 in discrete energy norm (see, e.g., Table 1). Moreover the similarity of curves erl2 and ergrad (Figure 8) confirm the closeness of their order of convergence. There is no discordance with the theoretical result where a linear convergence is announced (see Proposition 28). Indeed, the presence of discontinuities in the permeability tensor coefficients prevents the exact solution from being regular enough in the whole domain. More precisely in presence of such discontinuities, the exact pressure is in and never in no matter how regular may be the data and the domain boundary. So a much slower convergence should take place.

Competing Interests

The authors declare that they have no competing interests.