Fixed-point algorithms for stationary flow in porous media

COMPUTER METHODS NORTH-HOLLAND

IN APPLIED

MECHANICS

AND ENGINEERING

56 (1986) 17-45

FIXED-POINT ALGORITHMS FOR STATIONARY FLOW IN POROUS MEDIA

L.D. MARINI

and P. PIETRA

Istituto di Analisi Numerica de1 C.N. R., Corss C. Albert0 5, 27100 Pavia, Italy Received 13 February 1984 Revised manuscript received 7 June 1985

We study from the theoretical and experimental points of view some fixed-point algo~thms for solving steady-state filtration problems through porous media. The method used applies to a wide class of problems (no restrictions are made on the geometry of the domain or on the boundary conditions; capillarity effects can be treated, the medium can be inhomogeneous and anisotropic). It does not require any previous knowledge on the free boundary, it works on a fixed domain, and the computational cost is slightly more than the cost of solving a linear problem with the overrelaxation method. Moreover, it can be easily implemented in a finite element code. Many numerical tests are presented.

0. Intr~uction

We consider here the steady flow of an incompressible fluid through a porous medium. It is known that the mathematical formulation of the problem leads to fr~c-boundary problems for elliptic equations (see, e.g., [l, 21 for the mathematical modelling). These problems have been studied heuristically for many years by solving a sequence of problems, each of them on a different domain which is adjusted at each step (see, e.g., [3,4]). A great improvement to the theory was made by Baiocchi (see [S, 6]), who proved that in many cases these problems can be transformed into variational or quasi-variational inequalities (see also [7,8]). This approach has a rigorous mathematical basis and gives rise to efficient numerical algorithms (see [9]), However, its application requires some restrictions on the geometry of the domain; also, different geometry and boundary conditions may lead to different formulations. A new mathematical approach was recently introduced by Alt [lo] and by Brezis et al. [ll]. This formulation is more complicated, but allows the treatment of these problems in a very general framework with respect to the geometry of the domain, the boundary conditions, and the constitutive hypotheses on the medium which can be inhomogeneous and anisotropic. Moreover, with Alt’s formulation one can take into account capillarity effects. This kind of formulation has been the starting point for interesting nume~cal developments. In 1121, Alt studied a discretization of the problem and proved existence and convergence results. However, the concrete examples of decompositions are of the uniform type so that the scheme is bound to some kind of finite difference framework; in particular, an approximation of the domain is always made even in the rectangular case. Following Alt’s approach, Pietra presented in [13] another scheme using conforming finite elements coupled with upwinding techniques for the discretization of the nonlinear term (see [14]). This scheme permits much more general 00457825/861$3.50

0

1986, Elsevier Science Publishers

13.V (North-Holland)

18

L. D. Marini, P. Pietra, Algorithms for stationary flow in porous media

decompositions and lies in a true finite element framework. In particular, it adapts immediately to any polygonal domain without any changes in the geometry. In the present paper we follow this scheme which leads, like Alt’s, to a fixed-point problem for a suitable monotonic operator. We study, from the theoretical and experimental point of view, some fixed-point algorithms, all of the iterative type. For two of them we prove convergence; we note that they are of the monotonic type so that the exact solution (of the discrete problem) is the common limit of two sequences, one decreasing and one increasing. The other two, of the heuristic type, don’t have this nice property and their convergence is not proved.However, the numerous experimental results were very satisfactory and showed their practical efficiency. All the algorithms, tested in a large number of different physical situations, proved to be extremely effective in dealing with filtration problems. The results obtained show on one hand the goodness of the approximation and prove on the other hand the greater efficiency of the new algorithms proposed here with respect to the original one used in [12,13]. To conclude, we point out that a general program has been implemented, including all these algorithms and applicable in a very general framework. No restrictions are made on the geometry of the domain and on the boundary conditions. The boundary can be arbitrarily split into three parts (not necessarily connected): an impervious one, a part in contact with the air, and a part in contact with the reservoirs. (The only restriction is that this last part must be nonempty.) The program works for inhomogeneous and isotropic media and the implementation of the anisotropic case is presently being studied. The method works on a fixed domain and does not need any initial guess on the free boundary, Moreover, the method is cheap; it costs a little more than solving a linear problem with the overrelaxation method. Finally, its versatility makes it suitable for entering finite element codes; in particular, this program has already been integrated into the Modulef library [15,16]. The layout of the paper is the following: in Section 1 we briefly recall the continuous formulation and in Section 2 the discrete approach. Section 3 is devoted to the presentation of the algorithms, while in Section 4 a large number of numerical tests is presented.

1. The continuous problem We recall here, for the reader’s convenience, the continuous formulation. The cross-section of the porous medium is denoted by R. Its boundary an is divided into three parts: S,, S,, and S,; S, is the part in contact with the air, S, the part in contact with the reservoirs, and S, the impervious one (see Fig. 1). We look for the pressure u = u(x) of the fluid satisfying u = 0 in the (unknown) dry region. Let U > 0 denote the given pressure, hydrostatic plus atmospheric on S,, and atmospheric on S, . From the mathematical point of view, we assume that 0 is a bounded connected open set of R* with a Lipschitz-continuous boundary 80 and that S, has positive length. Consider K={vEH’(R)Iv=uonS,,v~uonS,}, r = { + E L”(a)

IOS I+!J S 1, a.e. in a} .

We shall use the following formulation

of the continuous

problem:

19

L. D. Marini, P. Pietra, Algorithms for stationary flow in porous media

s2 s,

Fig. 1. The continuous

(RESERVOIR)

(IMPERVIOUS

1

problem.

Find (u, r) E K x r such that u2 0 I

a.e. in R ,

y =

1

n(Vu+y*e)rV(u-u)dxaO

a.e. in {x 1U(X)> 0} ,

(1.1)

VUEK,

where x: = (x,, x2), e = (0,l) is the vertical unit vector, and 1 = {T,(X)} (i, j = 1,2) is a second-order symmetric tensor representing the permeability of the porous medium. We assume 3constanta>O: For results concerning [lo, 11,17,18].

[17,1[ja~]t[2 existence,

V[ER2.

uniqueness,

and regularity

of this problem

we refer to

REMARK 1.1. In formulation (l.l), u = U(X)is the pressure of the fluid and y = y(x) turns out to be, in most practical cases, the characteristic function of the set {x ( U(X)> O}. Nevertheless, in some situations the region {x I 0 < y(x) < 1) is not empty; this fact may be interpreted as the presence of unsaturated regions in the physical problem. (See also Section 4.) The choice of the permeability coefficients allows to treat the general case of an REMARK 1.2. inhomogeneous anisotropic medium. Moreover, formulation (1.1) permits to take into account capillarity effects by simply assuming the atmospheric pressure to be strictly positive. In order to neglect these effects, the atmospheric pressure must be assumed equal to zero, as usual.

2.The

discrete problem

For the discretization of problem (1.1) we follow the scheme introduced by Pietra [13], which we refer to for a mathematically rigorous formulation. In [13], the basic idea is to approximate the linear term by means of classical finite element methods, while for the discretization of the

20

L.D.

Mar%,

P. Pietra, Algorithms for stationary flow in porous media

nonlinear term the upwind techniques proposed by Tabata [14] are used. More precisely, we make the following assumptions: 0 is a polygon and { Th}h is a family of regular decompositions (see [19]) of 0 into triangles such that VTETh

and

Vangle8ofT:

0<8~in-8,,

(2.1)

where 8, is defined by: Oa&+n

,

cos 8, =

2 d2+h+l/A

with b, = centroidof

A = ]]~(b~)-~]] ]]~(b~)]] ,

Note that for an isotropic medium condition

T

.

(2.1) becomes (since A = 1 in this case):

Moreover, it is possible to treat problems with weak anisotropy (the limit case is A = 3, which implies 0 = i IT). Let N be the number of vertices of T,,. Define: 2 := {Pi E fi, i E I} = {set of vertices of T,,} ;

1:={1,2,...,N};

I, := {i E 2 ( Pi E S,} ;

Z,:={iEZlP,ES,}; ui:=u(Ppo,

iEZ,UZ,.

We denote by Tci, the upwind triangle, associated with Pi, defined by Pi is a vertex of T(,,, Tci, . Pi intersects the oriented half-line with end point Pi and direction r(Pi)e (see Fig. 2 for the isotropic case).

1

e direction

Fig. 2. Definition of Tci,.

(2.3)

L.D.

21

Marini, P. Pietra, Algorithms for stationary flow in porous media

We note that for some points Pi on the boundary do it may happen that the upwind triangle does not exist; on the other hand, a point Pi may have two upwind triangles. In such a case one of them can be arbitrarily chosen. Let 4’(x) (i E Z) be the continuous piecewise linear functions defined by function values at the triangles vertices. We set si:=supp~‘={xEn 2

aij :=

j&>O})

iEI,

j- V&b,)V+‘dx,

TET,,

(2.5)

i, ZEZ,

T

$ meas(Sj)(e~(Pi)V~i)]T~j, eij :=

(2.4)

,

if Tcj, exists ,

(2.6)

otherwise .

0,

We can now write the approximate

problem as follows:

Find (u, 7) E RN x R” such that ujao OSyiGl

C

ui = ui )

ViEZ; ViEZ;

UiSUi)

iEZ,;

iEZ,;

if ui > 0 ;

Yi = l9

(2.7)

(aij + e,)(u, - ui) 3 0

jEI

ViEZ, It

VuERv,

ifEZ,;

s.t.ui=zii,

can be easily proved that problem

vi --L ui )

iEZ,.

(2.7) can be written in the equivalent

form:

Find (u, y) E RN X RN such that uigO OGyiS1

ui = ui )

ViEZ; ViEZ;

1

yi=l,

iEZ,;

ui s ui )

iEZ,;

ifui>O;

(2.8)

C Caijuj + eijYj) { 2: I

/El

The coefficients aij and eij defined in (2.5), (2.6) are such that aii > 0 ,

aij so

eii B 0 ,

e, d 0

)

(2.9) (2.10)

(as proved in [13, Lemma 3.11). Conditions (2.9), (2.10) are sufficient to guarantee the existence of a solution (u, y) of problem (2.8), as proved in [12]: in [12], problem (2.8) is transformed into a fixed-point problem which has a solution if conditions (2.9), (2.10) are satisfied.

L. D. Mark,

22

P. Pietra, Algorithms for stationary flow in porous media

Since all the algorithms presented in Section 3 are fixed-point algorithms, we shall show here, for the reader’s convenience, how problem (2.8) can be reduced to a fixed-point problem. Three situations have to be considered. It follows from (2.8) that Case 1. iEZ~Z,~Z,. aiiui + eiiyi = Ci(u, y) ,

(2.11)

where Ci(u, y) = - C aiiuj - 2 eijyj j#i

We try to solve equation ‘i =

.

j#i

(2.11) by choosing yi = 1 so that we have

$ [‘i(‘Y

Y) -

eii]*

(2.12)

II

If ui > 0, the pair (ui, yi) satisfies (2.8); otherwise, that is if Ci(u, y) - eii < 0, it follows from (2.9), (2.10) that eii > O.since Ci(u, y) 3 0 so that we can solve (2.11) by putting ui = 0, ui=O,

(2.13)

and the pair (ui, yi) still fulfills conditions (2.8). Case 2. i E II. Again we try to solve (2.11) with the trial pair (ui, yi) defined in (2.12); if ui satisfies ui G zii, we have found a solution of (2.8), otherwise we take the pair ui = iii

(2.14)

Yiz17

)

which still satisfies (2.8). Case 3. i E Z2. In this case the condition ui = Ui must be verified so that we can choose the pair (2.14) as a solution of (2.8). To summarize all these steps, we are looking for a pair (u, y) such that: ui20

ViEZ;

OGyiGl

ViEZ;

Yiyi=l 7 if ui >O,

(2.15)

and such that (2.16)

(u, Y> = F(w Y)

with F defined in (2.12)-(2.14). The study of the existence of a fixed point for F can be carried out in an easier way by introducing a new vector w: w=u+y~O. It is easy to check that, if we know w, we can find the pair (u, y) which satisfies (2.15) uniquely: q(w) = max(w, - 1,0) ,

yi(w) = min(wi, 1))

i E I.

L. D. Mark,

P. Pietra, Algorithms for stationaryflow in porous media

23

We can then rewrite (2.16) in terms of w. Define

h q(w)

=‘.

if C,(w) 3 e,, 3 0,

u .

Problem (2.13,

[C,(w) - eii] + 1 ,

+_ w9 II

if C,(w) < e,, .

7

(2.16) reduces to the following fixed-point problem:

Find WE@:= w=Aw

{vEIR~~v~z=O

ViEI}suchthat

(2.17)

,

where the operator

A : rWt+ rWt is defined by

Bi(w)

(Aw)~ =

iEZLZ,.Z,,

9

Ui + 1 )

ieZ,,

min(Bi(w), Ui + 1))

iEZ,.

It can be checked that hypotheses

(2.18)

(2.9), (2.10) imply the monotonicity

of A, that is (2.19)

with obvious notation z2vWzi2ui

ViEZ.

The existence of fixed points for A is then guaranteed if there exist a supersolution subsolution for A, that is, two vectors wf, w! such that:

The existence of such vectors is discussed in [12], where explicit formulas for computing are given: w$=O

Vi,

and a

them

(2.20)

wti = 1 + Sup(x’e + Ui) - x’e , iEI

with xi = (xf, xi) = coordinates REMARK

of Pi, i E I.

2.1. In a general case the continuous

problem (1.1) does not have a unique solution

L. D. Marini, P. Pietra, Algorithms for stationary flow in porous media

24

(see [17,18]), and the same situation occurs for the discrete problem. We shall speak in general about a maximal solution wf and a minimal solution w- of the discrete problem. This is also confirmed by some of the numerical experiments of Section 4, although in most practical cases these two solutions coincide.

3. Numerical algorithms for solving the discrete problem (2.8) We shall present in this section four different algorithms for solving problem (2.8). They are all iterative schemes for searching for the fixed points of problem (2.17). For two of them we prove convergence theorems, while the last two are simply heuristic schemes; nevertheless, they turned out to be very effective in all the numerical tests. We also report, at first, for the reader’s convenience, the original algorithm used in [12,13], which is the simplest and the starting point for the faster algorithms of this section. 3.1. The ‘Jacobi’ algorithm The most natural way of solving (2.17) is to apply iteratively the operator A defined in (2.18). The scheme is the following: If wt, w! are given by (2.20)) then for IZa 0, w: and wlT being known, we compute “+I = Aw= , W+ Convergence

w,+*

(3.1)

= Awl.

of algorithm (3.1)

We can state the following theorem. THEOREM 3.1. Let {WI} and {w”_} be the sequences defined by (3.1). and a w- such that:

pm&w: =

W+~W->o,

w+ )

Then there exist a W’

lim wr = w- ,

(3.2)

n--rm

and w+ (w-) is the maximal (minimal)

solution of problem

PROOF.

(2.19) imply (since wt > WY) that:

The monotonicity w” s

w:+l:=

Aw:

properties

< w:+l:= Aw:
,

(2.8).

n>O.

Hence, ViEI the sequence {(wl),} is monotonic nondecreasing and bounded from above, and the sequence {(w :)i} is monotonic nonincreasing and bounded from below. Therefore, the result (3.2) holds. 3.2.

The ’ Gauss-Seidel’

The convergence

algorithm

of algorithm (3.1) is quite slow; by analogy with linear problems the first

L. D. Ma&i, P. Pietra, Algorithms

for stationaryflow

in porous media

25

idea in order to accelerate the convergence is to apply the scheme (3.1) substituting the ‘Jacobi’ operator A by the ‘Gauss-Seidel’ operator G defined by G:R!+iW;,

W+G=GW,

Gtfi= (Gw), := (A$)),

(3.3)

Vi E 1,

where zti) E OwTis given by p

=

w

*

,

fori>&

@j

zyi= i

w, ,,

I

j-=Ci, Pi,

Vjfl.

We can state the following theorem. THEOREM 3.2. Let G be the operator defined in (3.3), (3.4) and A be the operator defined in (2.18). Then the following properties hold:

ULV+GU~GV,

(3.5)

VVER;,

Av~v~Gu~Av~v,

(3.6)

VUE@,

Av~v=$Gv~A~vst~,

(3.7)

VVEIR~,

Av=v($Gv=v.

(3.8)

VU,UER;,

PROOF. Let us first prove (3.5), putting zi = Gu, ti = Cu. From the definition of G we have: Ui = (Ag’“$

,

z$ = (A~(“)~

Vi E I,

with g(l) = u, q(l) = v. Hence l(l) 2 q”‘, which implies for the monotonicity of A: Ul G3lj 1 *

From definition

(3.9)

properties (2.19)

(3.10)

(3.4) we have: (3.11)

Hence, r(i) B q(‘) Vi > 1, from (3.10)-(3.11) and a simple induction argument; therefore, (from (2.19), (3.9)) gi 3 v;-Vi > 1. In order to prove (3.6)) it suffices to note that if u is such that Au 2 u, then 6 = Gu satisfies @i= (AV(~))~ with q(l) = u. This implies U12 u1 and hence rl(0 b v Vi > 1. Thus we get Gj:= (Ar)“‘), 3 (Av)~ 3 vi Vi > 1. Similarly, one can prove (3.7). The result (3.8) is immediate: if Au = v, then from (3.6), (3.7) we have Gv = u; on the other hand, if Gv = v, the definition of G gives qti’ = v and hence

L.D.

26

(Au),

Mark,

= (Aq”‘),

P. Pietra, Algorithms for stationary flow in porous media

=: (Gv), = u,

Vi E I.

We can now write the scheme of the algorithm: If wt, w”_ are given by (2.20) , then for IZ2 0, w! and wl being known, we compute W-‘+’

w;+‘=Gw:, Convergence

= Gw”

(3.12)

.

of algorithm (3.12)

We can state the following theorem. THEOREM 3.3. result holds:

Let (WY } and { wl } be the sequences dejined in (3.12).

lim WY= w+ ,

n--tm

lim wl = w- , n+o:

(3.13)

where w+ and w- are the maximal and the minimal solution of problem convergence is faster than the one obtained with algorithm (3.1). PROOF.

Then the following

The proof is a simple consequence

3.3. The ‘modified Gauss-Seidel’

of Theorems

(2.8).

Moreover,

the

3.1 and 3.2.

algorithm

A faster algorithm can be constructed by analyzing the behaviour of the solution w = u + y of the discrete problem (2.8). In fact, consider the set D(w) = {Pi, i E Z 1 wi > l} = {Pi, i E Z 1u, >O and 7, = l} . At the interior points of D(w), that is, in the set g(w) = {P, E D(w) ) all their neighbours

are in D(w)} A da ,

we get from (2.8) Cazju,=-Ceij=:ei,

P,EZ$w),

jE1

JEJ

or equivalently,

in terms of w,

C aijw, = ei,

jEI

Pi E k(w),

(3.14)

since 2 JEJ

aiJwJ

=

pzc,

aij(uJ

+

‘1

3

ifP,Ed(w)

I

The idea is then to construct

and

taij=O. I

an iterative algorithm,

where at each step the ‘Gauss-Seidel

L. D. Marini, P. Pietra, Algorithm

27

for stationary flow in porous media

iterate Wn+l = Gw” is modified in fi(@““) in order to get a new iterate w”+‘, which satisfies equation (3.14). In doing this, we implicitly assume y”+’ = 1 in &Gn+‘), which implies wnfl(: = un+l + y”+l) > 1

in &G,+‘).

(3.15)

condition (3.15) would be automatically satisfied, if Since WY+’= $:+l> 1 for P, E c?&$““), all the ‘e,’ were positive. Unfortunately, this is not the case in general, so that we have somehow to project on the convex of vectors u 3 1. To make it more precise, let us introduce the operator R:

R:l?i+[W;, where z = Rw is evaluated We first compute

w--+z=Rw,

(3.16)

by the following steps. W = Gw, with G defined in (3.3)-(3.4)

(3.17)

.

Then we identify the set D(W)=

{P,,iE

(3.18)

III+,>l},

and its discrete interior d(w) = {Z’, E D(W) 1all their neighbours

(3.19)

are in D(W)} L an .

Finally, let z E rWTbe the solution of the following minimum problem with constraints: z E IwTsuch that (3.20)

We then assume: (3.21)

z=Rw. REMARK

3.4.The minimum

complementarity

problem

(3.20) can be written in an equivalent

form as a

system

z E lF8: such that zi 2 1

,

(3.22) (2,

-

z, =

l){z,

w,

This formulation

)

a,,z, - el} = 0 ,

VP, E D(w)

VP,@L@). will be used in the sequel.

’

28

L. D. Marini. P. Pietra, Algorithms for stationary flow in porous media

We can now write the scheme of the algorithm: If zy, z! are given by (2.20) , (3.23)

then for II b 0, z: and zl being known, we compute z:+l=

Rz:,

z:+'=Rz"

with R defined in (3.16)-(3.21). We shall prove that the sequences {z:} and {z”} converge to the maximal and the minimal solution of the discrete problem (2.8) faster than the sequences {WY}and {rvl} constructed with the algorithm (3.12). Proving convergence of algorithm (3.23) is more complicated than in the case of algorithms (3.1) or (3.12), because R is not monotonic. However, we shall prove that the sequence {z: } ({ zl }) itself is monotonic nonincreasing (nondecreasing). To this end we need some preliminary lemmata. Consider, for any set of indices K C I,the operator A, : Fly-, IRT defined by (3.24) with A defined in (2.18). It follows immediately that A, is monotonic, that is

usv~A,usA,u

VU,UE[W:, Another

easy consequence

of the definition

Vu E rWt such that

from the monotonicity

properties

(2.19) of A

VKCZ.

(3.25)

(3.24) of A, is the following property:

usAusw,

(3.26)

where w is a fixed point for A then (3.27)

VKCZ.

usA,v~Auav

Moreover, by applying A to the first inequality in (3.27) we get Au S AA Ku which, combined with AK u G Au, gives usA,u~AA,uav

(3.28)

VKCZ.

Respectively, VUEIW:

suchthat

(3.29)

WSAUGU,

we deduce (3.30)

w~Au~A,usu VKCZ, w~AA,u~A,ua.~ VKCZ.

(3.31)

For i E I, let us denote by A' the operator (i=O,..., N) be defined by

z(0) = u ,

z(i) = A'z(i - 1)

For any vector u E rWzlet z(i) A, with K = {i}.

Vi E I.

(3.32)

L.D.

Marini, P. Pietra, Algorithms for stationary flow in porous media

Then it easily follows from the definition

29

(3.3), (3.4) of G that

Gu = z(N) .

From this new form of G we can simply derive other useful properties

of G:

VUEIW~) w~Ausu~w~AGu~Gu,

(3.34)

where w is a fixed point for A. Property (3.33) ((3.34)) comes easily from (3.28) ((3.31)) applied recursively in (3.32), since z(0) = u and at each step z(m) satisfies inequality (3.26) ((3.29)). 3.5. Let A,

LEMMA

G, and R be the operators defined in (2.18), (3.3)-(3.4), and respectively, and let w be a fixed point for A. We have the following result:

(3.16)-(3.21),

VuEEq

such that

u s Av s w ,

(3.35)

then

,,~G,,=:~sRu=:z~Az~w. PROOF. It follows immediately

(3.36) from (3.35) and the properties

(3.6), (3.8) of G that

usvsw. Moreover,

(3.37)

from (3.33) we also get v~AAv
(3.38)

Hence, by applying the result (3.27) to V we obtain from (3.38) (3.39)

VKCI,

I%A,&ASW

where A, is defined in (3.24). Let us choose the set K as follows: K={kEZI

PkEd(v)})

with h(G) defined in (3.19) and consider the following sequence: Z

(0) =

6

,

Zcm) = A,Z@-‘)

The sequence z(“‘) is monotonic Z(m)< Z(m+l) < w

,

m 3

nondecreasing

vmao.

(3.40)

1 .

and bounded

from above, that is (3.41)

L. D. Marini, P. Pietra, Algorithms for stationary flow in porous media

30

In fact, from (3.39) we have z(O)G AKzcO) =: z(l) s w. Then, by induction arguments, (3.41) follows from the monotonicity properties (3.25) of AK and the obvious consideration that A,w = Aw = w. Hence, (3.41) implies that there exists a vector z(“) such that

(9

z(“) =

zCm) ;

lim rn-+rn

(ii)

z(“) = AKz(=) ;

(iii)

vSZ(mkZ(mkW

(3.42) vm>o.

We claim that

z(“) s

Ru

z =

(3.43)

.

As a matter of fact, z(“) b 6 implies zr’ > 1 for each i E K. Then, recalling the definition (3.19) of &fi), z(“) (=zP + 1 in K) is given by (see (3.24) and (2.18)) - 2 LX~,U:“’ + ei) + 1 ,

iE K ,

I+’

(3.44) iGK.

Moreover,

we have

c

a

p

‘I

=

I

$9

a

II1

+

2 ajju~“’+ c aij j#l

ICI

,

iE K .

(3.45)

I+'

Recalling that cj aij = 0 Vi E I, from (3.44) and (3.45) we deduce that z(“) satisfies c a,,zy’ = e, ,

iE K ,

ti”‘>1)

iEK,

i

p)

L

= 5.

17

(3.46)

ieK.

The uniqueness of the solution of the complementarity Combining (3.37), (3.42)(iii), and (3.43) we have

system (3.22) then implies (3.43).

It remains to prove that (3.47)

z
To this end we note that the sequence A K~(m)6 AA,z(“)

s w

(3.40) satisfies

Vm 2 0 ,

(3.48)

31

L. D. Marini, P. Pietra, Algorithms for stationary flow in porous media

by applying (3.28) recursively to z(“‘),and since z(O)= V satisfies (3.38). At the limit in (3.48) we obtain (3.47) and the proof is complete. We can now state the following theorem. THEOREM 3.6. Let {WY}, {z”} be the sequences defined in (3.12) and (3.23), and let w - be the minimal solution of problem (2.8). Then the following holds:

w:
respectively,

(3.49)

VnaO.

PROOF. Since z! = w! satisfies z! s AZ’_ s w- we can apply the result (3.36) of Lemma 3.5 recursively and obtain -fl+1

:=Gz~~t~+~
z” 62 _

VnaO.

(3.50)

By induction we can now deduce (3.49); in fact WY< z _’ s w- holds, and wl < z” d w- implies, using (3.50) and the monotonicity of G: n+l

W-

:=

Gw” s Gz” d z_n+l 6 w- .

REMARK 3.7. A quite natural consideration that follows from Theorem 3.6 is that for the construction of the sequence {zI> the projection on the convex of vectors v 2 1 is not necessary. In fact, from (3.50) we get zlf >i”Vn; _ in particular, since by definition we have ZY > 1 in fi(Zl), it follows that zlf automatically satisfies ~11> 1 in B(Y). Therefore, Vn the complementarity system (3.22) reduces to the linear system (3.46). The situation is different for the sequence {z:}. We shall prove that this sequence is monotonic nonincreasing, so that the analogue of inequality (3.50) will be of type z: s z’: Vn, which does not imply z$ 2 1. Hence, in the construction of the sequence {z:}, the use of the projection on the convex v 2 1 is, a priori, necessary. This implies that the previous techniques have to be adapted in order to prove the convergence of {z y } . Let K C I be a set of indices and let Ai : rWy+ rWTbe the operator

defined by (3.51)

withA defined in (2.18). It follows immediately from the monotonicity Ai is monotone, that is

Vu,uElR~, The following property VUE@

u~u+A;u~A;v is a trivial consequence

suchthat

WSAUSU,

VKCZ.

properties (2.19) of A that

(3.52)

of definition (3.51): (3.53)

we have w
VKCI,

(3.54)

32

L.D.

Marini, P. Pietra, Algorithms for stationary flow in porous media

where w is a fixed point for A. Moreover,

if we associate with any v verifying (3.53) the set:

K=K(u)={iEZ(u,>l},

(3.55)

we also get A&
(3.56)

In fact, by definition (A&), = vi Vi G K, while Vi E K either (A&), = 1
(3.57)

(3.57) can easily be proved by applying A in (3.56) and using (3.54).

Inequality

LEMMA 3.8. Let A, G, and R be the operators defined in (2.18), (3.3)-(3.4), and (3.16)-(3.21), respectively, and let w be a ftxed point for A. The following result then holds: suchthat

VUER~

WGAUSU,

(3.58)

then

wdAz
(3.59)

PROOF. The proof, which is quite technical, follows the same steps as that of Lemma 3.5. From (3.58) and properties (3.7), (3.8) of G we immediately have (3.60)

w
(3.61)

Let us choose the following set of indices: K={iEZj with d(v)

P,e&V)},

(3.62)

defined as in (3.19). From (3.61)-(3.62)

we get (see (3.53)-(3.56)) (3.63)

w~A&A;&rj. Consider now the sequence ztrn) defined by Z

(0) =

6

7

Z(m) = A;z’m-”

The sequence zcrn)is monotonic

,

m 3

nonincreasing

1 .

and bounded

(3.64)

from below, that is

L. D. Marini, P. Pietra, Cm)
wdz

Algorithms for stationary flow in porous media

33

(3.65)

Vmal,

as can be easily seen, by induction arguments, from the monotonicity and the fact that Aiw = Aw = w, and z(O) satisfies (see (3.63))

properties

(3.52) of Ai

Hence, (3.65) im pl ies that there exists a vector z(“) E IRT such that: (i) (ii) (iii)

z(“) = lim zcrn); = z(m) ; A ;zP) m-m

w S z(“;)d z (m)SV

(3.66) Vm.

We claim that z(p) E z = Rv .

(3.67)

Indeed, from the definition (3.51) of Ai we deduce that Vi E K 2,‘“’2 1 and zi”’ 2 (Az@))(. As for Lemma 3.5, we can deduce with similar arguments that z(“) satisfies c a,,zjm) > e, j p

I 2

1

1

ViEK,

(2~“’- 1)(x

aL,zi”’ - ei = 0 I

)

J

The uniqueness of the solution of the complementarity (3.60), (3.66)(iii), and (3.67) we then have

system (3.22) implies (3.67). Combining

It remains to prove that WGAZSZ.

(3.68)

To this end we note that the sequence w s A/@@” s A;z’“’

(3.64) satisfies

VmSO,

(3.69)

by applying (3.57) recursively to z@) and since z(O)= G satisfies (3.61). At the limit in (3.69) we obtain (3.68) and the proof is complete.

L.D.

34

Marini, P. Pietra.

Algorithms for stationary flow in porous media

We can now state the following theorem. THEOREM 3.9. Let {WY}, {z:} be the sequences defined in (3.12) and (3.23), and let wf be the maximal solution of problem (2.8). Then the following holds: W

+<“< --z+--w,

n

VnaO.

(3.70)

Since zt = WYverifies wf s AZ’s recursively and obtain

zy we can apply the result (3.59) of Lemma 3.8

PROOF.

w+ d AZ;+’ < z:+l

s

respectively,

,;+’ := Gz:

Q

z: Vn>O.

(3.71)

By induction we can now deduce (3.70). In fact, w+ < zy s WYholds and w+ d z: s w: implies, thanks to (3.71) and to the monotonicity of G, W+ < ‘Z+

n+l

6

Gz:

< G,,,”

+

=-

*

w:+’

Hence, (3.70) holds. We end with the following theorem. THEOREM result:

3.10.

pinz”

Let {z:},

=

w-

)

{z”} be the sequences defined by (3.23);

we then have the following (3.72)

lim z: = W+ ,

n+=

where wf and w- are the maximal and the minimal solution of problem convergence is faster than the one obtained with the algorithm (3.12). PROOF.

The proof is a simple consequence

(2.8).

Moreover,

the

of the results (3.49) and (3.70) of Theorems

3.6

and 3.9. REMARK

3.11.

The implementation of this algorithm demands the resolution, at each step, of the minimum problem (3.20) (or equivalently of the complementarity system (3.22)). The method we used is the overrelaxation method with projection (SORP) (see e.g. [20,21]); at each fixed-point step the initial guess for the SORP is the just computed vector r? : (W1). (See (3.17).) 3.4.

Some heuristic algorithms

A crucial step in the algorithm (3.16)-(3.21) is the solution at each fixed-point step of the minimum problem (3.20). This can make the global algorithm quite expensive, even if a good solver is used. On the other hand, the solution of (3.20) is an intermediate step within the fixed-point algorithm so that a precise computation is actually unnecessary. These considerations led us to build up two other algorithms, of the heuristic type, which are variants of algorithm (3.16)-(3.21) and are based on the use of the SORP.

L. D. Mark,

P. Pietra,

Algorithms for stationary flow in porous media

35

The first consists in applying the scheme of algorithm (3.16)-(3.21), where, instead of solving problem (3.20), just one SORP iteration is carried out. Of course, no proofs of convergence can be given in this case; nevertheless, numerical experiments show that the scheme is good and converges faster than the previous ones, as we shall see in the next section. Another efficient algorithm, called ‘the incorporated SORP’, emerged from similar considerations after observing the satisfactory performance of the l-SORP iteration algorithm. The idea is basically the same and consists in modifying directly, component by component, the intermediate Gauss-Seidel iterate via SORP. More precisely, the iterative scheme can be described as follows: For 123 0, WV” being known, we compute

ivn+’ by n+l

,;+I

=

(AZ”‘),

with

z(l) = w I’+’andfori>l,

zig)=

wJ

w; 3

-

W;+l =

n+l

WE

Y

max{l, oWl+’ + (1 - w)w~} ,

if ,:+l

’

jci’

jai,

(3.73)

sl+h,

if +r+’ > 1 + h .

Algorithm (3.73) is a kind of global SORP with a relaxation parameter o = 1 in the dry region and o > 1 in the wet region. The classical SORP, with a fixed relaxation parameter o > 1 in the whole domain, does not work very well, neither in theory (no monotonicity) nor in practice, as we shall see in the next section.

4. Numerical tests A general program has been implemented to solve problem (1.1) with the discretization of Section 2. This program includes all the algorithms presented above, and works for general domains 0, general boundary conditions, and inhomogeneous isotropic media. Automatic decompositions can be used, provided condition (2.2) is verified. All the algorithms were first tested on a classical example of a rectangular dam with two levels and impervious bottom (see Fig. 3); the domain was decomposed into 448 triangles and 255 nodes. In the comparative Table 1 the number of fixed-point iterations and the CPU time are reported. Since the problem has a unique solution, the stopping test for all the methods was MaxI(

- (wY)~[ < lop3 - (max wet level) .

(4.1)

As already stated, in the ‘modified Gauss-Seidel’ method, the initial guess for the SORP at the nth fixed-point step was X(O)= WY (x(O)= w ? ) for the supersolution sequence (subsolution sequence). The stopping test for the SORP (separately for supersolutions and subsolutions) was

36

L. D. Marini, P. Pietra,

Algorithms for stationary j?ow in porous media

Fig. 3. Dam with two levels and impervious

Table 1 Comparison

bottom.

of algorithms of Section 3

Method Jacobi Gauss-Seidel Modified Gauss-Seidel l-SORP iteration Incorporated SORP

Fixed-point iterations

CPU time

398 198 70 73 58

2’ 33” 92 1’ lo” 02 1’ 10” 56 48” 96 31”68

In Figs. 4-15 we report the results obtained on many different physical problems. The boundary conditions are plotted as in Fig. 3. In all the pictures, except for Figs. 14 and 15, the level lines of the piezometric head U(X) + e. x = u(x) + x2 are drawn; these are orthogonal to the stream lines. Defining the set L!l = { T E Fh 13Pi E T where U, > 0)) as a discrete ‘free boundary’ the line ~9.0l L ~9flhas been chosen and drawn. Note that, even when in the continuous problem y is the characteristic function of the set {U(X) > 0)) in the discrete problems x can take values in the interval IO, l[. In the pictures small vertical lines indicate the nodes where 0 < yh < 1. In Figs. 4-10 the medium is isotropic and homogeneous and the capillarity effects are neglected. An interesting situation is shown in Figs. 7-9: small changes in the data modify the number of connected components of the wet region. The same example of Fig. 4, but with a capillarity fringe, leads to the solution of Fig. 13. Note that in this case convergence results are not proved for the scheme (2.7) (see [13]), but the method still works in practice. An example of inhomogeneous medium is considered in Figs. 11-12; the boundary conditions are the same as in Fig. 10, while the permeability coefficient is 1 in the upper half and 3 (or 9) in the lower half of the domain. In Fig. 12 the solution has two saturated regions, separated by an unsaturated zone. (We refer the reader to [lo] for a theoretical discussion.) In the last example the solution is not unique; Figs. 14-15 plot the computed minimal and maximal solution, respectively. Between these two ones there are many other solutions, corresponding to different levels of the groundwater lake. In these figures the pressure u is plotted; Fig. 15 shows that the pressure is hydrostatic in the groundwater reservoir.

L.D. Marini, P. Pietra, Algorithms for stationaryflow in porous media

Fig. 4. Mesh of 448 elements

and 255 nodes: R = {(x, Y>10. cx~l.,

O.~ya2.};

31

levelsy,=1.6,y,=0.6.

v _

7

Fig. 5. Mesh of 800 elements

and 439 nodes: fi = {(x, y) IO. s x C 2., 0. cy s 1.8); levels y, = 1.8, y, = 0.4.

38

L.D.

Ma&i,

P. Pietra. Algorithms

Fig. 6. Mesh of 980 elements

Fig. 7. Mesh of 400 elements

for stationary jlow

in porolts media

and 534 nodes: levels y, = 1.15, yz = 0.3.

and 221 nodes: G = {(x, y) 10. G x s l., 0. G y c 1.); levels y, = y, = 1.

L.D. Marini, P. Pietra, Algorithms for stationaryJ’?OWin porous media

.

.

Fig. 8. Mesh of 400 elements and 231 nodes: 0 = {(x, y) ) 0. s x G 1.6, 0. s y G l.}; levels y, = y, = 1.

Fig. 9. Mesh of 400 elements and 231 nodes: R = {(x, y) IO. s x =G1.6, 0. s y s 1.); levels Y, = ,y2 = 1.

39

40

L.D. Murini, P. Pietru,

Fig. 10. Mesh of 808 elements

Algorithms for stationary flow in porous media

and 443 nodes: a = {(x, y) IO. 4 n 4 2., 0. 4 y 4 1.8); levels y1 = 1.8, yz = 0.

Fig. 11, Mesh of 808 elements and 443 nodes: a= {(x, y) IO. c x G 2., 0. S y 6 1.8); levels y, = 1.8, y2 = 0 Permeability coefficient is 1 in the upper half, 3 in the lower half.

L. D. Marini, P. Pietra, Algorithms for stationaryflow in porous media

41

-

/

Fig. 12. Mesh of 808 elements and 443 nodes: fi = {(x, y) IO. c x c 2., 0. ==ysl.B}; Permeability coefficient is 1 in the upper half, 9 in the lower half.

levels y, = 1.8, y, = 0.

Fig. 13. Mesh of 448 elements and 255 nodes: 0 = {(x, y) IO. ax~1.,0.sy~2.};1eve1sy1=1.6,y,=0.6.Same problem as Fig. 4 with capillarity effects.

42

L. D. Marini, P. Pietra,

Algorithms for stationary flow in porous media

Fig. 14. Mesh of 294 elements and 180 nodes: levels y, = 0.9, y, = 0. Minimal solution.

Fig. 15. Mesh of 294 elements and 180 nodes: levels y, = 0.9, y, = 0. Maximal solution.

L.D.

Marini, P. Pietra,

Algorithms for stationary flow in porous media

43

For all the examples the ‘incorporated SORP” algorithm was used with the stopping test (4.1) for the examples of Figs. 4-13, where the solution is unique. For the last example the stopping test was: MFX l(w;), - (wy- ‘),I d 10e3 . (max wet level) , Max I

I(w:),

- (WY-'),I

c 10e3 - (max wet levei) .

For the decomposition of the domain into triangles we used programs of automatic decomposition from the library Modulef [22]. An automatic triangulation does not satisfy in general condition (2.2); however, an easy and cheap postprocessing module can be implemented, as described in [16), in order to obtain a new decomposition satisfying (2.2). An example of triangulation is shown in Fig. 16. Note that the method does not require any a priori guess on the position of the free-boundary Iine, so that no re~neme~t of the mesh has ever been used for any of the tests presented. However, for other problems, a refinement around the expected free boundary could possibly accelerate convergence. (Many other examples of trianguIations can be found in [16], which we refer to for a detailed description of the programs as well.) As already mentioned, the global SORP method is not reliable, in the sense that its behaviour is unpredictable. Tests were performed on three different problems, corresponding to the examples of Figs. 4, 9, and 14, Xn Fig. 17 the number of iterations is plotted against w for the problems of Figs. 4 and 9. In both cases there exists a value w0 of the relaxation parameter (depending on the physical problem and on the number of unknowns) such that, for o s w0 the algorithm converges quite fast while for o > o0 it may either converge very slowly or even run for over 3000 iterations without converging. An example of the behaviour of the method is illustrated in Fig. 18 (problem of Fig. 4, o = 1.75, convergence after 770 iterations). The figure shows, in the x-y plane, the points Q, with coordinates x = w:(P), y = w:“‘(P), where P is a

Fig. 16. Example of triangulation

used for the problem of Fig. 6.

44

L. D. Marini. P. Pietra, Algorithms for stationary flow in porous media

3oooL 800

700

PROELE OF FIG.,7?

_

Fig. 17. Number of iterations

Fig. 18. The behaviour

against w.

of the SORP: problem of Fig. 4, w = 1.75.

L. D. Marini,

P. Pietra,

Algorithms

for stationary

flow

in porous

media

45

given point in the wet region. The segments joining Q, to Q,,, are plotted for the last 270 iterations (with the last 70 magnified). The sequence keeps oscillating in a neighbourhood of a point different from the ‘limit’ point Qx. = (w+(P), w+(P)) and finally reaches Qz only at the last iterations. In general, the lack of monotonicity has serious consequences for ‘big’ values of w: after a small number of iterations the subsolution is greater than the supersolution. In the case of nonuniqueness of the solution (see Fig. 14), both sequences { w’J}, {WY } converged to the maximal solution, and the minimal was never achieved. All the numerical tests were performed on the DPS8/44 Honeywell system of the Centro di Calcoli Numerici of the University of Pavia. References [l] J. Bear, Dynamics of Fluids in Porous Media (American Elsevier, New York, 1972). [2] M.E. Harr, Groundwater and Seepage (McGraw-Hill, New York, 1962). [3] C.W. Cryer, On the approximate solution of the free boundary problems using finite differences, J. ACM 17 (1970) 397-411. [4] M. Muskat, The Flow of Homogeneous Fluids Through Porous Media (McGraw-Hill, New York, 1937). [5] C. Baiocchi, Su un problema di frontiera libera connesso a questioni di idraulica, Ann. Mat. Pura Appl. 92 (4) (1972) 107-127. [6] C. Baiocchi, Studio di un problema quasi-variazionale connesso a problemi di frontiera libera. Boll. Un. Mat. Ital. Suppl. 11 (4) (1975) 589-613. [7] C. Baiocchi, V Comincioli, E. Magenes and G.A. Pozzi, Free-boundary problems in the theory of fluid flow through porous media: existence and uniqueness theorems, Ann. Mat. Pura Appl. 97 (4) (1973) l-82. [8] C. Baiocchi and A. Capelo, Variational and Quasi-variational Inequalities; Applications to Free-Boundary Problems (Wiley, Chichester, 1983). [9] C. Baiocchi, V. Comincioli. L. Guerri and G. Volpi, Free boundary problems in the theory of fluid through porous media: a numerical approach, Calcolo 10 (1973) l-86. [lo] H.W. Alt, Stromungen durch inhomogene porose Medien mit freiem Rand, J. Reine Angew. Math. 305 (1979) 89-115. [ll] H. Brezis, D. Kinderleherer and G. Stampacchia. Sur une nouvelle formulation du probleme de l’ecoulement a travers une digue, CR. Acad. Sci. Paris 287 (1978) 711-714. [12] H.W. Alt, Numerical solution of steady-state porous flow free boundary problems, Numer. Math. 36 (1980) 73-98. [13] P. Pietra, An up-wind method for a filtration problem, RAIRO Anal. Numer. 16 (4) (1982) 463-481. [14] M. Tabata, Uniform convergence of the up-wind finite element approximation for semilinear parabolic problems, J. Math. Kyoto Univ. 18 (1978) 327-351. [15] D. Begis and A. Perronnet, The Club Modulef. A library of computer procedures for finite element analysis, Report Modulef No. 73, INRIA, Roquencourt, France, 1982. [16] L.D. Marini and P. Pietra, Solution of filtration problems through porous media. The programs PIGRA, DAMIAN. PLOTDAM, Rept. IAN-CNR No. 395, Pavia, Italy, 1984. [17] H.W. Ah and G. Gilardi, The behaviour of the free boundary for the dam problem, Ann. Sci. Norm. Sup. Pisa 9 (4) (1982) 571-626. [18] M. Chipot, Probleme de l’ecoulement a travers une digue, Doctorat d’Etat. Universite Pierre et Marie Curie, Paris, 1981. [19] P.G. Ciarlet, The Finite Element Method for Elliptic Problems (North-Holland, Amsterdam, 1978). [20] R. Glowinski, J.L. Lions and R. Tremolieres, Numerical Analysis of Variational Inequalities (North-Holland, Amsterdam, 1981). [21] J.C. Miellou, Methodes de Jacobi, Gauss-Seidel, sur (sous) relaxation par blocs, appliquees a une classe de problemes non lineaires, C.R. Acad. Sci. Paris 273 (1971) 1257-1260. du Club Modulef, [22] 0. Koutchmy, P. Joly and A. Perronnet, Les modules de maillage bidimensionelles Rept. No. 4, INRIA, Roquencourt, 1977.