Optimized interface conditions in domain decomposition methods for problems with extreme contrasts in the coefficients

Journal of Computational and Applied Mathematics 189 (2006) 539–554www.elsevier.com/locate/cam

Optimized interface conditions in domain decomposition methodsfor problems with extreme contrasts in the coefficients

E. Flaurauda, F. Nataf b,∗, F. Williena

aIFP, DTIMA (Division Technologie, Informatique et Mathématiques Appliquées), 1 & 4 Av. de Bois Préau,92852 Rueil Malmaison, France

bCNRS, UMR 7641, CMAP, Ecole Polytechnique, 91128 Palaiseau, France

Received 28 September 2004

Abstract

When the coefficients of a problem have jumps of several orders of magnitude and are anisotropic, many precon-ditioners and domain decomposition methods (DDM) suffer from plateaus in the convergence due to the presenceof very small isolated eigenvalues in the spectrum of the preconditioned linear system. One way to improve thepreconditioner is to use a linear algebra technique called deflation, or very similarly coarse grid corrections. Inboth cases, it is necessary to identify and compute, at least approximately, all the eigenvectors corresponding tothe “bad” eigenvalues. In the framework of DDM, we propose a way to design interface conditions so that conver-gence is fast and does not have any plateau. The method relies only on the knowledge of the smallest and largesteigenvalues of an auxiliary matrix. The eigenvectors are not used. The method relies on van der Sluis’ result on aquasi-optimal diagonal preconditioner for a symmetric positive definite matrix. It is then possible to design Robininterface conditions using only one real parameter for the entire interface. By adding a second real parameter andmore general interface conditions, it is possible to take into account highly heterogeneous and anisotropic media.Numerical results are given and compared with other approaches.© 2005 Elsevier B.V. All rights reserved.

Keywords: Domain decomposition method; Optimized interface conditions; Deflation; Porous media flow; Discontinuous andanisotropic coefficients

∗ Corresponding author. Tel.: +33 1 69 33 46 48; fax: +33 1 69 33 30 11.E-mail address: [email protected] (F. Nataf).

0377-0427/$ - see front matter © 2005 Elsevier B.V. All rights reserved.doi:10.1016/j.cam.2005.05.019

http://www.elsevier.com/locate/cam

mailto:[email protected]

540 E. Flauraud et al. / Journal of Computational and Applied Mathematics 189 (2006) 539–554

1. Introduction

The classical Schwarz method is based on Dirichlet boundary conditions. Overlapping subdomainsare necessary to ensure convergence. It has been proposed independently in [18,25] to use more generalinterface conditions in order to accelerate the convergence and to permit nonoverlapping decompositions.In [18], exact absorbing conditions are used in domain decomposition methods (DDMs). They are op-timal in terms of iteration counts [33] but are practically very difficult to use. In [25], Robin interfaceconditions are proposed. These seminal papers have been the basis for many other works: [2–4] or [12] forHelmholtz and Maxwell problems. The idea to design the interface conditions by solving an optimizationproblem related to the convergence rate of the DDM was apparently first raised in [39]. This optimizationproved to be difficult. By using the relation between interface conditions in DDMs and exact absorbingboundary conditions, the optimization becomes tractable and has been the subject of many works: seee.g., [1,7,8,11,20,26,46]. Such transmission conditions are essential for evolution [13,43] as well and forsystems of equations, see [5] for the Euler equations. The approach in these papers consists in choosinga frozen coefficients approach either at the continuous level and then discretized (see e.g., [11,12,32]), orat the discrete level (see e.g., or [14]). See also [26,35,40] for other approaches. In any case, parametershave to be computed at each interface node.

We propose in this paper to use a novel approach where only one or two real parameters have to be chosenfor the entire interface. The method relies on van der Sluis’ result on a quasi-optimal diagonal precondi-tioner for a symmetric positive definite matrix, see [41]. It is then possible to design Robin interface con-ditions using only one real parameter for the entire interface. By adding a second real parameter and moregeneral interface conditions (similar to the optimized of order two interface conditions [1,20]), it is possi-ble to take into account highly heterogeneous and anisotropic media. This kind of optimization would behazardous with a frozen coefficient approach where a discontinuity can not be taken into account properly.

The typical equation we have in mind is

�

�tP − div(�∇P) = f (1)

with � a highly heterogeneous and anisotropic tensor. As an example, this equation arises in porousmedia flow simulations through Darcy’s law. Typically, P is the pressure, � is the compressibility of theporous medium, �t is the time step in an implicit scheme, � is the intrinsic permeability tensor of theporous media which depends heavily on the lithology under consideration. The contrast in the lithologiescan induce a discontinuity of the permeability tensor of typically four orders of magnitude. Moreover,the tensor is highly anisotropic. The ratio of the permeability in the vertical direction to the one in thehorizontal directions may be very large as well. We shall consider anisotropy ratios of up to four ordersof magnitude. For such difficult problems, it is common to have plateaus in the convergence of iterativesolvers even when using preconditioned Krylov-type methods. This is the case for DDMs and for otheriterative methods as well. This phenomenon can be related to the existence of a few very low eigenvaluesin the spectrum of the preconditioned system, see [42]. One possible remedy is to use the deflation method,see e.g., [28,29,34,38,44]. It is then necessary to know all the eigenvectors corresponding to the “bad”eigenvalues. In the method we propose, only the knowledge of two extreme eigenvalues of some auxiliarymatrix is needed.

More precisely, in Section 2 we define the model problem under study. In Section 3, we substruc-ture the DDM. In Section 4, we introduce the Robin interface condition. In Section 5, we optimize a

E. Flauraud et al. / Journal of Computational and Applied Mathematics 189 (2006) 539–554 541

two parameter family of interface conditions. In Section 6 we show numerical results. We conclude inSection 7.

2. Setting of the semi-discrete problem

We consider a model problem set in an infinite tube � = R × � where � is some bounded open set ofRp for some p�1. A point in � will be denoted by (x, y). Let

L := − �

�xc(y)

�

�x+ B(y), (2)

where c is a positive real-valued function and B is a symmetric positive definite operator independent ofthe variable x. For instance, if p = 2 one might think of

B := �(y, z) −(

�

�y�y(y, z)

�

�y+ �

�z�z(y, z)

�

�z

)(3)

with homogeneous Dirichlet boundary conditions and ��0, c, �y, �z > 0 are given real-valued functionsand (y, z) ∈ �.

We want to solve the following problem:

L(u) = f in �,

u = 0 on ��

by a DDM. The domain is decomposed into two non-overlapping half-tubes �1 = (−∞, 0) × � and�2 = (0, ∞) × �. The problem can be considered at the continuous level and then discretized (see e.g.,[11,12,32]), or at the discrete level (see e.g., [26,35] or [14]). We choose here a semi-discrete approachwhere only the tangential directions to the interface x = 0 are discretized whereas the normal directionx is kept continuous.

We therefore consider a discretization in the tangential directions which leads to

Lh := − �

�xC

�

�x+ B, (4)

where B and C are symmetric positive matrices of order n where n is the number of discretization pointsof the open set � ⊂ Rp. For instance if we take B to be defined as in (3), B may be obtained via a finitevolume or finite element discretization of (3) on a given mesh or triangulation of � ⊂ R2.

We consider a DDM based on arbitrary interface conditions Q1 and Q2. The corresponding OptimizedSchwarz method (OSM) reads:

Lh(un+11 ) = f in �1, Lh(u

n+12 ) = f in �2,

Q1(un+11 ) = Q1(u

n2) on �, Q2(u

n+12 ) = Q2(u

n1) on �,

(5)

where � is the interface x = 0. It is possible to both increase the robustness of the method and itsconvergence speed by replacing the above fixed point iterative solver by a Krylov-type method. This ismade possible by substructuring the algorithm in terms of interface unknowns

H1 = Q1(u2)(0, .) and H2 = Q2(u1)(0, .).


Let us define the operator

T : H1, H2, f −→ (Q2(v1)Q1(v2)),

where vi , i = 1, 2 solves

Lh(vi) = f in �i ,

Qi(vi) = Hi on �.(6)

The substructured problem is obtained, see [30], by matching the interface conditions on the interfaceand reads(

H1H2

)− �T(H1, H2, 0) = �T(0, 0, f ), (7)

where � is the swap operator on the interfaces: �((H1 H2)T) = (H2 H1)

T.At this point, it should be noted that the analysis of the present paper is restricted to rather idealistic

geometries. However, the same formalism can be used for a domain decomposition into an arbitrarynumber of subdomains [12]. It has also been found there that the convergence estimates provided in thissimple geometry predict very accurately the ones observed in practice even for complicated interfaceboundaries.

3. The substructured problem

The convergence rate of (5) and the spectra of (7) depend on the choice of the interface conditions Q1,2.In order to design an efficient method, we need to have a formula for the substructured problem and sofirst for the solution to (6) with f = 0. An essential tool will be the Dirichlet to Neumann map whosesymbol is obtained here via a factorization of the operator Lh.

3.1. Semi-continuous factorization

The factorization can be sought in this form where � is a SPD matrix of order n

Lh =(

− �

�xC. + �

)C−1

(C

�

�x. + �

)

= − �

�xC

�

�x− �

�x� + �

�

�x+ �C−1�

= − �

�xC

�

�x+ �C−1�.

It is thus necessary to have �C−1� = B. We have

� = C1/2A1/2C1/2, (8)

where A := C−1/2BC−1/2. Finally, we have the double equality

Lh =(

− �

�xC. + �

)C−1

(C

�

�x. + �

)=(

�

�xC. + �

)C−1

(−C

�

�x. + �

). (9)


Taking

Q1 =(

C�

�x+ �

)and Q2 =

(−C

�

�x+ �

)

leads to a convergence in two steps of (5), see [33] or [30]. This result is optimal in terms of iterationcounts. But, the matrix � is a priori a full matrix of order n costly to compute and use. Instead, we willuse approximations of it in terms of sparse matrices denoted �ap. We lose convergence in two steps. Inorder to have the best convergence rate, we design in Sections 4 and 5 optimized sparse approximationsto the matrix � w.r.t. the DDM.

3.2. Spectra of the substructured problem

We substructure in terms of

(H1H2

)=

⎛⎜⎜⎝(

C�

�x+ �ap

)(u)(

−C�

�x+ �ap

)(u)

⎞⎟⎟⎠ .

We need to compute T(H1, H2, 0) for arbitrary vectors H1, H2 ∈ Rn. From (9), the solution v2 toproblem (6) has the general following form:

v2 = exp(−C−1�x)(�) + exp(C−1�x)(�)

for some �, � ∈ Rn. Since the solution has to be bounded as x goes to infinity, we have � ≡ 0. Theboundary condition on � yields (� + �ap)(�) = H2 so that

v2 = exp(−C−1�x)(� + �ap)−1(H2).

Thus, Q1(v2)|x=0 = (−�(� + �ap)−1 + �ap(� + �ap)

−1)(H2). A similar formula holds for Q2(v1)|x=0.The substructured problem (7) has thus following form:

(I − �T(., ., 0))

(H1H2

)= G, (10)

where T(., ., 0) has the following expression:

T(., ., 0) = −(

(� − �ap)(� + �ap)−1 0

0 (� − �ap)(� + �ap)−1

)(11)

and

G = �T(0, 0, f ).

We have a first result relating the spectra of the substructured problem to the convergence rate of theadditive Schwarz method. For any matrix M with real eigenvalues, let �eff(M) be the ratio of its largesteigenvalue to its smallest one. The convergence of a Krylov method applied to a diagonalizable matrixM depends on �eff(M), see [36]. We have


Lemma 1. We assume that �ap is a SPD matrix of order n. Let Sc(�ap) be the convergence rate of theSchwarz algorithm, i.e., Sc(�ap) = max{||\ ∈ Sp((� − �ap)(� + �ap)

−1)}. We have that

Sc(�ap) < 1.

Moreover, the matrix Sub(�ap) := I − �T(., ., 0) has real eigenvalues in (0, 2) symmetric w.r.t. one and

�eff(Sub(�ap)) = 1 + Sc(�ap)

1 − Sc(�ap).

Proof. Let (v, ) be an eigenvector, eigenvalue of (� − �ap)(� + �ap)−1. We first prove that is real

and belongs to (−1, 1). Indeed, w = (� + �ap)−1(v) satisfies

(� − �ap)(w) = (� + �ap)(w).

We take the scalar product with w and get

= (�w, w) − (�apw, w)

(�w, w) + (�apw, w).

Since � and �ap are SPD matrices, we have proved that is real and belongs to (−1, 1). As for the secondpart of the proof, we have that (v, v)T, 1 − and (v, −v)T, 1 + are eigenmodes of Sub(�ap). Let usnotice that a very similar result may be found in [14]. �

Minimizing the effective condition number is thus equivalent to minimizing the convergence rate ofthe Schwarz algorithm. We now give a partial optimality result:

Lemma 2. Let �ap be a SPD matrix. Then,

min�∈R

�eff(Sub(��ap)) = �eff(Sub(�opt�ap)) = �(�−1/2ap ��−1/2

ap )1/2,

where

�opt = (�min(�−1ap �)�max(�

−1ap �))1/2.

Proof. We have

Sc(��ap) = max�∈Sp((��ap)

−1�)

∣∣∣∣1 − �

1 + �

∣∣∣∣= max

(∣∣∣∣∣1 − �min((��ap)−1�)

1 + �min((��ap)−1�)

∣∣∣∣∣ ,∣∣∣∣∣1 − �max((��ap)

−1�)

1 + �max((��ap)−1�)

∣∣∣∣∣)

.

This expression is minimized by taking � = �opt as defined in Lemma 2. In that case, we get

Sc(�opt�ap) = 1 − �

1 + �,


where

� :=√

�min(�−1ap �)/�max(�−1

ap �) = �(�−1/2ap ��−1/2

ap )−1/2.

Thus, we have (recalling that minimizing the convergence rate of the Schwarz method is equivalent tominimizing the effective condition number of the substructured problem)

min�∈R

�eff(Sub(��ap)) = �eff(Sub(�opt�ap)) = 1/� = �(�−1/2ap ��−1/2

ap )1/2. �

The above lemma shows that finding optimized sparse approximations to � w.r.t. the DDM reducesto finding optimal sparse preconditioners w.r.t. the condition number. In Section 4, we consider diagonalapproximations. In Section 5, we consider an approximation with a sparsity equals to that of matrix A.

4. Robin interface conditions

Notation: Consider the largest (resp., smallest) eigenvalue denoted by �Max(M) (resp., �min(M)) forany matrix M .

We consider the case where �ap is a diagonal matrix.1 From Lemma 2, we see that it is sufficient tofind an optimal diagonal preconditioner to matrix �. We shall use Theorem 4 (van der Sluis) which statesthat the diagonal of a SPD matrix is a quasi-optimal diagonal preconditioner. More precisely, we provea condition number estimate for the following choice:

�q−optap := �opt0C

1/2diag(A)1/2C1/2, (12)

where �opt0 is defined in Lemma 2. We have

Theorem 3.

�eff(Sub(�q−optap )�m1/4. min

D∈D �(D−1AD−1)1/4,

where D= {positive def inite diagonal matrices} and m is the maximum number of nonzeros in anyrow of A.

As an example, for a standard finite volume discretization for a three dimensional problem m = 5 andm1/4 = 1.49.

The sequel of the section is devoted to the proof of the theorem. We first give a series of results oflinear algebra. The basis for the proof is

Theorem 4 (van der Sluis). If F is a SPD matrix, then

minD∈D �(D−1/2FD−1/2)��(diag(F )−1/2Fdiag(F )−1/2)�m. min

D∈D �(D−1/2FD−1/2),

where D= {positive def inite diagonal matrices} and m is the maximum number of nonzeros in anyrow of F .

1 The authors thank Olivier Dubois, McGill University for his kind contribution to this section.


Details on this theorem can be found in [17,41].

Lemma 5. Let L be a nonsingular matrix with positive real eigenvalues, then

�(L) = �(LTL)1/2 ��Max(L)

�min(L).

Proof. See [16] �

Lemma 6. Let E and F be SPD matrices. Then,

�(E−1/4F 1/2E−1/4)2 ��(E−1/2FE−1/2).

Proof. Let E and F be any SPD matrices. Let us define L := F 1/2E−1/2. We have by Lemma 5,

�(E−1/2FE−1/2)��Max(F

1/2E−1/2)2

�min(F 1/2E−1/2)2 .

The spectrum of F 1/2E−1/2 is the same as the spectrum of F 1/4E−1/2F 1/4 which is symmetric

�(E−1/2FE−1/2)��Max(E

−1/4F 1/2E−1/4)2

�min(E−1/4F 1/2E−1/4)2 = �(E−1/4F 1/2E−1/4)2.

Proof of Theorem 3 is now easy. Indeed, by applying successively Lemmas 2, 6, Theorem 4, we have

�(Sub(�q−optap )) = �((�

q−optap )−1�)1/2 ��(diag(A)−1/2Adiag(A)−1/2)−1/4

�m1/4 minD∈D �(D−1/2AD−1/2)1/4. �

5. Two parameter interface condition

In the previous section, the interface condition is a Robin interface condition which reads fordomain �1

C�

�x+ �optC

1/2DC1/2,

where D = diag(A)1/2, see (12). In this section, we want to design more efficient interface conditions byconsidering more general interface conditions than Robin interface conditions.

Inspired by Higdon’s trick for absorbing boundary conditions [19] (see also [12]), we first consider aninterface condition of the form

Q :=(

C�

�x+ �1C

1/2DC1/2)(

C�

�x+ �2C

1/2DC1/2)

for some positive parameters �1, �2 and D is an invertible matrix not necessarily equal to diag(A)1/2.

This product yields a second-order derivative w.r.t. x the normal direction

Q := C�

�x

(C

�

�x

)+ (�1 + �2)C

1/2DC1/2C�

�x+ �1�2C

1/2DCDC1/2.


By using the operator Lh, this second-order term can be replaced by CB so that condition Q is equivalentto

Q := CB + (�1 + �2)C1/2DC1/2C

�

�x+ �1�2C

1/2DCDC1/2.

In order to be able to apply results of Section 3, we still have to write this condition in the form

C�

�x+ �ap,2

for some operator �ap,2. Since interface conditions are equivalent up to the left composition with anyinvertible operator acting along the interface, we obtain an equivalent condition R by left multiplying Q

by the inverse of (�1 + �2)C1/2DC1/2:

R := C�

�x+ C1/2 D−1A + �1�2D

�1 + �2C1/2. (13)

In other words, we choose to approximate � by

�ap,�1,�2 := C1/2 D−1A + �1�2D

�1 + �2C1/2 (14)

with �1, �2 > 0 . Let us notice that

(1) If D = diag(A)1/2, D−1/2AD−1/2 is another approximation to A1/2 that is consistent with approxi-mating A1/2 by D.

(2) The matrix A may be seen as a discretization matrix of a second-order partial differential operatorin the tangential directions to the interface. It is thus related to the optimized of order two interfaceconditions [1,20].

As in Section 4, we have to find the best parameters �1, �2 in (14). If D and A1/2 commute, it is easy toshow (see [12]) that the convergence rate is

Sc(�ap,�1,�2) = max�∈Sp(D−1A1/2)

∣∣∣∣� − �1

� + �1

∣∣∣∣∣∣∣∣� − �2

� + �2

∣∣∣∣ .

The spectrum of the matrix D−1A1/2 is discrete. If it is replaced by the segment [�m, �M ], minimizingSc(�ap,�1,�2) w.r.t. to the parameters �1 and �2 reduces to the optimization solved by Wachspress forADI methods [45] and whose solution is given in

Theorem 7. Suppose matrices D and A1/2 commute. Let �m:=�min(D−1A1/2) and �M :=�max(D

−1A1/2).The choice

�1,opt�2,opt = �m �M , (15)

�1,opt + �2,opt =(

2√

�M�m (�m + �M))1/2

(16)


yields a bound on the condition number

�eff(Sub(�ap,�1,opt,�2,opt))�1√2

(√�M

�m

+√

�m

�M

)1/2

.

We have considered in Theorems 3 and 7 optimized interface conditions for a decomposition into twosubdomains. For an arbitrary domain decomposition, the optimization problem becomes very difficultbut is not mandatory. Indeed, in practice, for an arbitrary domain decomposition, the optimized interfaceconditions are designed using two-domain decomposition results. Numerical results for various equationshave shown the adequacy of this approach, see [21] for the convection–diffusion equation, [2] or [24] forthe Helmholtz equation and the end of Section 6.3 in our case. Moreover, in [31] a theoretical result inthe constant coefficient case establishes a formula linking the convergence rate for an arbitrary numberof subdomains to the one of the two-subdomain case.

6. Numerical results for the semi-discrete problem

In this section, we test various interface conditions and algorithms in the semi-continuous frameworkof the previous sections. More precisely, except for Section 6.3, we work in 2D on the infinite tube� = R × (0, 1) and consider the operator

L = − �

�xc(y)

�

�x+ �(y) − �

�y�(y)

�

�y(17)

with Dirichlet boundary condition at the bottom and a Neumann boundary condition at the top in order toshow that the method is not limited to a specific boundary condition. We use a finite volume discretizationof the operator in the y direction which yields a tridiagonal matrix B of order ny . It is then possible toform the matrices of the substructured problems (10) for various interface conditions and study theirspectra. We either plot the spectra or give in the tables thee ratio of the largest norm of the eigenvaluesof the substructured matrix over its smallest real part. When the eigenvalues are real, this is the effectivecondition number. It corresponds to the case when both subdomains have the same coefficients. When thesubdomains are different, it is possible to have nonreal eigenvalues. Then, the ratio of the largest norm ofthe eigenvalues over its smallest real part is a good indicator of the performance of a GMRES algorithm[36]. We also give iteration counts (#iter in the tables) corresponding to solving Eq. (10) by a GMRESalgorithm [37] with a random right-hand side G. The stopping criterion is a reduction of the residualby a factor 10−6. Although we do not consider a discretization in the x direction, the results are a goodindication of what would happen in the corresponding fully discrete computations. The fully discretizedequations are considered in Section 6.3.

For the optimized interface conditions denoted by opt0 and opt2 below, computing parameters � byLemma 2 or Theorem 7 demands the computation of the square root of a matrix. In order to avoid this extracost, the extreme eigenvalues of the matrix product D−1A1/2 are approximated by the square root of theeigenvalues of D−1AD−1. Thus, there is no need to compute the square root of matrix A. Comparisonswith the exact formula have shown little difference in terms of iteration counts. Typical variations are of2 or 3 iterations.


Table 1Results for highly heterogeneous problems

ny 10 20 40 80 160

(opt0) #iter 11 17 22 28 37|�|max/real(�)min 6.8 31.4 48.8 71.9 1.1e+02

(opt2) #iter 9 11 15 17 18|�|max/real(�)min 1.8 3.8 4.9 5.9 7.2

(noprec) #iter 10 22 61 136 320�max/�min 7.3e+02 1.1e+04 2.5e+04 5.3e+04 1.1e+05

(diagprec) #iter 7 17 27 42 64�max/�min 42.7 1.1e+03 2.4e+03 5.1e+03 1.1e+04

We now define more precisely the names written in the tables and corresponding to the various DDMswhich have been tested

• opt0: The interface condition is the one studied in Section 4 except for the simplification mentionedabove in the previous paragraph.

• opt2: The interface condition is given by formula (13) where D = diag(A)1/2 and �1, �2 are given byformulas (15) and (16) except for the simplification mentioned above in the previous paragraph.

• noprec: The conjugate gradient is applied to the substructured system �(u)=G which corresponds toa Schur-type method without preconditioner. In this case, G is the jump of the fluxes at the interfacewhen homogeneous Dirichlet conditions are imposed at the interface.

• diagprec: The above system is preconditioned by its diagonal.

Let us remark that matrices D and A1/2 commute when the coefficients are constant in each subdomain.Although they do not commute in the general case, the computation of the parameters �1, �2 in opt2 isbased on Theorem 7. Results of Tables 1 and 2 show that the condition numbers are lower than the onesgiven by the other approaches although the coefficients are very are highly discontinuous and anisotropic.Moreover, we have performed a numerical optimization of the coefficients parameters which show indeedthat our choice of the coefficients is not far from the optimal ones. For instance, in the case ny = 80 ofTable 1, a numerical optimization yields a condition number of 1.41 instead of 5.9 when using Theorem7. The improvement is noticeable but anyway both figures are lower.

6.1. Highly heterogeneous problems

In Table 1, the diffusion coefficients are highly heterogeneous: c(y) = �(y) = val([10y]), where [ ] isthe integer part function and val is the vector val=[adabababab] where a =1.e4, b=1.e0 and d =1.e2.We have � = 0.

6.2. Different subdomains

In the above case, by symmetry of the problem w.r.t. the interface, a Neumann–Neumann or FETIalgorithm would give convergence in one iteration. In this section, we compare the optimized interface


Table 2Results for highly heterogeneous problems and different subdomains

ny 10 20 40 80 160 320

(opt0) #subdom. solves 8 22 32 40 48 56|�|max/real(�)min 1.9 25.6 43.5 65.1 94.1 1.3e+2

(opt2) #subdom. solves 8 11 13 15 15 16|�|max/real(�)min 7.6 3.6 4.6 5.7 6.8 8.2

(Neumann– #subdom. solves 12 18 24 28 32 32Kappa) �max/�min 22.1 31.9 35.6 40.7 47.8 59.7(Neumann– #subdom. solves 10 18 24 24 24 28MatKappa) �max/�min 1.9 2.2e+2 3.0e+2 4.2e+2 6.2e+2 9.6e+2

conditions approach developed so far to these algorithms when the operators in domains �1 and �2 arenot the same, i.e., coefficients C jump across the interface.

We now define more precisely the names written in the tables and corresponding to the various DDMswhich have been tested: opt0, opt2, NeumannKappa and NeumannMatKappa.

opt0 and opt2: In both cases, matrices �ap,i , i = 1, 2 are built separately as in Section 6. Theseapproximations do not take into account the fact they are used in a domain decomposition in which nowoperators vary from one domain to the other. Numerical results show that for opt2 iteration counts arestill good.

NeumannKappa: This corresponds to a Neumann–Neumann algorithm. The interface problem issymmetric. It is thus better to use a conjugate gradient algorithm applied to the substructured problem(�1 +�2)(u)=G preconditioned by w1�

−11 w1 +w2�−1

2 w2 with wi =C1/C1 +C2 and �i the DtN mapof domain �i , i = 1, 2.

NeumannMatKappa: The same as above except that the weights in the preconditioner come from thediscretization matrix wi is the diagonal of the discretization matrix of the problem.

For these last two methods, one iteration consists in solving a Dirichlet and a Neumann boundaryvalue problem in each subdomain. In Table 2, we consider a highly heterogeneous case: �1,2 = 0, c1(y)=val1([10y]) and val1 is the vector.

val1 = [bdbababbdb], where a = 1.e4, b = 1.e0 and d = 1.e2, �1(y) = val2([10y]) and val2 is thevector.

val2 = [babadabbeb], where a = 1.e4, b = 1.e0, d = 1.e2 and e = 1.e3, c2(y) = val3([10y]) and val3is the vector.

val3 = [abagbbagab], where a = 1.e4, b = 1.e0 and g = 1.e2 and �2(y) = val4([10y]) and val4 is thevector.

val4 = [badabaaadb], where a = 1.e0, b = 1.e4 and d = 1.e2.

Iteration counts for opt0 are significantly higher than in Table 1. Whereas, the interface conditionsopt2 are quite insensitive to the fact that operators are not the same in the subdomains. From the theoryfor Neumann–Neumann or FETI method (see e.g., [23,27] or [22] and references herein) the conditionnumber can be made independent of the jumps of the coefficients across the interface if these coefficientsare smooth in each subdomains. This is not the case here where coefficients jump both along and acrossthe interface. Therefore, results of Table 2 do not contradict these theoretical results.


Run: startFile.cere

Age: -2 Ma

X: Length in m

Y: Depth in m

Lithology

Fault

Seal

Shale

Sand

Source

startFile.cere - SnapShot #1 - Lithology - -2 Ma - ( m , m )

0 400 800 1200 1600 2000 2400 2800 3200 3600 4000

0

400

800

1200

1600

2000

0 400 800 1200 1600 2000 2400 2800 3200 3600 4000

0

400

800

1200

1600

2000

Fig. 1. Lithology for a bounded domain simulation. The interface is located between the two geological blocks. The computationhas been obtained with a prototype code developed at IFP.

6.3. Fully discretized equations

Previous results correspond to a semi-continuous problem on an infinite tube decomposed into twohalf tubes. A complete analysis of the fully discrete case is made in [10]. This analysis is necessary whenthe mesh size in the x direction is large. The analysis is similar to the one made here except that thefactorization (9) has to be replaced by a LDU factorization of the matrix discretizing the operator L. Theformula for the matrix � is more complex and hence its approximation by a sparse matrix is more complexas well. Here, we give results for a finite volume simulation performed on a domain bounded in both x

and y directions. We have only tested Robin and optimized of order 2 (opt2) interface conditions. For atwo-domain decomposition, the global computational domain is the rectangle [0, 8000] × [0, 2000] with160×40 discretization points. In Fig. 1, we plot the corresponding eigenvalues of the substructured systemfor the opt2 interface conditions for the computational domain and for a smaller one ([0, 4000]×[0, 2000]with 80 × 40 discretization points). The domain is composed of multiple layers with two lithologies:


2

1.8

1.6

1.4

1.2

1

0.6

0.8

0.4

0.2

00 10 20 30 40 50 60 70 80

Eigenvalues of the substructured problem

vp

Fig. 2. Eigenvalues of the substructured problem for various interface conditions: triangle: opt2 (smaller domain) , square: opt2.

� = 0.00788918 and � = 3.15567, see Fig. 2. The iteration counts are: 11 iterations for the opt2 interfaceconditions and 21 for Robin interface conditions for the two-domain decomposition. For the same kindof lithologies, we have tested decompositions in more than two subdomains. The iteration counts are: forthree subdomains: 15 iterations for opt2 and 47 for Robin and for five subdomains: 24 for opt2 and 43for Robin.

7. Conclusion

For problems with highly anisotropic and discontinuous coefficients, plateaus in the convergence ofKrylov methods exist even when using “good” preconditioners. A classical remedy is to use deflatedKrylov methods. We have developed in this paper a new algebraic approach in the DDM framework.We propose a way to compute optimized interface conditions for domain decomposition methods forsymmetric positive definite equations. Compared to deflation, only two extreme eigenvalues have to becomputed. Numerical results show that the approach is efficient and robust even with highly discontinuouscoefficients both across and inside subdomains. The nonsymmetric case is under study, see [15]. Theoptimization of the interface condition is then much more difficult. This was already the case for acontinuous approach of the problem, see [6]. This is even more the case at the discrete level where wehave no explicit formula for the discrete Steklov–Poincaré operator. Let us mention that such interfaceconditions can be used on nonmatching grids, see [9].


References

[1] Y. Achdou, C. Japhet, P. Le Tallec, F. Nataf, F. Rogier, M. Vidrascu, Domain decomposition methods for non-symmetricproblems, in: C.-H. Lai, P.E. BjZrstad, M. Cross, O.B. Widlund (Eds.), Eleventh International Conference on DomainDecomposition Methods, Domain Decomposition Press, Bergen, 1999, pp. 3–17.

[2] J.-D. Benamou, B. Deprés, A domain decomposition method for the helmholtz equation and related optimal controlproblems, J. Comput. Phys. 136 (1997) 68–82.

[3] X.C. Cai, M.A. Casarin, F.W. Elliott Jr., O.B. Widlund, Overlapping Schwarz algorithms for solving Helmholtz’s equation,in: J. Mandel, C. Farhat, X.C. Cai (Eds.), Domain Decomposition Methods, vol. 10 (Boulder, CO, 1997), AmericanMathematical Society, Providence, RI, 1998, pp. 391–399.

[4] B. Després, P. Joly, J.E. Roberts, A domain decomposition method for the harmonic Maxwell equations, in: IterativeMethods in Linear Algebra (Brussels 1991), North-Holland, Amsterdam, 1992, pp. 475–484.

[5] V. Dolean, S. Lanteri, F. Nataf, Optimized interface conditions for domain decomposition methods in fluid dynamics,Internat. J. Numer. Methods Fluids 40 (2002) 1539–1550.

[6] O. Dubois, Optimized schwarz methods for the advection–diffusion equation, Master’s Thesis, Mcgill University, 2003.[7] B. Engquist, H.-K. Zhao, Absorbing boundary conditions for domain decomposition, Appl. Numer. Math. 27 (4) (1998)

341–365.[8] I. Faille, E. Flauraud, F. Nataf, F. Schneider, F. Willien, Optimized interface conditions for sedimentary basin modeling, in:

I. Herrera, D.E. Keyes, O.B. Widlund, R.Yates (Eds.), 13th International Conference on Domain Decomposition Methods,Lyon, 2000.

[9] I. Faille, F. Nataf, L. Saas, F. Willien, Finite volume methods on non-matching grids with arbitrary interface conditionsand highly heterogeneous media, in: R. Kornhuber, R. Hoppe, J. Priaux, O. Pironneau, O. Widlund, J. Xu (Eds.), DomainDecomposition Methods in Science and Engineering Series, Lecture Notes in Computational Science and Engineering,vol. 40, 2004.

[10] E. Flauraud, F. Nataf, Optimized interface conditions in domain decomposition methods, Application at the semi-discreteand at the algebraic level to problems with extreme contrasts in the coefficients, Technical Report RI 524, CMAP, EcolePolytechnique, 2004.

[11] M.J. Gander, G.H. Golub, A non-overlapping optimized Schwarz method which converges with an arbitrarily weakdependence on h, in: Fourteenth International Conference on Domain Decomposition Methods, 2002.

[12] M.J. Gander, F. Magoulès, F. Nataf, Optimized Schwarz methods without overlap for the Helmholtz equation, SIAM J.Sci. Comput. 24(1) (2002) 38–60 (electronic).

[13] M.J. Gander, C. Rohde, Overlapping Schwarz waveform relaxation for convection dominated nonlinear conservation laws,Preprint 13-2003, Mathematisches Institut, Albert-Ludwigs-Universität Freiburg, 2003.

[14] M. Genseberger, Domain decomposition in the Jacobi–Davidson method for eigenproblems, Ph.D. Thesis, UtrechtUniversity, September 2001.

[15] L.G. Giorda, F. Nataf, Optimized schwarz methods for unsymmetric layered problems with strongly discontinuous andanisotropic coefficients, Technical Report 561, CMAP, CNRS UMR 7641, Ecole Polytechnique, France, 2004, submittedfor publication.

[16] G. Golub, C.F. Van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore, 1996.[17] A. Greenbaum, Iterative methods for solving linear systems, Frontiers inApplied Mathematics, vol. 17, Society for Industrial

and Applied Mathematics (SIAM), Philadelphia, PA, 1997.[18] T. Hagstrom, R.P. Tewarson, A. Jazcilevich, Numerical experiments on a domain decomposition algorithm for nonlinear

elliptic boundary value problems, Appl. Math. Lett. 1 (3) (1988).[19] R.L. Higdon, Absorbing boundary conditions for difference approximations to the multi-dimensional wave equations,

Math. Comput. 47 (176) (1986) 437–459.[20] C. Japhet, Conditions aux limites artificielles et décomposition de domaine: Méthode oo2 (optimisé d’ordre 2), application

à la résolution de problèmes en mécanique des fluides, Technical Report 373, CMAP (Ecole Polytechnique), 1997.[21] C. Japhet, F. Nataf, F. Rogier, The optimized order 2 method, application to convection-diffusion problems, Future

Generation Comput. Systems FUTURE 18 (2001).[22] A. Klawonn, O.B. Widlund, M. Dryja, Dual-Primal FETI methods for three-dimensional elliptic problems with

heterogeneous coefficients, SIAM J. Numer. Anal. 40 (2002) 159–179.


[23] P. Le Tallec, Domain decomposition methods in computational mechanics, in: J. Tinsley Oden (Ed.), ComputationalMechanics Advances, vol. 1(2), North-Holland, Amsterdam, 1994, pp. 121–220.

[24] S.C. Lee, M.N.Vouvakis, J.F. Lee,A non-overlapping domain decomposition method with non-matching grids for modelinglarge finite antenna arrays, J. Comput. Phys., 2005.

[25] P.L. Lions, On the Schwarz alternating method, III: a variant for nonoverlapping subdomains, in: T.F. Chan, R. Glowinski,J. Périaux, O. Widlund (Eds.), Third International Symposium on Domain Decomposition Methods for Partial DifferentialEquations, held in Houston, Texas, March 20–22, 1989, SIAM, Philadelphia, PA, 1990.

[26] G. Lube, L. Mueller, H. Mueller, A new non-overlapping domain decomposition method for stabilized finite elementmethods applied to the nonstationary Navier–Stokes equations, Numer. Linear Algebra Appl. 7 (2000) 449–472.

[27] J. Mandel, M. Bresina, Balancing domain decomposition for problems with large jumps in coefficients, Math. Comput. 65(216) (1996) 1387–1401.

[28] R.B. Morgan, GMRES with deflated restarting, SIAM J. Sci. Comput. 24 (2002).[29] R. Nabben, C. Vuik, A comparison of deflation and coarse grid correction applied to porous media flow, Technical Report

R03-10, Delft University of Technology, 2003.[30] F. Nataf, Interface connections in domain decomposition methods, in: Modern Methods in Scientific Computing and

Applications, NATO Advanced Study Institute, Université de Montréal, NATO Science Series II, vol. 75, 2001.[31] F. Nataf, F. Nier, Convergence rate of some domain decomposition methods for overlapping and nonoverlapping

subdomains, Numer. Math. 75 (3) (1997) 357–377.[32] F. Nataf, F. Rogier, Factorization of the convection–diffusion operator and the Schwarz algorithm, M3AS 5 (1) (1995)

67–93.[33] F. Nataf, F. Rogier, E. de Sturler, Optimal interface conditions for domain decomposition methods, Technical Report 301,

CMAP (Ecole Polytechnique), 1994.[34] R.A. Nicolaides, Deflation of conjugate gradients with applications to boundary value problems, SIAM J. Numer. Anal.

24 (2) (1987) 355–365.[35] F.X. Roux, F. Magoulès, S. Salmon, L. Series, Optimization of interface operator based on algebraic approach, in:

I. Herrera, D.E. Keyes, O.B. Widlund, R.Yates (Eds.), 14th International Conference on Domain Decomposition Methods,Cocoyoc, Mexico, 2002.

[36] Y. Saad, Iterative Methods for Sparse Linear Systems, PWS Publishing Company, 1996.[37] Y. Saad, M.H. Schultz, GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM

J. Sci. Statist. Comput. 7 (1986) 856–869.[38] Y. Saad, M. Yeung, J. Erhel, F. Guyomarc’h, A deflated version of the conjugate gradient algorithm, SIAM J. Sci. Comput.

21(5) (2000) 1909–1926 (electronic); Iterative Methods for Solving Systems of Algebraic Equations, Copper Mountain,CO, 1998.

[39] H. Kian, Tan, M.J.A. Borsboom, On generalized Schwarz coupling applied to advection-dominated problems, in: D.E.Keyes, J. Xu (Eds.), Seventh International Conference of Domain Decomposition Methods in Scientific and EngineeringComputing, held at Penn State University, October 27–30, 1993, American Mathematic Society, 1994, pp. 125–130.

[40] W.P. Tang, Generalized Schwarz splittings, SIAM J. Sci. Statist. Comput. 13 (2) (1992) 573–595.[41] A. van der Sluis, Condition numbers and equilibration matrices, Numer. Math. 14 (1969) 14–23.[42] A. van der Sluis, H. van der Vorst, The rate of convergence of conjugate gradients, Numer. Math. 48 (1986) 543–560.[43] S.Vandewalle, M. J. Gander, Optimized overlapping schwarz methods for parabolic pdes with time-delay, in: R. Kornhuber,

R. Hoppe, J. Périaux, O. Pironneau, O. Widlund, J. Xu (Eds.), Domain Decomposition Methods in Science and EngineeringSeries, Lecture Notes in Computational Science and Engineering, vol. 40, 2004.

[44] C. Vuik, A. Segal, J.A. Meijerink, An efficient preconditioned cg method for the solution of a class of layered problemswith extreme contrasts in the coefficients, J. Comput. Phys. 152 (1999) 385–403.

[45] L.W. Eugene, Optimum alternating-direction-implicit iteration parameters for a model problem, J. Soc. Ind. Appl. Math.10 (1962) 339–350.

[46] F. Willien, I. Faille, F. Nataf, Frédéric Schneider, Domain decomposition methods for fluid flow in porous medium, in:Sixth European Conference on the Mathematics of Oil Recovery, September 1998.

Documents

Optimized interface conditions in domain decomposition methods for problems with extreme contrasts in the coefficients