Optimized algebraic interface conditions in domain ... - Semantic Scholar

Optimized algebraic interface conditions in domain decomposition methods for strongly heterogeneous unsymmetric problems Luca Gerardo-Giorda1 and Frédéric Nataf2 1

2

Dipartimento di Matematica, Universit` a di Trento - Italy. (This author’s work was supported by the HPMI-GH-99-00012-05 Marie Curie Industry Fellowship at IFP - France) [email protected] ´ CNRS, UMR 7641, CMAP, Ecole Polytechnique - France [email protected]

1 Introduction Let Ω = R × Q, where Q is a bounded domain of R2 , and consider the elliptic PDE of advection-diffusion-reaction type given by −div (c∇u) + div (bu) + ηu = f in Ω Bu = g on R × ∂Q,

(1)

with the additional requirement that the solutions be bounded at infinity. After a finite element, finite differences or finite volume discretization, we obtain a large sparse system of linear equations, given by A w = f.

(2)

Under classical assumptions on the coefficients of the problem (e.g. η − 1 2 div b > 0 a.e. in Ω) the matrix A in (2) is definite positive. We solve problem (2) by means of an Optimized Schwarz Method: such methods have been introduced at the continuous level in [4], and at the discrete level in [5]. We design optimized interface conditions directly at the algebraic level, in order to guarantee robustness with respect to heterogeneities in the coefficients.

2 LDU factorization and absorbing boundary conditions In this section we illuminate the link between an LDU factorization of a matrix and the construction of absorbing conditions on the boundary of a domain (see [1]). As it is well known in domain decomposition literature, such conditions

2

Luca Gerardo-Giorda and Frédéric Nataf

e ∈ R3 be a can provide exact interface transmission operators. Let then Ω bounded polyedral domain. We assume that the underlying grid is obtained as a deformation of a Cartesian grid on the unit cube, so that for suitable integers Nx , Ny , and Nz , w ∈ RNx ×Ny ×Nz . If the unknowns are numbered lexicographically, the vector w is a collection of Nx sub-vectors wi ∈ RNy ×Nz , i.e. T T w = (w1T , . . . , wN ) . (3) x e reads From (3), the discrete problem in Ω B w = g,

(4)

where g = (g1 , .., gNx )T , each gi being a Ny ×Nz vector, and where the matrix B of the discrete problem has a block tri-diagonal structure   D1 U1    L1 D2 . . .  , B= (5)   .. ..  . . UNx −1  LNx −1 DNx where each block is a matrix of order Ny × Nz . An exact block factorization of the matrix B defined in (5) is given by B = (L + T)T−1 (U + T), where





0

  L1 . . . L=  ..  .

..

. LNx −1 0

    

(6)



 0 U1  .. ..    . . , U=   ..  . UNx −1  0

while T is a block-diagonal matrix whose nonzero entries are the blocks Ti defined recursively as  for i = 1  D1 Ti =  −1 Di − Li−1 Ti−1 Ui−1 for 1 < i ≤ Nx . At this time, we can give here the algebraic counterpart of absorbing boundary conditions. Assume g = (0, .., 0, gp+1 , .., gNx ), and let Np = Nx − p + 1. To reduce the size of the problem, we look for a block matrix K ∈ (RNy ×Nz )Np , ˜= each entry of which is a Ny × Nz matrix, such that the solution of Kv = g (0, gp+1 , .., gNx )T satisfies vk = wk+p−1 for k = 1, ..Np . The rows 2 through Np in the matrix K coincide with the last Np − 1 rows of the original matrix B. To identify the first row, which corresponds to the absorbing boundary

Algebraic OSM for strongly heterogeneous unsymmetric problems

3

condition, take as a right hand side in (4) the vector g = (0, .., 0, gp+1 , .., gNx ), and, owing to (6), consider the first p rows of the factorized problem      −1    T1 T1 U1 w1 T1 0 −1        L1 T2 . T2 U2 T2    ..   ..      =  . .    .. .. .. .. ..   wp     . . . . . 0 Tp Up wp+1 Lp−1 Tp Tp−1 The first two are p × p square invertible block matrices, so we need to consider only the third one, a rectangular p × (p + 1) matrix: from the last row we get Tp wp + Up wp+1 = 0,

(7)

which, identifying v1 = wp and v2 = wp+1 , provides the first row in matrix K. Assume now that g = (g1 , .., gq−1 , 0, .., 0)T . A similar procedure can be developed to reduce the size of the problem, by starting the recurrence in the factorization (6) from DNx , as  −1 Li for 1 ≤ i < Nx  Di − Ui Ti+1 e Ti =  DNx for i = Nx , and we can easily obtain the equation for the last row in the reduced equation as Lq wq−1 + Teq wq = 0. (8)

3 Optimal interface conditions for an infinite layered domain In this section we go back to problem (1), where the domain Ω is infinite ¯ =Ω ¯1 ∪ Ω ¯2 , in the x direction. We consider a two domain decomposition Ω Ω1 ∩ Ω2 = ∅, where Ω1 = R− × Q,

Ω2 = R+ × Q,

and we denote with Γ = ∂Ω1 ∩ ∂Ω2 the common interface of the two subdomains. We assume that the viscosity coefficients are layered (i.e. they do not depend on the x variable), and consider a discretization on a uniform grid via a finite volume scheme with an upwind treatment of the advective flux. The resulting linear system is given by       A11 A1Γ 0 w1 f1  AΓ 1 AΓ Γ AΓ 2   wΓ  =  fΓ  (9) 0 A2Γ A22 w2 f2

4


where wi is the vector of the internal unknowns in domain Ωi (i = 1, 2), and wΓ is the vector of interface unknowns. In order to guarantee the conservativity of the finite volume scheme, the vector of interface unknown consists of two sets of variables, wΓ = (wΓ , wλ )T , the first one expressing the continuity of the diffusive flux, the second expressing the continuity of the advective one. If the unknowns are numbered lexicographically, the matrix A is given by   .. .. .. .. . . . .   0   L D U 0 1 1 1     L1 D1Γ U1Γ    A =  · · · · · · 0 L1Γ DΓ Γ U2Γ 0 · · · · · ·  (10) ,   L D U 2Γ 2Γ 2     0 L2 D2 U2 0   .. .. .. .. . . . . where the block DΓ Γ is square, whereas the blocks LiΓ , and UiΓ (i = 1, 2) are rectangular. By duplicating the interface variables wΓ into wΓ,1 and wΓ,2 , we can define a Schwarz algorithm directly at the algebraic level, as k+1 v1 f1 A11 A1Γ = k+1 k − AΓ 2 v2k f Γ + (T1 − DΓ Γ ) vΓ,2 A Γ 1 T1 vΓ,1 (11) k+1 v2 f2 A22 A2Γ . = k+1 k − AΓ 1 v1k f Γ + (T2 − DΓ Γ ) vΓ,1 A Γ 2 T2 vΓ,2 As it is well known in literature, if we take T1 = AΓ Γ − AΓ 2 A−1 22 A2Γ

T2 = AΓ Γ − AΓ 1 A−1 11 A1Γ ,

the algorithm (11) converges in two iterations. We are in the position to give the following result, the proof of which will be given in [2]. Lemma 1. Let A be the matrix defined in (9), and let T1,∞ and T2∞ be such −1 −1 that T1,∞ = D1 − L1 T1,∞ U1 and T2,∞ = D2 − U2 T2,∞ L2 . We have −1 AΓ 1 A−1 11 A1Γ = L1Γ D1Γ − L1 T1,∞ U1

−1

U1Γ

−1

L2Γ .

−1 AΓ 2 A−1 22 A2Γ = U2Γ D2Γ − U2 T2,∞ L2

Noticing that AΓ Γ = DΓ Γ , the optimal interface operators are given by −1 −1 Tex U1Γ 1 = DΓ Γ − L1Γ D1Γ − L1 T1,∞ U1 (12) −1 −1 ex T2 = DΓ Γ − U2Γ D2Γ − U2 T2,∞ L2 L2Γ .


5

4 Optimized algebraic interface conditions for a non-overlapping Schwarz method ex The lack of sparsity of the matrices Tex 1 and T2 in (12), make them unsuitable in practice. Therefore we choose for T1 and T2 in (11) two suitable ex approximations of Tex 1 and T2 , respectively. At the cost of enlarging the size of the interface problem, we choose Tapp and 1 Tapp defined as follows: 2

−1 app −1 Tapp = DΓ Γ − L1Γ D1Γ − L1 (T1,∞ ) U1 U 1 −1 1Γ app −1 app L2Γ , T2 = DΓ Γ − U2Γ D2Γ − U2 (T2,∞ ) L2

(13)

app app where T1,∞ and T2,∞ are suitable sparse approximations of T1,∞ and T2,∞ , respectively. The most natural choice would be to take their diagonals, but, in order to have a usable condition, we wish to avoid the computation of both T1,∞ and T2,∞ , which is too costly. Notice that if Dj , Lj , and Uj (j = 1, 2) were all diagonal matrices the same would hold also for Tj,∞ . Moreover, if all the matrices involved commute, or if Lj = UjT , we would have

T1,∞

D1 + = 2

r

(−L1 )1/2 D1 (−U1 )−1/2 (−L1 )−1/2 D1 (−U1 )1/2 − L1 U1 . 4

and a similar formula holds for T2,∞ , with the roles of L2 and U2 exchanged. These considerations have led us to consider the following approximations of T1,∞ and T2,∞ . Let dj , lj , and uj be the diagonals of Dj , Lj and Uj , respectively. Robin: We choose in (13) D1 + α1opt D1 , 2

app T1,∞ =

√ where D1 = diag

d21 −4l1 u1 2

. The optimized parameter is given by

(α1opt )2 = max

q

r12 + I12 ,

q

r1 R1 − I12 ,

(14)

where we have set r1 := min Re λ, R1 := max Re λ, and I1 := max Im λ, √ −2 d21 −4l1 u1 (−L1 )1/2 D1 (−U1 )−1/2 (−L1 )−1/2 D1 (−U1 )1/2 λ ∈ σ − L1 U1 diag , 4 2 app with a similar formula for T2,∞ .

Order 2: This condition is obtained by blending together two first order approximations, and we have −1 app e1 , L1 ] + (α1 + α2 )L1 e12 + (α1 + α2 )D e1 + α1 α2 Id − L1 U1 T1,∞ = L1 [D D

6


e1 = where [., .] is the Lie bracket, where D and where

D1−1 D1 , 2

L1 = D1−1 L1 , U1 = D1−1 U1 ,

q p (α1 + α2 ) = 2 (r1 + R1 ) r1 R1 ,

2

2

(α1 α2 ) = r1 R1

(15)

r1 and R1 being defined as before. The tuning of the optimized parameters for both conditions can be found in [3], and a more exhaustive presentation of the construction of interface conditions and of the numerical tests will be given in a forthcoming paper [2]. The proposed interface conditions are built directly at the algebraic level, and are easy to implement. However, they rely heavily on the approximation of the Schur complement and, if on one hand the extension to a decomposition into strips appears quite straightforward, on the other hand further work needs to be done in order to analyse their scalability for an arbitrary decomposition of the computational domain. Finally, it is easy to prove the following result (see [2]). Lemma 2. The Schwarz algorithm k+1 v1 f1 A11 A1Γ = app app f Γ + (T2 − DΓ Γ ) vkΓ,2 − AΓ 2 vk2 AΓ 1 T2 vk+1 Γ,1

A22 A2Γ AΓ 2 Tapp 1

v2k+1 vk+1 Γ,2

=

f2 k k f Γ + (Tapp − D Γ Γ ) vΓ,1 − AΓ 1 v1 1

.

converges to the solution to problem (9). 4.1 Substructuring The iterative method can be substructured in order to use a Krylov type method and speed up the convergence. We introduce the auxiliary variables h1 = (Tapp − DΓ Γ ) vΓ,2 − AΓ 2 v2 , 2

h2 = −AΓ 1 v1 + (Tapp − DΓ Γ ) vΓ,1 , 1

and we define the interface operator Th     h1 −AΓ 1 v1 + (Tapp − DΓ Γ ) vΓ,1 1  Th :  h2  7−→  app f (T2 − DΓ Γ ) vΓ,2 − AΓ 2 v2 where f = (f1 , fΓ , f2 )T , whereas (v1 , vΓ,1 ) and (v2 , vΓ,2 ) are the solutions of A11 A1Γ v1 f1 = AΓ 1 Tapp vΓ,1 f Γ + h1 2 and


A22 A2Γ AΓ 2 Tapp 1

v2 vΓ,2

=

f2 f Γ + h2

7

.

So far, the substructuring operator is obtained simply by matching the conditions on the interface, and reads in matrix form T Id − ΠTh (h1 , h2 ) = F, (16) where Π is the swap operator on the interface, where F = ΠTh (0, 0, f ), and where the matrix Th is given in the following lemma (for a proof see [2]). Lemma 3. The matrix Th in (16) is given by   app app ex (T1 − Tex − DΓ Γ )−1 0 1 ) (T1 + T2 .  app ex ex −1 0 (Tapp − T ) (T + T − D ) ΓΓ 2 2 2 1

5 Numerical Results We consider problem (1) in Ω = R × (0, 1), with Dirichlet boundary conditions at the bottom and a Neumann boundary condition on the top. We use a finite volume discretization with an upwind scheme for the advective term. We build the matrices of the substructured problem for various interface conditions and we study their spectra. We give in the tables the iteration counts corresponding to the solution of the substructured problem by a GMRES algorithm with a random right hand side G, and the ratio of the largest modulus of the eigenvalues over the smallest real part. The stopping criterion for the GMRES algorithm is a reduction of the residual by a factor 10−10 . We consider both advection dominated and diffusion dominated flows, and different kind of heterogeneities. We report here the results for three different test cases. Test 1: the flow is advection dominated, the viscosity coefficients are layered, and the subdomains are symmetric with respect to the interface. Test 2: the flow is diffusion dominated, the viscosity coefficients are layered, but are not symmetric with respect to the interface. Test 3: the flow is diffusion dominated, the viscosity coefficients are layered, non symmetric w.r.t. the interface, and anisotropic, with an anisotropy ratio up to order 104 . The velocity field is diagonal with respect to the interface and constant. The R 6.1. A more detailed descripnumerical tests are performed with MATLAB tion of the test cases as well as futher numerical results can be found in a forthcoming paper [2]. Both conditions perform fairly well, in both terms of iteration counts and conditioning of the substructured problem, especially for the second order conditions, that show a good scalability with respect to the mesh size.

8

Luca Gerardo-Giorda and Frédéric Nataf p = q = 10 ny Test 1 iter Robin Order 2 cond Robin Order 2 Test 2 iter Robin Order 2 cond Robin Order 2 Test 3 iter Robin Order 2 cond Robin Order 2

10 4 4 1.05 1.01 7 6 1.61 1.21 9 7 5.42 1.54

20 6 5 1.25 1.02 10 6 1.83 1.26 17 10 18.27 2.75

40 8 6 1.68 1.14 13 8 2.59 1.30 27 14 24.75 4.48

80 11 8 3.27 1.34 16 11 3.52 1.83 35 16 31.04 5.92

160 16 9 6.57 1.61 19 15 3.94 2.76 42 19 38.32 6.32

320 23 10 13.51 1.92 21 19 4.12 3.68 47 21 47.29 6.86

Table 1. Iteration counts and condition number for the substructured problem in Tests 1-3

6 Conclusions We have proposed two kinds of algebraic interface conditions for unsymmetric elliptic problem, which appear to be very efficient and robust in term of iteration counts and conditioning of the problem with respect to the mesh size and the heterogeneities in the viscosity coefficients.

References 1. B. Engquist and A. Majda. Absorbing boundary conditions for the numerical simulation of waves. Math. Comp., 31(139):629–651, 1977. 2. Luca Gerardo Giorda and Frédéric Nataf. Optimized Algebraic Schwarz Methods for strongly heterogeneous and anisotropic layered problems. Technical Report 575, CMAP (Ecole Polytechnique), 2005. 3. Luca Gerardo Giorda and Frédéric Nataf. Optimized Schwarz Methods for unsymmetric layered problems with strongly discontinuous and anisotropic coefficients. J. Num. Math., 13(4):265–294, 2005. 4. Pierre-Louis Lions. On the Schwarz alternating method. III: a variant for nonoverlapping subdomains. In Tony F. Chan, Roland Glowinski, Jacques Périaux, and Olof Widlund, editors, Third International Symposium on Domain Decomposition Methods for Partial Differential Equations , held in Houston, Texas, March 20-22, 1989, Philadelphia, PA, 1990. SIAM. 5. F.-X. Roux, F. Magloulès, S. Salmon, and L. Series. Optimization of interface operator based on algebraic approach. In Ismael Herrera, David Keyes, Olof B. Widlund, and Robert Yates, editors, Domain Decomposition Methods in Sciences and Engineering, pages 297–304. UNAM, 2003. Proceedings from the Fourteenth International Conference, January 2002, Cocoyoc, Mexico.