Benchmark on discretization schemes

for anisotropic diffusion problems

on general grids (December 10th)
Raphale Herbin and Florence Hubert
Laboratoire dAnalyse, Topologie et Probabilits, UMR 6632
Universit de Marseille
39 rue Joliot Curie 13453 Marseille
herbin@[Link], fhubert@[Link]
We present here a number of test cases and meshes which were designed to form
a benchmark for finite volume schemes. We address a two-dimensional anisotropic diffusion
problem, which is discretized on general, possibly nonconforming meshes. In all cases, the
diffusion tensor is taken to be anisotropic, and at times heterogenous and/or discontinuous.
The meshes are either triangular or quadrangular. The results which are expected from the
participants to the benchmark range from the number of unknowns, the errors on the fluxes or
the minimum and maximum values, to the order of convergence (when available).

ABSTRACT.

KEYWORDS:

Anisotropic medium, diffusion process, finite volume schemes, benchmark

1. Introduction
The aim of this benchmark is to provide a number of test cases in order to compare the properties (convergence, robustness...) of existing discretization schemes for
anisotropic diffusion problems using general grids.
In all test cases except test 8, the domain is the unit square. The boundary of the
domain is divided into = D N where Dirichlet (resp. Neumann) boundary
conditions are given on D (resp. on N ).
The considered diffusion problem is formulated as:

(Ku) = f on ,

u=u
on D ,

Ku n = g on N ,

(1)

where K : IR 22 is the diffusion (or permeability) tensor, f the source term,

u
and g the Dirichlet and Neumann boundary conditions, and n denotes the outward
unit normal vector to N .
For each test case, we propose some meshes which will be used for the comparison
between the various schemes. The corresponding data files are given in the different
formats which are explained in the README file of the web site. For any related
question to these meshes or formats please get in touch with one of the organizers
(herbin [Link], fhubert [Link]).
The scheme which is used should be described in the introduction of the paper,
along with the known (mathematically proven) results of convergence, stability, or
error estimates.
In order to facilitate the programming, some FORTRAN subroutines giving the
source terms, the diffusion tensors, the exact solutions and their derivatives (when
available), are given in the file sour es.f90 available on this web site.
One may submit a benchmark paper even if only a partial number of test cases and
meshes are performed. The file [Link] should be used for submission and
display of the numerical results.
Please make sure to register on the web site for the benchmark if you intend to
participate, since all updates will be sent to the benchmark mailing list produced by
the registrations.

2. The tests
Test 1: Mild anisotropy
A homogeneous anisotrotic tensor is considered:
K=

1.5 0.5
0.5 1.5

Test 1.1: mesh1 (triangular mesh), mesh4 (distorted quadrangular mesh) This
first solution is very regular, and is tested first on a regular" triangular mesh and
then on a distorted quadrangular mesh. On this latter mesh, we wish to see whether
oscillations appear and whether the approximate solution remains within the bounds
of the exact solution.

u(x, y) = 16x(1 x)y(1 y), f = (Ku)

D = , N = ,

u
= u|

 Test 1.2: mesh1 (triangular mesh), mesh3 (locally refined nonconforming rectangular mesh) This solution increases at the origin, and therefore we use an nonconforming rectangular mesh to see how the schemes behaves.

u(x, y) = sin ((1 x)(1 y)) + (1 x)3 (1 y)2 , f = (Ku)

D = , N = ,

u
= u|
Test 2: Numerical locking [BAB 92, MAN 07]

K=

1 0
0

Meshes : mesh1 (triangular mesh)

Values of the parameter : 105 , 106 .

2 1/y

u(x,
y)
=
sin(2x)e
, f = (Ku)

D = , N = ,

g = (Ku n)| ,

u dx = 0.

Note that the maximum and the minimum of the solution are located on the boundary, and are more difficult to obtain with the Neumann boundary conditions imposed
here. Since is large, the solution is almost constant in the y variable.

Test 3: Oblique flow

This test case represents a flow with boundary conditions such that the pressure
driven flow would like to go from vertex (0,0) to vertex (1,1), but is impeded by a
heterogeneous anisotropic tensor with high permeability in a direction at 40 degrees
from the horizontal and low permeability in the orthogonal direction. This test case is
inspired by a talk given by I. Aavatsmark in Paris in December 2006 at GDR MOMAS.
After the first publication of this benchmark on the web, I. Aavatsmark told us that in
fact, there are more severe test cases for monotony, which he generously handed out
to us. They are described in Tests 8 and 9 below.
K = R

1
0

R1 ,

where R is the rotation of angle = 40 degrees and = 103 .

The shape of the solution is depicted in Figure 1; it was obtained by a computation

by a hybrid finite volume scheme" (see [EYM 07]) on a fine grid.

Figure 1. Approximate solution on a fine grid for Test 3, oblique flow

We wish to see how the schemes respect the maximum principle. Hence the results
should show the maximum and miminum values of the approximate solution. Since
the exact solution is not known, it is difficult to measure the precision of the scheme
with respect to the values of the approximate solution; hence the outward fluxes should
also be given to compare the various schemes, along with two computations of the
energy given by the discrete counterpart of the formulae:
Z
Z
Ku n udx
(2)
Ku udx, E2 =
E1 =

Note that E1 may only be computed for those methods which include a discrete gradient while E2 can be computed with the boundary outward normal fluxes only. Even
though E1 and E2 should converge to the same value on fine grids, there could be a
noticeable difference between E1 and E2 on the coarsest meshes, and the authors are
encouraged to comment on this difference.
Meshes: mesh2 (uniform rectangular mesh) and a reference mesh.

D = , N = , f = 0,

u
is continuous and piecewise linear on and such that

1 on ((0, .2) {0.} {0.} (0, .2)

0
on ((.8, 1.) {1.} {1.} (.8, 1.)

u
(x, y) =
1

on ((.3, 1.) {0} {0} (.3, 1.)

12
on ((0., .7) {1.} {1.} (0., 0.7)
2

Test 4: Vertical fault

The medium considered here is a pile of anisotropic layers with a fault in the
middle, which leads to a discontinuity of the layers at x = .5. Each geological layer
is meshed with one layer of discretization cells only. A Dirichlet boundary condition
is imposed.
The domain may be decomposed as = 1 2 , with 2 = \ 1 , with
1 = 1 r1 , and
1

= (0.; .5]

r1 = (.5; 1)

4
[

k=0
4
[

k=0

[.05 + 2k .1; .05 + (2k + 1) .1) ,

!

[2k .1; (2k + 1) .1) .

It is described in Figure 2 where 1 is in black and 2 in white.

(0,1)

(0,0)

(1,1)

(1,0)

Figure 2. The computational domain and approximate solution on a fine grid (320
320) for Test 4, vertical fault
As in the case of Test 3, the exact solution is not known, and the expected results
are the same as those of test 3, namely minimum and maximum values, outward fluxes
and the energies E1 and E2 .
The diffusion tensor K is anisotropic and heterogenous, and is given by:
 



102

=
on 1 ,




10
0
K=
, with

 

0

102

=
on 2

103

 Meshes: mesh5 (nonconforming rectangular mesh) see Figure 8, the square

mesh 20 20 denoted by mesh5reg and a reference mesh, for instance the square
mesh 320 320 called mesh5ref .
Boundary conditions:

D = , N = , f = 0,
u(x, y) = 1 x

Test 5: Heterogeneous rotating anisotropy

This test is inspired from [AND 07, LEP 05], and induces numerical locking for
some schemes.
K=

1
2
(x + y 2 )

103 x2 + y 2
(103 1)xy

(103 1)xy
x2 + 103 y 2

u(x, y) = sin x sin y, f = (Ku),

Meshes : mesh2 (uniform rectangular meshes)

Test Boundary conditions:
(
D = , N = ,

u
(x, y) = sin x sin y

Test 6: Oblique drain

This test case represents a situation which is encountered in underground flow engineering where an oblique drain consisting in a very permeable layer concentrates
most part of the flow; this drain is meshed with only one layer of discretization cells.
in the case of a pressure gradient driven transport, as often described in reservoir engineering, it seems important that the discretization cells consist in only one homoneneous material: numerical experiments show that otherwise the solution may be badly
approximated. Here we consider the steady case, but wish to verify that the outward
fluxes are as close as possible to the exact values for the meshes considered here, both
for the conforming and nonconforming meshes.
The domain is composed of 3 subdomains:

1 = {(x, y) ; 1 (x, y) < 0},

2 = {(x, y) ; 1 (x, y) > 0, 2 (x, y) < 0},

3 = {(x, y) ; 2 (x, y) > 0},

with

1 (x, y) = y (x .5) .475,

2 (x, y) = 1 (x, y) 0.05.

We take the slope of the drain = 0.2 and define the exact solution and the source
term by:
u(x, y) = x y, on , f = (Ku),

where the permeability tensor K is such that its principal axes are parallel and perpendicular to the drain:


0
R1 ,
K = R
0
with such that = tan and :

 


102

=
, on 2 ,

10

 


=
, on 1 3 .

101

Meshes : mesh6, mesh7

Boundary conditions:

D = , N = ,
u
(x, y) = x y

Test 7: Oblique barrier

This test case is similar to the Test 6, except that we now have to deal with a barrier,
and the aim is that the scheme should respect this barrier as well as the outward fluxes.
We take the same geometry as test 6 above, with the slope of the drain = 0.2.
We take the exact solution to be

1 (x, y) on 1 ,

1 (x, y)/102 on 2 ,
u(x, y) =

2 (x, y) 0.05/102 on 3 ,

and f = (Ku), where the permeability tensor K is heterogeneous and isotropic:



0
K=
,
0
with :
=

1 on 1 ,

102 on 2 ,

1 on 3 .

X
Figure 3. Parallelogram-shaped domain showing the distances X and Y and the
angle .

Meshes : mesh6
Boundary conditions:

D = , N = ,

1 (x, y) on 1 ,

1 (x, y)/102 on 2 ,
u
(x, y) =

2 (x, y) 0.05/102 on 3 .
Test 8: Perturbed parallelograms [AAV 07]
This test case was given to us by I. Aavatsmark [AAV 07], and is meant to test the
schemes for the violation of the maximum principle within the domain. The domain
is parallelogram shaped, as shown in figure 3. The parameters shown in figure 3
are X = 1, Y = 1/30 and = 30 . The medium is homogeneous and isotropic with
K = Id.
Mesh: mesh8 (perturbed parallelogram mesh) see Figure 10.
Boundary conditions and right hand side:
Z

D = , N = , f = 0 in all cells except cell (6, 6) where

f (x) dx = 1.

cell(6,6)

u
(x, y) = 0 on .

Note that the solution u of this problem should be a function with a maximum in
cell (6, 6), decreasing smoothly to zero towards the boundary. If u shows internal
oscillations or if u < 0, Hopfs first lemma is violated.
Note that cell (i, j) is numbered by i + 11(j 1) in the data mesh8.

Test 9: Anisotropy and wells [AAV 07]

Here, is again the square unit domain = (0, 1) (0, 1). The medium is
homogeneous and anisotropic with




1
0
cos
sin
K = M ()
M
(),
M
()
=
,
(3)
0 103
sin cos
where = 67.5.
Mesh: mesh9, the grid is a square uniform grid with 11 11 cells (Figure 11).
Boundary conditions and right hand side:
- The source density f is zero in all cells.
- The pressure is fixed in two cells, approximating a sink and a source with
fixed pressure:
u = 0 in cell (4, 6),
u = 1 in cell (8, 6).

(4)

- Homogeneous Neumann conditions apply at the outer boundary:

Ku n = 0 on .

(5)

The solution u of this problem should satisfy u [0, 1]. If u has extrema on the
no-flow boundary with u 6 [0, 1], Hopfs second lemma is violated.
Note that cell (i, j) is numbered by i + 11(j 1) in the data mesh9.

3. Expected results
When refined the meshes are numbered i = 1 to ngrid, from coarsest to finest.
The structure of the expected results is given in the file [Link]. For each
value of i, one should provide:
For all runs:
nunkw number of unknowns
nnmat number of nonzero terms in the matrix
sumflux the discrete flux balance, that is:

sumflux= flux0+flux1+fluy0+fluy1-sumf,
where flux0, flux1, fluy0, fluy1 are the outward fluxes at the boundaries x =
0, x = 1, y = 0, y = 1, for example
Z
flux0 is an approximation of
Ku n ds
x=0

P
and sumf = KT |K|f (xK ) where xK denotes some point (which should be precised) of the control volume K.
umin: value of the minimum of the approximate solution.

umax: value of the maximum of the approximate solution.

When the analytical solution is known and the mesh refined:
Let us denote by u the exact solution, by T the mesh and by uT = (uK )KT the
piecewise constant approximate solution.
erl2, relative discrete L2 norm of the error erl2 where:
X

|K|(u(xK ) uK )2

KT
X
erl2 =

KT

|K|u(xK )

12

where xK denotes some point (which should be precised) of the control volume K,
(or a variant of such a norm, to be precised).
ergrad relative L2 norm of the error on the gradient, if available (give the definition of the discrete gradient)
ratiol2: for i 2,
ratiol2(i) = 2

ln(erl2(i)) ln(erl2(i 1))

ln(nunkw(i)) ln(nunkw(i 1))

ratiograd, for i 2,
ratiograd(i) = 2

ln(ergrad(i)) ln(ergrad(i 1))

.
ln(nunkw(i)) ln(nunkw(i 1))

erflx0,erflx1, erfly0,erfly1 relative error between flux0, flux1,

fluy0, fluy1 and the corresponding flux of the exact solution:
R


 flux0 + x=0 Ku n 


R
erflx0 = 

Ku n
x=0

(except for the fluxes at y=0 and y=1 for case test 2 -Numerical locking- because these
are zero,give the value of the approximate fluxes only in this case).
erflm L norm of the error on the meanvalue of the flux through the edges of
the mesh, if available (give the definition of numerical flux (Ku n)T )



Z

 1
(Ku n (Ku n)T ) , edges of T
erflm = max 
||

10

 o vl2 order of convergence of the method in the L2 norm of the solution as

defined by exrl2 with respect to the mesh size:
ocvl2 =

ln(erl2(imax)) ln(erl2(imax 1))

ln(h(imax)) ln(h(imax 1))

where h is the maximum of the diameter of the control volume

o vgradl2 order of convergence of the method in the L2 norm of the gradient
as defined by ergradl2 with respect to the mesh size:
ocvgrad =

ln(ergrad(imax)) ln(ergrad(imax 1))

ln(h(imax)) ln(h(imax 1))

For tests 3, 4, 8,9, the exact solution is not known, the maximum principle and the
overall precision of the scheme can be tested by computing:
umin: value of the minimum of the approximate solution.
umax: value of the maximum of the approximate solution.
flux0, flux1, fluy0, fluy1 outward normal fluxes to the boundaries x =
0, x = 1, y = 0 and y = 1 (useless in case of test 9).
The same values should be computed on a reference fine grid so as to be able to
compare the results.
For tests 3, 4, since f = 0, we can compute :
ener1, ener2: values of the discrete computations of the energies given by the
formulae (2) (if a discrete gradient is available in the case of ener1).
eren: relative error between ener1, ener2:
eren =

|ener1 ener2|
.
max(ener1, ener2)

4. The meshes
Figures of the meshes are given after the references. We provide two different
formats for the meshes:

*.typ1 and *.typ2, which are described in the file README in the directory
Meshes.
The size steps of the meshes are given in the following table:

11

i
mesh1
mesh2
mesh3
mesh 4
mesh 5
mesh6
mesh7
mesh8
mesh9

1
2.50E-01
3.54E-01
3.54E-01
3.29E-01
1.41E-01
1.25E-01
1.25E-01
1.24E-01
1.29E-01

2
1.25E-01
1.77E-01
1.77E-01
1.70E-01

3
6.25E-02
8.84E-02
8.84E-02

4
3.12E-02
4.42E-02
4.42E-02

5
1.56E-02
2.21E-02
2.21E-02

6
7.81E-03
1.10E-02

7
3.91E-03
5.52E-03

Please use the same format as above when entering your results in the tables
of [Link], that is the format: sign integer dot integer ingeger E
sign integer integer. This format is given for instance by the FORTRAN format ES9.2

5. References
[AAV 07] A AVATSMARK I., Tests cases for violation of monotonicity, Communication to
FVCA5 organizers, 2007.
[AND 07] A NDREIANOV B., B OYER F., H UBERT F., Discrete duality finite volume schemes
for Leray-Lions-type elliptic problems on general 2D meshes, Numer. Methods Partial
Differential Equations, vol. 23, num. 1, 2007, p. 145195.
[BAB 92] BABUKA I., S URI M., On locking and robustness in the finite element method,
SIAM J. Numer. Anal., vol. 29, num. 5, 1992, p. 12611293, Society for Industrial and
Applied Mathematics.
[EYM 07] E YMARD R., G ALLOUT T., H ERBIN R., A new finite volume scheme for
anisotropic diffusion problems on general grids: convergence analysis., C. R., Math.,
Acad. Sci. Paris, vol. 344, num. 6, 2007, p. 403-406.
[LEP 05] L E P OTIER C., Schma volumes finis pour des oprateurs de diffusion fortement
anisotropes sur des maillages non structurs, C. R. Math. Acad. Sci. Paris, vol. 340,
num. 12, 2005, p. 921926.
[MAN 07] M ANZINI G., P UTTI M., Mesh locking effects in the finite volume solution of 2D anisotropic diffusion equations, J. Comput. Phys., vol. 220, num. 2, 2007, p. 751771,
Academic Press Professional, Inc.

12

Figure 4. Triangular mesh with acute angles: meshes mesh1_1 (left) and mesh1_4
(right)

Figure 5. Uniform rectangular mesh: meshes mesh2_1 (left) and mesh2_4 (right)

13

Figure 6. Locally refined non conforming rectangular mesh: meshes mesh3_1 (left)
and mesh3_4 (right)

Figure 7. Conforming distorted quadrangular mesh: meshes mesh4_1 (left) and

mesh4_2 (right)

14

Figure 8. Non conforming regular rectangular mesh mesh5

Figure 9. Coarse oblique mesh: mesh6 (Left).Fine oblique mesh for the oblique barrier and drain tests: mesh7 (Right)

Figure 10. Perturbed parallelogram grid mesh8 with 11 11 cells. To visualize the
grid, on this picture the height is 3 times the real height of the grid described in the
test case.

15

Figure 11. square uniform grid mesh9 with 11 11 cells.

16