Convergence Analysis of Domain Decomposition ... - CiteSeerX
INSTITUT NATIONAL DE RECHERCHE EN INFORMATIQUE ET AUTOMATIQUE
Convergence Analysis of Domain Decomposition algorithms with full overlapping for the advection-diffusion problems. Patrick Le Tallec , Moulay D. Tidriri
N˚ 2435 Octobre 1994
PROGRAMME 6 Calcul scientifique, modélisation et logiciel numérique
ISSN 0249-6399
apport de recherche
1994
Convergence Analysis of Domain Decomposition algorithms with full overlapping for the advection-di usion problems. Patrick Le Tallec , Moulay D. Tidriri Programme 6 | Calcul scienti que, modelisation et logiciel numerique Projet MENUSIN Rapport de recherche n 2435 | Octobre 1994 | 57 pages
Abstract: The aim of this paper is to study the convergence properties of a Time
Marching Algorithm solving Advection-Di usion problems on two domains using incompatible discretizations. The basic algorithm is rst presented, and theoretical or numerical results illustrate its convergence properties. This study is based on spectral theory, a priori estimates and a Di-Giorgi-Nash maximum principle . Key-words: non symmetric elliptic operators, domain decomposition, incompatible grids, advection-di usion, time marching algorithms, overlapping, maximun principle. (Resume : tsvp) Universite Paris-Dauphine and INRIA, Domaine de Voluceau- Rocquencourt- B.P. 105- Le
Chesnay Cedex (France) email: [email protected] Department of Mechanical Engineering, Yale University, New Haven, CT 06520 email: [email protected]. This work has been done in 1991 at INRIA and has been supported by the Hermes Research program under grant number RDAN 86.1/3. The author is presently supported by the National Science Foundation under contract number ECS-8957475 and by the United Technologies Research Center.
Unité de recherche INRIA Rocquencourt Domaine de Voluceau, Rocquencourt, BP 105, 78153 LE CHESNAY Cedex (France) Téléphone : (33 1) 39 63 55 11 – Télécopie : (33 1) 39 63 53 30
Analyse de Convergence des Methodes de Decomposition de Domaines avec Recouvrement total Resume : Le but de l'article propose est d'etudier les proprietes de convergence d'un algorithme de marche en temps applique a la solution d'un probleme d'advection di usion pose sur deux sous-domaines recouvrants. Les discretisations utilisees sur chaque sous-domaine peuvent ^etre totalement independantes. L'article rappelle les principes de l'algorithme utilise, en presente une analyse theorique, et une illustration numerique. L'etude theorique s'appuie sur des estimations apriori, sur le principe du maximun de Di Giorgi-Nash, et sur une analyse spectrale unidimensionnelle. Mots-cle : operateurs elliptiques nonsymetriques, decomposition de domaines, grilles incompatibles, advection-di usion, algorithmes de marche en temps, recouvrement, principe du maximum.
Convergence of Domain Decomposition Algorithms
3
1 Introduction Domain decomposition methods have recently appear as an ecient strategy for solving large scale problems on parallel computers ([1], [2], [3], [4], [5], [6]). Nevertheless, they can also be used in order to couple di erent models [11], [17]. In this paper we will examine a domain decomposition strategy which can be applied to such situations. This approach was introduced in order to solve several diculties that occur in
uid mechanics. In particular, our aim is to introduce several sub-domains in order to do one of the following :
Solve di erent problems on each subdomain. Use di erent kinds of approximation methods on each subdomain [7]. Use \local re nement techniques" or \mesh adaptive techniques", locally, per subdomain ([10]).
The subdomains fully overlap and the coupling is achieved through \friction" forces acting on the internal boundary of the domain, these friction forces being updated by an explicit time marching algorithm. The theoretical study of our method will be done on an Advection-Di usion problem, which will serve as our model problem from now on. The analysis will be made at the continuous level, independently of any discretization strategy, which means that the derived results will be mesh independent. In the next section we will describe this model problem. In the third section we will present our algorithm for some special cases. The fourth section will treat the one-dimensional stationary problem. We will show from this study that we can improve the convergence of this method by introducing a relaxation parameter (see [5]) and by replacing one of the Dirichlet condition by a Neumann type condition. The fth section will focus on the resulting multidimensional case. We rst study the proposed strategy for the time implicit problem. We then establish a priori estimates and show from these estimates exponential convergence of our algorithm. A detailed study of the semi-discrete algorithm will be done in section 6. In the last section we study the numerical stability of the explicit algorithm. Practical applications of the proposed algorithm to real life CFD problems can be found in [14], [18], [19], and [20].
RR n 2435
4
Patrick Le Tallec , Moulay D. Tidriri
Γe Ω
Γι Γb
ΩV
Figure 1: description of computational domain.
2 The model problem Consider a bounded domain, of IRn such that its boundary @ is lipschitzian and
loc = v a connected domain of IRn with loc ( g. 1). Let : ?b = @ \ @ loc ; ( internal boundary) ?i = @ loc \ ; ( interface) ?1 = @ n?b :(far eld boundary) We denote by n the external unit normal vector to @ or @ loc . Let v a velocity eld given such that :
8 < divv = 0 in ; : v:n = 0 on ?b:
(1)
The model problem we want to solve is the following:
Inria
Convergence of Domain Decomposition Algorithms
5
Find ', a real valued function, de ned on and satisfying :
8 > div(v') ? ' = 0 in ; > < ' = '1 on ?1 ; > > :
(2)
' = 0 on ?b :
Above v is the ow velocity and is the di usion coecient. Most CFD algorithms will in fact consider the solution of this problem as the stationary solution of the evolution problem (3) described below : Find : (0; T ) ! IR such that,
8 @ > > @t + div(v ) ? > > < > > > > : (0)
= 0 in (0; T ); = 1 on ?1 (0; T ); = 0 on ?b (0; T );
(3)
= 0 in : The general CFD algorithm consists then in integrating (3) in time until reaching a stationary solution.
3 General Algorithm 3.1 Time-Continuous case
Let us introduce the local sub-domain loc (see g. 1) which has as external boundary ?i , and let us consider the trace loc of on the sub-domain loc , as an independent variable, to which we associate an arbitrary independent initial value ol 6= o j loc . We now replace the evolution problem (3) by the following evolution system : Find (resp. loc ) : ! IR (resp. loc ! IR) satisfying :
8 @ > + div(v ) ? = 0 in (0; T ); > > < @t = 1 on ?1 (0; T ); > > > : @ = @ loc on ?b (0; T ); @n
RR n 2435
@n
(4)
6
Patrick Le Tallec , Moulay D. Tidriri
8 @ > loc + div(v ) ? > loc loc = 0 in loc (0; T ); > < @t loc = 0 on ?b (0; T ); > > > : loc = on ?i (0; T );
(5)
(0) = o in ; loc(0) = ol in loc :
(6)
Remark 3.1 The global problem (4)-(6) has no no-slip boundary condition. This
suppresses the boundary layer which appears at low viscosity and facilitates the numerical solution of this problem. The boundary layers are modeled by the local problems (5) (6) which are only to be solved on a small domain loc , with a very ne discretisation if needed. The two problems are only coupled by their boundary conditions and not by volumic interpolation.
Remark 3.2 The \Neumann" version of this method consists in replacing in (5) the boundary condition on ?i by the corresponding Neumann condition :
@ loc = @ on ? : i @n @n
3.2 Time Discrete case
The general algorithm that we propose for the solution of our model problem (2) is as usual to integrate in time the evolution problem (4)-(5)-(6) until we reach a stationary solution. This integration in time is then achieved by the following uncoupled semi-explicit algorithm, where the operators are treated implicitely inside each subdomain and where the coupling boundary conditions are treated explicitely :
set 0loc = ol and 0 = 0; then, for n 0; nloc and n being known, solve successively 8 n+1 n loc ? loc + div(v n+1 ) ? n+1 = 0 in ; > > loc loc loc > t < nloc+1 = n on ?i ; > > : n+1 = 0 on ?b ;
(7)
loc
Inria
7
Convergence of Domain Decomposition Algorithms
8 n+1 n > ? + div(v n+1 ) ? n+1 = 0 in ; > > > < t n+1 = 1 on ?1 ; > > n+1 n+1 > > : @ = @ loc on ?b : @n
(8)
@n
Remark 3.3 We have a full decoupling between (7) and (8). They can be discretized and solved by two independent solution techniques.
Remark 3.4 The "Neumann" version of this method consists in replacing the Dirichlet boundary condition on ?i by the corresponding Neumann condition :
@ nloc+1 = @ n on ? : i @n @n
Remark 3.5 The fully implicit version of this method consists in replacing the condition :
nloc+1 = n on ?i
@ nloc+1 = @ n on ? ) by the condition : i @n @n nloc+1 = n+1 on ?i ; @ n+1 @ n+1 on ? ). (resp. loc = i @n @n (resp.
The two subproblems are then coupled at each time step.
Remark 3.6 If we replace in (8) by E de ned as follows :
E = n loc; and ?b by ?i , and if we set t = 1, we obtain a nonoverlapping version of our strategy, which is a standard Dirichlet-Neumann algorithm [15], [16].
Remark 3.7 The initial condition ol is not assumed to be equal to 0 on the local
subdomain loc because in most cases this condition is impossible to impose at the discrete level since the grid used on loc will be in general di erent from the grid used on . In addition, even if we assume ol = 0 , we will not have nloc = n on
loc unless we use the fully implicit algorithm on compatible grids.
RR n 2435
8
Patrick Le Tallec , Moulay D. Tidriri
4 Stationary one-dimensional case
For t = +1, the above algorithm can be written :
set 0loc = 0 and 0 = 0, then, for n 0; nloc and n being known, solve 8 > div(v nloc+1 ) ? nloc+1 = 0 in loc ; > < nloc+1 = n on ?i; > > : n+1 loc
(9)
= 0 on ?b ;
8 > div(v n+1 ) ? n+1 = 0 in ; > > < n+1 = 1 on ?1 ; > > n+1 n+1 > loc on ?b : : @ @n = @ @n
(10)
In one space dimension, if we take as domain the interval ]0; 1[ of IR decomposed into two sub-domains =]0; 1[ and loc =]h2; 1[ then we obtain
0 < h2 < 1;
8 (n)0 > v'2 ? '(2n)" > > < '(2n) (h2) > > > : '(2n) (1) 8 (n)0 > v'1 ? '(1n)" > > < '(1n) (0) > > 0 > : '1(n) (1)
(11)
= 0 on ]h2 ; 1[; = '(1n?1)(h2 );
(12)
= b; = 0 on ]0; 1[; (13)
= a; 0
= '2(n) (1): By introducing relaxation parameters, we can also introduce the following variants of the above algorithm :
Inria
9
Convergence of Domain Decomposition Algorithms
8 (n)0 > v'2 ? '(2n)" = 0 on ]h2 ; 1[; > > < '(2n)(1) = b; > > > : '(2n) (h2 ) = 2 '(1n?1)(h2) + (1 ? 2 )'(2n?1) (h2); 8 (n)0 > v'1 ? '(1n)" = 0 on ]0; 1[; > > < '(1n) (0) = a; > > 0 0 0 : '1(n) (1) = 1'2(n) (1) + (1 ? 1 )'1(n?1) (1):
(14)
(15)
We will now try to produce conditions under which algorithm (15)-(14) converges, and those for which this convergence is optimal. For this purpose, we write the interface solution under the form 0 '1(n) (1) = 0(1) + n; '(2n) (h2) = (h2) + n ;
(16) (17) where is the converged solution we are looking for. Using the analytical solutions of problems (14) (15), we obtain the following induction formula
with
n
n
!
= MIN
(n?1)
(n?1)
!
:
0 1 ?( v ) v h 1 ? 2 2 e (e 2 ? 1) B CC v ? v v MIN = B B @ 1v(1 ? 2)( v ) e( )v12(e( )h2 ? 1) + (1 ? ) CA 1 e( )(h2?1) ? 1 (e( )(h2 ?1) ? 1)
(18)
(19)
This iterative process converges if the spectral radius of matrix MIN is less than 1. A direct but tedious calculation then yields :
Lemma 4.1 The spectral radius of the transfer matrix of algorithm (15)-(14) is : p (MIN ) = max[ 12 jD D2 ? 4Rj] (20) with
RR n 2435
10
Patrick Le Tallec , Moulay D. Tidriri
D = 2 ? (1 + 2 ) + 12 e(?v=) (e(v=)h2 ? 1) e?v= ev(1h2 =) ? 1
(21)
R = (1 ? 1 )(1 ? 2 ):
(22)
From this calculation, we can get the following results :
i) When h2 goes to 1 (nonoverlapping), D goes to +1, and then, (MIN ) goes to +1. There is no-convergence at this limit. ii) The optimal convergence is obtained in the case where all the eigenvalues of the matrix MIN are zero, i.e., when : D = 0 and R = 0. The latter conditions imply in particular
1 = 1 or 2 = 1: If we choose, in addition, 1 = 2 , the condition D = 0 implies h2 = 0. In this case the sub-domain loc is equal to the whole domain, and the associate algorithm has no interest. iii) The convergence of the method depends symmetrically on both relaxation parameters. From ii), it is reasonable to take one of the i equal to 1 and call the other . By setting : (?v= ) (v= )h2 A = 1 ? e e(v=()(eh2 ?1) ??1 1) (23) we have then
(MIN ) = j1 ? Aj:
In this case, we get the following convergence results :
Theorem 4.1 1) There is optimal convergence (in 1 iteration) if (?v= ) (v= )h2 = opt = f1 ? e e(v=()(eh2 ?1) ??1 1) g?1 < 1:
(24)
(25)
2) The algorithm converges for all in ]0; A2 [.
Inria
11
Convergence of Domain Decomposition Algorithms
Corollary 4.1 1) The case without relaxation ( = 1) converges only if : 2 1; A i.e., by setting d = 1 ? h2 (overlapping length), only if : d v Log (1 + e2?v= ) (stability condition). 2) When v goes to zero, we must have d 21 . Remark 4.1 This theorem states that the application of algorithm (15)-(14) to the steady case converges only if the overlapping d is suciently large. In the same situation, we will see that the unsteady approach will converge to the same steady solution but with less restrictions on d. This motivates the introduction of the time marching algorithm of section 3. Moreover, this time marching technique is well adapted to nonlinear problems such as those encountered in uid mechanics.
5 Implicit time discretization
5.1 The general algorithm
This section deals with the convergence analysis of the proposed algorithm in multiple dimension when one uses the fully implicit version of our strategy (4)-(6) : Set 0loc = ol and 0 = 0; then, for n 0; nloc and n being known, solve
8 n+1 n ? + div(v n+1 ) ? n+1 = > > > > < t n+1 = > > n+1 > > : @ @n = 8 n+1 n loc ? loc + div(v n+1 ) ? n+1 > > loc loc > t < nloc+1 > > > : n+1 loc
RR n 2435
0 in ;
1 on ?1 ; n+1 loc on ?b ; @ @n = 0 in loc ; = n+1 on ?i ; = 0 on ?b :
(26)
(27)
12
Patrick Le Tallec , Moulay D. Tidriri
5.2 Convergence analysis
We have the following convergence result : Theorem 5.1 The solution of algorithm (26)-(27) converges linearly in H 1( ) towards the solution of the stationary problem (2), for all values of t and all choices of loc .
Proof Step 1: Local L2 estimates of ? loc
Substracting (27) from (26), multiplying the result by n+1 ? nloc+1 and integrating by parts, we obtain the following relation : Z 1 Z 1 2 n +1 n +1 ( ? loc ) ? ( n ? nloc )( n+1 ? nloc+1 )
loc t
locZ t (28) n +1 2 n +1 + jr( ? loc )j = 0:
loc
By using Cauchy-Schwarz, we deduce 1 k n+1 ? n+1 k2 n+1 ? n+1 j2 loc 1;2; loc loc 0;2; loc + j 2t (29) 1 n n 2 2t k ? lock0;2; loc : From the Poincare inequality, it follows that (30) k n+1 ? nloc+1k20;2; loc 1 + 21 tc k n ? nlock20;2; loc where c is the Poincare constant. By setting co = k o ? oloc k20;2; loc , we get by induction our basic L2 estimate k n+1 ? nloc+1 k20;2; loc ( 1 + 21 tc )n+1co: (31) By summing (29), we also obtain
k n+1 ? nloc+1k20;2; loc + 2 t
n X i=p
j i+1 ? iloc+1j21;2; loc
k p ? ploc k20;2; loc :
(32)
Inria
13
Convergence of Domain Decomposition Algorithms
In particular, we obtain for p = 0
k n+1 ? nloc+1k20;2; loc + 2 t
n X i=0
j i+1 ? iloc+1j21;2; loc
k o ? oloc k20;2; loc :
(33)
Step 2: Local H 1 estimates of n ? nloc Substracting the two rst equations in (26) and (27), multiplying the result by
xn = (
n+1 ? n+1 ) ? ( n ? n ) loc ; loc
and integrating on loc , we obtain 0 =
Z Z
loc
+
?
Z
jxnj2 +
loc
Z
t
(34)
div(v ( n+1 ? nloc+1 ))xn
(35)
loc ( n+1 n+1 )
r
?
loc
rxn
@ ( n+1 ? n+1 )xn: loc @ loc @n
The use of the third relation in (26) and the second relation in (27) then yields
Z
loc
jxnj2 +
Z
loc
div(v ( n+1 ? nloc+1 ))xn +
Z
loc
r( n+1 ? nloc+1 )rxn = 0:
By using now Cauchy-Schwarz, we obtain
Z
loc
jxnj2 +
Z
2 r( n+1 ? nloc+1)rxn kv2k1 j n+1 ? nloc+1 j21;2; loc
loc + 1 kxn k20;2; loc : 2
By construction of xn and Cauchy-Schwartz, the last inequality becomes
RR n 2435
14
Patrick Le Tallec , Moulay D. Tidriri
kvk21 j n+1 ? n+1j2 1 kxn k2 + j n ? nloc j21;2; loc 1 ; 2 ;
0 ; 2 ;
loc loc loc 2 2 2t
(36)
? 2 t j n+1 ? nloc+1j21;2; loc :
By using again relation (29), it follows that 2 kxnk20;2; loc k2vk1t (k n ? nlock20;2; loc ? k n+1 ? nloc+1 k20;2; loc ) + t (j n ? nloc j21;2; loc ? j n+1 ? nloc+1 j21;2; loc ) t (G(n) ? G(n + 1))
(37)
with G de ned by
2 G(n) = k2vk21 k n ? nloc k20;2; loc + j n ? nloc j21;2; loc : From (37), G is a decreasing function. This last property then yields
(n + 2 ? p)G(n + 1)
nX +1 i=p nX +1 i=p
(38)
G(i)
X
2 n+1
j i ? ilocj21;2; loc + k2vk21
i=p
k i ? ilock20;2; loc (39)
By using (29) again, we obtain
nX +1 +1 2 nX (n + 2 ? p)G(n + 1) 21t k i?1 ? iloc?1 k20;2; loc + k2vk21 k i ? iloc k20;2; loc i=p
i=p
But since k i ? iloc k0;2; loc is a decreasing function, we nally have (n + 2 ? p)G(n + 1) 2 (n + 2 ? p)( 21t + k2vk21 )k p?1 ? ploc?1 k20;2; loc ;
(40)
Inria
Convergence of Domain Decomposition Algorithms
15
that is 2 G(n + 1) ( 21t + k2vk21 )k p?1 ? ploc?1 k20;2; loc ; 8p n: Step 3: Global L2 estimate of n ? sta
(41)
We now de ne n+1 by
8 n+1 < ? sta in n loc; n+1 = : n+1 ? sta in loc; loc
(42)
with sta the solution of the stationary problem (2). By construction, n+1 satis es the following equations :
8 n+1 n ? + div(v n+1 ) ? n+1 = 0 in [ ( n ); > > loc loc > < t n+1 = 0 on @ ; > > > : n+1 continuous across ?i:
(43)
By multiplying this equation by n+1 and integrating by parts over loc and
n loc , we get Z n+1 ? n Z n+1 + jr n+1 j2
Z t1 @ n+1 ) n +1 2 (44) + ( 2 ( ) v:n ? n+1 @n Z@ loc 1 + ( ( n+1 )2(v:n) ? n+1 @ n+1 ) = 0: @n @ ( n loc ) 2 By taking into account the boundary conditions in (43) we obtain the following relation : Z n+1 ? n Z n+1 + jr n+1 j2
(45) Z @t n +1 n +1 n +1 ? @n ( loc ? ) = 0: ?i
RR n 2435
16
Patrick Le Tallec , Moulay D. Tidriri
But on loc ; nloc+1 ? n+1 satis es the equation ( nloc+1 ? n+1 ) ? ( nloc ? n ) + div[v ( nloc+1 ? n+1 )] ? ( nloc+1 ? n+1 ) = 0: (46) t Therefore, after multiplication by n+1 , integration by parts and from CauchySchwarz, we obtain
j
Z @ ( nloc+1 ? n+1 ) n+1 j 12 kxn k20;2; loc + 1 c21 k n+1 k20;2; loc @n 2c1 2 ?i 2 + 1 c2k n+1 k2 + kv k1 j n+1 ? n+1 j2 2c21
loc
1;2; loc
2
1
+ j nloc+1 ? n+1 j21;2; loc + j n+1 j21;2; loc ; 2 2 with c1 > 0 arbitrary. This latter inequality combined with (45) yields
0;2; loc
(47)
Z n+1 ? n n+1 + 2 j n+1 j21;2; ? c21k n+1 k20;2; loc t
21c2 kxnk20;2; loc + ( k2vck21 + 2 )j nloc+1 ? n+1 j21;2; loc : 2
1
1
By using Cauchy-Schwarz again we deduce : 1 n 2 1 n 2 1 ? c2 )k n+1 k2 + j n+1 j2 ( 2 0;2; 2 1;2; 2t k k0;2; + 2c2 kx k0;2; loc + t 1 1 2 k v k ( 2c21 + 2 )j nloc+1 ? n+1 j21;2; loc 1 (48) By using now Poincare inequality, we have (1 + (c ? 2c21)t)k n+1 k20;2; k n k20;2; + c2t kxn k20;2; loc 1
2 +t( kv k1 + )j n+1 ? n+1 j2
c21
loc
1;2; loc :
(49)
Inria
Convergence of Domain Decomposition Algorithms
17
This writes
k n+1k20;2; Ak nk20;2; + B1kxnk20;2; loc +B2 j nloc+1 ? n+1 j21;2; loc
with
(50)
A = 1 + (c ?1 2c2 )t ; B1 = c2t A;
1
1 v 21 c21
B2 = t( k k + )A:
Step 4: Final result
Let us choose c1 such that c ? 2c21 > 0. By induction, it follows ?1 Ai (B kxn?i k2 k n+1 k20;2; Apk n+1?p k20;2; + Ppi=0 1 0;2; loc +
B2j nloc+1?i ? n+1?i j21;2; loc ): Since A < 1 by assumption on c1, this implies
n X k n+1k20;2; Apk n+1?pk20;2; + A(B1 kxik20;2; loc + i=n+1?p n X j iloc+1 ? i+1j21;2; loc ): B2
(51)
(52)
i=n+1?p
Now, by using (37) and (32), we obtain
k n+1k20;2; Apk n+1?p k20;2; + A(B1 t (G(n + 1 ? p) ? G(n + 1))+ B2 21t k nloc+1?p ? n+1?p k20;2; loc ):
The same relation written between 0 and n + 1 ? p yields
RR n 2435
(53)
18
Patrick Le Tallec , Moulay D. Tidriri
k n+1?p k20;2; An+1?p k ok20;2; loc + A(B1 t (G(0) ? G(n ? p + 1)) +B2 21t k 0 ? 0loc k20;2; loc ):
By combining this relation with (53), we nally obtain k n+1 k20;2; An+1 k 0k20;2;
+ Ap+1 (B1 t G(0) + B2 21t k 0 ? 0loc k20;2; loc )
+ A(B1 t G(n + 1 ? p) + 2B2 t k nloc+1?p ? n+1?p k20;2; loc ):
By taking p such that n = 2p + q and by using (41) we conclude that
k n+1 k20;2; An+1C2 + Ap+1C3 + C4k ploc ? pk20;2; loc which implies from (31) the linear convergence of k n+1 k20;2; towards 0.
On the other hand by using (37) in (48) we obtain 1 ? c2 )k n+1 k2 + j n+1 j2 1 k n k2 ( 2 0;2; 2 1;2;
0;2;
t 1 2t 2 k v k + 2 (G(n) ? G(n + 1)) + ( 21 + 2 )j nloc+1 ? n+1 j21;2; loc : 2c1t 2c1 Therefore by using (29) we obtain 1 ? c2)k n+1 k2 + j n+1 j2 1 k n k2 ( 2 0;2; 2 1;2;
t 1 2t n 0;2; n 2 2 k ? k k v k + 2c2t (G(n) ? G(n + 1)) + ( 2c21 + 2 )( loc 2 t0;2; loc ): 1
(54)
(55)
(56)
1
Our result follows then from (41), (31) and the linear convergence of k n k0;2; .
6 Theoretical analysis of the explicit method
6.1 Local estimates
In this section, we will try to establish local estimates and a maximum principle for an arbitrary elliptic operator of second order. These tools are then needed in the convergence analysis of the explicit version of our coupling method.
Inria
19
Convergence of Domain Decomposition Algorithms
Let L be an operator written under the form
Lu = aij (x)Dij u + bi(x)Diu + c(x)u; for any u in W 2;n ( ), with a bounded domain of IRn . The coecients aij ; bi and c; i; j = 1; :::; n are de ned on . As usual, the repeated indices indicate a summation from 1 to n. We suppose that the operator L is strictly elliptic in in the sense that the matrix A of coecients [aij ] is strictly positive everywhere in . Let and denote respectively the smallest and the biggest eigenvalue of A. Let D denote the determinant of the matrix A and D = D1=n . We have 0 < D : We suppose in addition that the coecients aij ; bi and c are bounded in , and that there exists two positive real numbers and such that : = ;
(L is uniformly elliptic)
(57)
(jbj=)2 : (58) Now, we are in a position to state the principal result of this section, proved in annex. Theorem 6.1 Let u 2 W 2;n( ) and suppose that Lu f with f 2 Ln( ) and c 0. Then for all spheres B = B2R (y) of center y and radius 2R included in
and for all p > 0, we have : Z sup u ctef( 1 (u+)p)1=p + R kf k g; (59) B R (y )
jBj
B
n;B
where the constant cte depends on (n; ; R2; p), but is independent of c. Above u+ = max(u; 0).
Remark 6.1 The statement of the same theorem can be found in [12], under the assumption
jcj= :
(60) So, there the constant cte depends indirectly on c through . That is exactly what we want to avoid, since we need this constant to be independent of c, that is of t.
RR n 2435
20
Patrick Le Tallec , Moulay D. Tidriri
6.2 Technical results
In this section, we will try to apply the \maximum principle" of the previous section to study the model subproblems used in our algorithm (7)-(8).
6.2.1 Passing from the local to the global solution Let W be the sub-space of H 1( ) de ned by
W = fw 2 H 1( ); w = 0 on ?1 g: We then de ne the bilinear forms on W a(v; w) =
Z
rvrw +
(v; w) =
Z
Z
div(V v )w;
vw:
(61) (62) (63)
In the above equation, V is the velocity eld de ned by the relation (1), i.e., V satis es : 8 < divV = 0 in ; (64) : V:n = 0 on ?b : The rst basic problem associated to (8), can be written as follows : Find v 2 W satisfying :
a(v; w) + (1= )(v; w) =
Z
?b
gwd?; 8w 2 W:
(65)
Here the function g is given in H ?1=2(?b ) and the coecient is strictly positive. Moreover, from now on, we suppose that the coecients and satisfy the following relation : 1: (66) This hypothesis is not necessary but simpli es the proofs to come. Moreover, it is not restrictive, since we want the convergence for small corresponding to t small. Let d denotes the overlapping distance as in gure 2. Let then be a real number such that p 0 < < 3 =d;
Inria
21
Convergence of Domain Decomposition Algorithms
and set
p
k = =( ):
We then have the following result :
Theorem 6.2 Let v the solution of problem (65). If is suciently small, we have
q 1=2 1 C1 C2=d 1 + kV k1 d=C2 q
kvk1=2;?i
(1= )exp(?kd2=36)kg k?1=2;?b ;
(67)
where C1 and C2 are constants, with C1 depending only on d; ; and , but not on .
The proof of this theorem will be done in ve steps and is based on the local estimates described in Theorem 6.1 and on the maximum principle.
Step 1. Global estimates
By using relation (64), we have the equality :
Z
Z
vdiv(V v) = 1=2 div(V v 2)
Z
= 1=2 V:nv 2 ?
= 0; 8 v 2 W: By taking w = v in (65), we then obtain
Z
f jrvj
2 + (1= )v 2
g=
Z ?b
gv:
(68)
From this equality we deduce the following bound :
kvk21;2; kgk?1=2;?b kvk1=2;?b : Now by using the trace theorem, we obtain :
kvk1;2; (co= )kgk?1=2;?b; RR n 2435
(69)
22
Patrick Le Tallec , Moulay D. Tidriri
which implies in particular
kvk0;2; (co= )kgk?1=2;?b: Step 2. Nearby local estimates
(70)
Let L be the second order operator associate to problem (65) :
Lv = ? v + V:rv + (1= )v:
The operator
(71)
L0 = ?L
satis es the assumptions of Theorem 6.1, with c = ?1= ; f = 0. Applying then Theorem 6.1 for p = 2; y 2 i (the subdomain of width 3d externally delimited by ?i ) ( g. 2) and Ky = Bd=3 (y ), the sphere of center y and radius d=3 we obtain
kvk1;Ky c1kvk0;2;B2d=3(y): Therefore
kvk1;Ky c1kvk0;2; ;
where c1 is a constant depending only on ; ; d2 and (3=2d)n=2 . Moreover, there exist y1 ; : : :; yl belonging to i such that
(72)
2i = [y2 i B d6 (y ) [lj =1 Kyj = K with
Kyj = Bd=3 (yj ): By applying the relation (72) to each Kyj we obtain :
kvk1;K sup c1j kvk0;2; : j =1;:::;l
By setting c1 = sup c1j , we nally get j =1;:::;l
kvk1;K c1kvk0;2; : Step 3. Use of maximum principle
(73)
For any Mi in i ; we introduce (see g. (2)) :
Inria
Convergence of Domain Decomposition Algorithms
23
Bi = the ball centered on Mi of radius d=6, vi = exp[k(r2 ? d2=36)]kvk1;@Bi: The operator L applied to vi , can be written in polar coordinates (with r =
Mi M ) :
Therefore We set then : with
Lvi = 4(?k2 r2 ? k + k2 V:err + 41 )vi: Lvi 4(?k2r2 ? k2 jV:erjr + ( 41 ? k ))vi:
(74)
'(r; k) = a(k)r2 + b(k)r + c(k);
(75)
a(k) = ?k2 b(k) = ? k jV:e j
(76)
r
(77)
c(k) = 41 ? k:
(78)
2
We seek to satisfy the following relation :
0 inf '(r; k) for 0 r d6 :
As '(r; k) decreases on R+ , this will be satis ed if and only if '(d=6) 0; i.e. if and only if 2 2 kV k + 1 ? k 0: ? k 36d ? kd12 4
We replace k by its value. We therefore have to satisfy 2 d2 ? dkp V k 1 ? (36 ) 12 + 4 ? p 0: RR n 2435
24
Patrick Le Tallec , Moulay D. Tidriri
p
Multiplying by , it remains 1 (1 ? 2 d2 ) (1 + dkV k ): p 4 9 12 p From the constraint < 3 =d, we nally obtain after division
kV k ][1 ? 2d2 ]?1: '(r; k) 0 i 4p1 [1 + d12 9
(79)
p
From the relation (74) and the previous calculation, we deduce that for < 3 =d and satisfying the above inequality, we have In addition, by construction
Lvi 0 = Lv: vi v on @Bi :
Consequently, by using the maximum principle we end up at the following relation :
v vi on Bi :
In particular
v(Mi) exp(?kd2=36)kvk1;@Bi :
We do the same for ?v , and nally we have
jv(Mi)j exp(?kd2=36)kvk1;@Bi ; 8Mi 2 i: Step 4 : H 1 estimate Let 2 H 1( ) be such that : 8 < = 1 in 1 = ? loc ; : supp i [ 1 : We have from (65),
Z
(? v + div(V v ) + v= ) 2v = 0:
(80)
(81)
By using the Green's formula we deduce :
Inria
25
Convergence of Domain Decomposition Algorithms
Z
Z Z ? v 2 v = jr(v)j2 ? jrj2v 2:
On the other hand, we have :
Z
Z
div(V v ) 2v =
div(V 2 v 2=2) ?
Z
With the relations (82)-(83), (81) becomes
Z
0 =
Z
=
(83)
Z
Z
=
V:rv 2:
(82)
( jr(v )j2 + div(V 2v 2 =2) + 2 v 2= ) ? (v 2jr j2 + v 2V r )
Z
(jr(v)j2 + jvj2) + (1= ? ) 2 v 2 ? (v 2 jrj2 + v 2 V:r)
1
?
Z
(jrvj2 + jvj2) +
Z
i
Z
i
Z
(jr(v)j2 + jvj2) + (1= ? ) 2v2
(v 2jr j2 + v 2V:r ):
Hence, we obtain :
kvk
2 1;2; 1 +
Z
i
(jr
j
(v ) 2 +
jvj Z
i
2) +
Z
(1= ? ) 2v 2 =
(v 2jr j2 + v 2V r ):
With the relation (66), it follows
kvk21;2; 1[ i
Z
i
(v 2 jr j2 + v 2V:r )
kvk21; i
Z
i
( jr j2 + V:r )
(84) (85)
kvk21; i ( jj21;2; i + kk0;2; i kV kjj1;2; i ) kvk21; i jj21;2; i ( + kk0;2; i kV k=jj1;2; i ): RR n 2435
(86)
26
Patrick Le Tallec , Moulay D. Tidriri
If we take such that
jj0;2; i 1; jj21;2; i = c2=d;
where c2 is a constant, (86) then becomes
q q 1=2 kvk1;2; 1 kvk1; i c2=d 1 + kV k d=c2 : Step 5: Getting the nal estimate.
(87)
We have
kvk1;@Bi kvk1;K ;
(88)
kvk1;@Bi c1kvk0;2; :
(89)
since @Bi K . From (73), this implies
By using the relations (80) and (89) it follows :
jv(Mi)j exp(?kd2=36)c1kvk0;2; ; 8Mi 2 i:
Consequently we have
kvk1; i exp(?kd2=36)c1kvk0;2; :
(90)
kvk1; i c1co exp(?kd2=36)kgk?1=2;?b :
(91)
By using now the relation (70), we obtain :
Now, by using the relation (87) we obtain :
q 1=2 q kvk1;2; 1 coc1 c2=d 1 + kV k d=c2
(92)
(1= )exp(?kd2=36)kg k?1=2;?b:
Inria
27
Convergence of Domain Decomposition Algorithms
To conclude we use the trace theorem which yields
kvk1=2;?i c3kvk1;2; 1 : Consequently, we arrive at the following relation :
q q kvk1=2;?i coc1c3 c2=d(1 + kV k d=c2)1=2
(93)
(1= )exp(?kd2=36)kg k?1=2;?b : This is our theorem with C1 = coc1 c3 and C2 = c2.
6.2.2 Passing from the global to the local solution
Let W denote now the sub-space of H 1 ( loc ) de ned by W = fw 2 H 1( loc )=w = 0 on ?b g: We then set for any v and w in W
a(v; w) =
Z
loc
rv:rw + Z
(v; w) =
loc
Z
loc
div(V v )w;
vw;
(94) (95)
where V is the velocity eld de ned by the relation (64). The second basic problem associated to (7) can be written as : Find v 2 W such that
a(v; w) + (1= )(v; w) =
Z
?i
vj?i = h;
@v w; 8w 2 W; @n
where h is given in H 1=2(?i ). We rst have the following lemma ; Lemma 6.1 For suciently small, we have a(w; w) + (1= )(w; w) (=2)kwk21;2; loc ; 8w 2 W:
RR n 2435
(96) (97)
28
Patrick Le Tallec , Moulay D. Tidriri
Γi
Ωi
Γ1
y Ky
d/3
d/3
Γb
d/3
d/3
Ω 2ι
Mi Bi
δΒι
Figure 2: Description of the Domain loc and of the splitting used in the majoration of the local solution.
Inria
29
Convergence of Domain Decomposition Algorithms
Proof : Under the hypothesis 1= =2 + (1=2 )kV k21 , and using the Cauchy-Schwarz inequality, we obtain :
a(v; v) + (1= )(v; v) =
Z
loc
rv:rv +
Z
loc
V:rvv + (1= )
Z
loc
v2
krvk20;2 + (1= )kvk20;2 ? kV k1krvk0;2kvk0;2 krvk20;2 + (1= )kvk20;2 ? (=2)krvk20;2 ?(1=2)kV k21kvk20;2 (=2)(krvk20;2 + kvk20;2): We will also make the simplifying assumption (66). We have then the following result : Theorem 6.3 Under the notation of Theorem 6.2 and for suciently small, the solution of problem (96) is bounded by
k@v=@nk?1=2;?b
q
2 C1 C2=d 1 + 1 + kV2 k1
!
q 1=2 k V k 1 d=C 1+
2
(1 + 1= 2)exp(?kd2=36)khk1=2;?i : Here C1 and C2 are constants with C1 depending on d; v; and .
Proof
To prove this theorem, we will proceed as in the proof of Theorem 6.2.
Step 1 : Global estimates. RR n 2435
(98)
30
Patrick Le Tallec , Moulay D. Tidriri
By taking w = v in (96), we obtain :
Z
loc
jrvj2 +
Z
loc
(div(V v )v + (1= )v 2) =
Z ?i
Using the above Lemma and setting = =2, this implies
@v h: @n
(99)
(100) kvk21;2; loc k@v=@nk?1=2;?i khk1=2;?i : We will try now to get estimates on k@v=@nk?1=2;?i . From (96) and (64) we obtain : Z @v Z 1 vw): w = ( r v r w + (1 = ) V: r vw + ?i @n
loc Therefore, for any w in W , we have
j
Z @v @n wj krvk0;2; loc krwk0;2; loc + (1= )kV k1 krvk0;2; loc kwk0;2; loc ?i
+ 1 kv k0;2; loc kwk0;2; loc
(krvk20;2; loc + (1= 2)kV k21krvk20;2; loc + (1= 2)kvk20;2; loc)1=2 (krwk20;2; loc + kwk20;2; loc + (1= 2)kwk20;2; loc )1=2
[1 + (1 + kV2 k1) ]1=2kvk1;2; loc (1 + ?2 )1=2kwk1;2; loc : 2
From the trace theorem, this yields
k@v=@nk?1=2;?i (1 + ?2)1=2 1 + 1 + kV2 k1 2
!1=2
From the relations (100)-(101), we end up with
kvk1;2; loc (1 + ?2 )1=2 1 + 1 + kV2 k1 2
kvk1;2; loc: (101)
!1=2
khk1=2;?i
(102)
Inria
ΓV K
Convergence of Domain Decomposition Algorithms
31
Γi
Ωb
d/4
d/6
Bi
d/6
d/6
d/4
Mi d/6
Γb
Ω i
Figure 3: Description of the local domain loc and of the splitting used in the majoration of the global solution. and hence in particular
kvk0;2; loc (1 + ?2 )1=2 1 + 1 + kV2 k1 2
!1=2
khk1=2;?i :
(103)
Step 2 : Nearby local estimates
Let Ky = Bd=4 (y ) be the sphere centered on y and of radius d=4, with y belonging to ?V (see g. 3). By construction, ?V is the center surface of loc and i is subdomain of width 6d centered on ?V . Then by following the same argument as in step 2 of the proof of Theorem 6.2, we obtain : kvk1;Ky c1kvk0;2; loc ; (104) where c1 is a constant depending only on d; ; and . Let us consider now y1 ; : : :; yl in i such that
RR n 2435
32
Patrick Le Tallec , Moulay D. Tidriri
2i =
[ y2 i
B d6 (y)
[l j =1
Kyj = K:
By applying the relation (104) to each Kyj , we obtain
kvk1;K sup c1j kvk0;2; loc = c1kvk0;2; loc: j =1;:::;l
(105)
Step 3: Using the maximum principle For any Mi 2 i , we introduce (see g. 3) : a ball Bi centered on Mi and of radius d=6, the function vi = exp[k(r2 ? d2=36)]kvk1;@Bi. By construction of k (see previous x), '(r; k) is positive for all r 2 [0; d=6]. Then
by following the same argument as in step 3 of the proof of Theorem 6.2, we obtain the following inequality :
jv(Mi)j exp[?kd2=36]kvk1;@Bi: Step 4 : H 1 estimates Let us consider 2 H 1 ( loc), such that
(106)
= 1 in b ; (Figure 3) supp b [ i:
By choosing w = 2 v in (96) we obtain :
Z
loc
(? v + div(V v ) + (v= )) 2v = 0:
(107)
By using the Green's formula, as in step 4 of the previous theorem, we obtain :
kvk1;2; b[ i Then, if we take such that
Z
i
(v 2jr j2 + v 2V:r ):
kk0;2; i 1 Inria
33
Convergence of Domain Decomposition Algorithms
and
jj21;2; i = c2=d;
we obtain nally as in the proof of Theorem 6.2 :
q q 1=2 kvk1;2; b[ i kvk1; i c2=d 1 + kVk1 d=c2 : Step 5 : Getting the result Since @Bi K by construction, (105) and (106) imply kvk1; i exp(?kd2=36)c1kvk0;2; loc:
(108)
(109)
Furthermore by using the relation (103), it follows :
!
2 1=2 kvk1; i 1 + 1 + kV2 k1 c1(1 + 1= 2)1=2exp(?kd2 =36)khk1=2;?i:
By using the relation (108), we then obtain :
!
2 1=2 1 + k V k kvk1;2; b[ i 1 + 2 1 q q 1=2 c1 c2=d 1 + kVk1 d=c2 (1 + 1= 2)1=2exp(?kd2 =36)khk1=2;?i :
Before concluding we will establish an estimate of the quantity
k@v=@nk?1=2;?b : From (96), and by choosing
w 2 H 1 ( loc); with w = 0 on @ b \ @ i ; we obtain :
Z
b
(? v + div(V v ) + v= )w = 0:
After application of the Green's formula and from (64), we obtain :
RR n 2435
(110)
(111)
34
Patrick Le Tallec , Moulay D. Tidriri
Z @v Z w = (rv rw + (1= )V:rvw + 1 vw): @n ?b
b
By proceeding as in step 1 of the current proof, we obtain the following estimate :
! 2 1=2 1 + k V k k@v=@nk?1=2;?b (1 + 1= 2)1=2 1 + 2 1 kvk1;2; b:
(112)
The relations (111)-(112) nally yield the desired result.
6.3 Convergence of the explicit algorithm with t = 1 and small. Consider the following elliptic problem :
8 > > > < + V:r ? = 0 in ; = 1 on ?1 ; > > > : = 0 on ?b ;
(113)
that we want to solve by the fundamental algorithm (7)-(8) with "t = +1": The proposed algorithm can be written in this case as
set oloc = ol and o = o. then, for n 0; nloc and n being known, solve
8 n+1 > n+1 n+1 loc > > < + V:r loc ? loc = 0 in loc; nloc+1 = n on ?i ; > > > : nloc+1 = 0 on ?b ;
8 n+1 > + V:r n+1 ? n+1 = 0 in ; > > > < n+1 = 1 on ?1 ; > > n+1 > @ n+1 > : @ = loc on ?b : @n
(114)
(115)
@n
Inria
35
Convergence of Domain Decomposition Algorithms
We will show in this section that the algorithm described above converges. More precisely we have the following theorem.
Theorem 6.4 For suciently small, being the solution of the stationary pro-
blem (113), we have :
n+1 @ in H ?1=2(? ), i) @ @nloc converges towards @n b
ii) n+1 converges towards in H 1=2(?i ), iii) n+1 converges towards in H 1( ), iv) nloc+1 converges towards in H 1( loc). Proof:
In the rest of this section, we will take = . By the transformation n+1 ! reduce the problem to the case 1 = 0. After multiplication and integration by parts in (114), we obtain that nloc+1 = n on ?i and satis es :
n+1 ? with the solution of the stationary problem, we will
nloc+1 w + Z V:r n+1w + Z r n+1 rw loc loc
loc
loc
loc
Z
=
Z @ n+1 loc @n w; 8w 2 W: ?i
(116)
We can apply then Theorem 6.3 and obtain
q0 1 n+1 @ 0 2 loc k @n k?1=2;?b c1 c2=d 1 + 2 (1 + kV k1)
q 0 1=2 1 1 + kV k d=c
1
2
(1 + 1= 2)exp(?kd2 =36)k n k1=2;?i :
RR n 2435
(117)
36
Patrick Le Tallec , Moulay D. Tidriri
On the other hand, after multiplication by w and integration by parts in (115), we obtain the equality
Z n+1 Z Z Z @ n+1 n+1w + r n+1 rw = loc w + V: r @n w
?b
(118)
with w 2 H 1( ) and w = 0 on ?1 . By applying Theorem 6.2, we obtain :
k n+1 k
1=2;?i
q 1=2 1 c1 c2=d 1 + kV k1 d=c2 q
exp(?kd2=36)k@ nloc+1 =@nk?1=2;?b :
(119)
By using (117) and (119), we then get :
q k@ nloc+1=@nk?1=2;?b c01 c02=d(1 + 12 1 + kV k21) 1 q 1=2 1 + kV k1 d=c02 q 1=2 q 1 c1 c2 =d 1 + kV k1 d=c2 ! 2 d (1 + 1= 2 )exp ?k 18 k@ nloc =@nk?1=2;?b ; . Therefore for suciently small, the coecient of reduction will be with k = p dominated by the exponential term and will be then strictly less than 1, implying the linear convergence towards zero of
k@ nloc+1=@nk?1=2;?b : This is exactly assertion (i). By using (119), it also follows that n+1 tends to 0 in H 1=2(?i ). From (69) applied with g = @ nloc+1 =@n, we have in addition
k n+1 k1;2; co k@ nloc+1=@nk?1=2;?b ;
Inria
Convergence of Domain Decomposition Algorithms
37
and therefore k n+1 k1;2; converges to zero at the speed of k@ nloc+1 =@nk?1=2;?b . From (102) we also have
k n+1 k
(1 + 1= 2)1=2
1=2 1 2 1 + 2 (1 + kV k1 ) k n k1=2;?i
2 loc 1;2; loc and then k n+1 k1;2; also converges to zero at the speed of k n k1=2;?i .
6.4 Convergence of a xed point method for the implicit scheme.
The implicit scheme of section (5) couples the global and the local problem. To uncouple them, it is advisable to use the xed point algorithm below :
set oloc;0 = ol and o = o, being known, then for k 0; nkj+1 ?i solve
8 n+1 ? n > loc + div(v n+1 ) ? n+1 loc;k+1 > loc;k+1 = 0 in loc ; loc;k+1 > < t n+1 on ? ; (120) +1 nloc;k i > +1 = k > > : +1 nloc;k +1 = 0 on ?b ; 8 n+1 ? n > k+1 n+1 ) ? n+1 = 0 in ; > + div( v k+1 k+1 > < t +1 = (121) nk+1 1 on ?1 ; > > > : +1 =@n = @ n+1 =@n on ? : @ nk+1 b loc;k+1 We will study now the algorithm (120)-(121). By setting
RR n 2435
n+1 n+1 loc;k;q = loc;k+1 ? loc;q+1 ;
(122)
n+1 k;q = ( k ? nq +1 );
(123)
38
Patrick Le Tallec , Moulay D. Tidriri
we have that loc;k;q and k;q verify the following equations :
8 > > < > > :
loc;k;q =t + div (v loc;k;q) ? loc;k;q
8 > > > < > > > :
= 0 in loc ; loc;k;q = k?1;q?1 on ?i ; loc;k;q = 0 on ?b ;
(124)
k;q =t + div (v k;q) ? k;q
= 0 in ; k;q = 0 on ?1 ;
loc;k;q on ? : @@nk;q = @ @n b
(125)
If t is suciently small, we have from the previous section that k;q and loc;k;q +1 are Cauchy sequences converge linearly to zero. Hence the sequences nk +1 and nloc;k which converge linearly towards the unique solutions n+1 and nloc+1 of the implicit scheme. This guarantees the convergence of the above xed point algorithm.
Inria
Convergence of Domain Decomposition Algorithms
39
7 Numerical analysis of the stability of the algorithm (7)-(8) In this section we focus our attention on the numerical solution of the steady problem (2), using now the explicit time marching algorithm (7)-(8) studied in the previous paragraphs. We rst assume here that the boundary condition on ?b in (8) is explicit n+1
n
( @ @n = @ @nloc ) so that the resulting algorithm is parallel (Jacobi type). Here, denotes the domain surrounding the obstacle (an ellipse in our numerical study) as described in Figure 1. The global and local domains are discretized by fully overlapping compatible nite element grids. The global mesh contains 1378 nodes and 2662 elements (see gure 4). Further the time marching algorithm is being initialized by setting 0 to zero.
RR n 2435
40
Patrick Le Tallec , Moulay D. Tidriri
Figure 4: Description of the nite element mesh and of the local subdomain.
Inria
Convergence of Domain Decomposition Algorithms
41
In a rst step, the velocity eld is obtained by solving the following problem : div v = curl v = 0 (inviscid incompressible ow), v1 = (1; 0); v:n = 0 on the body?b; with a rst order nite element method using the same global mesh.
RR n 2435
42
Patrick Le Tallec , Moulay D. Tidriri
If we set v = 0, the algorithm may or may not converge depending on the values of t. More precisely, we observe that the algorithm converges linearly when t < 0 and is ndivergent otherwise. This is graphically shown in gures (4-7) where +1 ? n k k the values of tk 0 k1 are plotted versus the iteration count n for t taken on the values 10?6, 10?1, 1, and 10. Further, when the velocity is taken suciently large, the algorithm becomes unconditionnally stable. In particular, initializing our algorithm by 0 = 0 with kv1 k = 1, = 0:1 and t = 100 results in a converging scheme ( g. 8). By intuition such a behavior seems natural. An overestimation of the solution n at the interface ?i implies an overestimation of the friction forces on ?b . For suciently small time steps, this overestimation will not a ect the value of n+1 on ?i and can therefore be ignored at the next time step. If the Reynolds is suciently large, this error will only a ect the wake region but will not have any in uence at the interface ?i . To the contrary, for large t and small , the introduced error does a ect the value of n+1 on ?i . The in uence of the error on n+1 may be ampli ed throughout the iteration process.
Inria
Convergence of Domain Decomposition Algorithms
43
Another variant of the algorithm consists of replacing the explicit Dirichlet condition
nloc+1 = n on ?i in the algorithm (7)-(8)
by the following semi-implicit condition
nloc+1 = n+1 on ?i:
In fact, this implies replacing the previously parallel algorithm (Jacobi like ) by the sequential algorithm (Gauss-Seidel like ). When we solve the pure di usion problem (i.e. with ow velocity v = 0) with = 1 and t = 1 (respectively t = 2) we obtain a better convergence history : the speed of the new algorithm is linear and clearly faster than the parallel algorithm. the domain of convergence is moderately larger (see table 1). To study experimentally in more details the convergence behavior of both algorithms we assume that we have a linear behavior of our algorithm, and hence that the error at the iteration n will satisfy the following inequality
k n+1 ? nk1 K nk 1 ? 0k1:
The algorithm converges if K < 1. An estimate for K can be found by considering as in table 1 the ratio n+1 n ? n1 log k k 1 ?? 0k k1 = ?log K
1
which is displayed as a function of ( t) for n = 14 and di erent values of V = v . A negative value of this ratio means divergence of the algorithm. As expected, this ratio is positive for suciently small values of t and converges towards zero as t goes to zero. In this table, we observe that for V = 0; t < 0 2, the algorithm converges. However the convergence is slow since the minimal contraction constant Kmin (for the optimal value of t) is close to one (see table 2). For V = 10, then the algorithm appears to converge for a much larger range of values of t and the optimal contraction constant is much smaller. This is summarized on table 2 where we have displayed the best possible contraction constants for each of the coupling algorithms and di erent values of the Reynolds V = v .
RR n 2435
44
Patrick Le Tallec , Moulay D. Tidriri
t 1/1000 1/10 1/2 2 Gauss-Seidel 0.06 0.1 0.22 0.5 V =0 Jacobi 0.1 0.22 0 V =0 Gauss-Seidel 0.03 0.25 1.46 2.12 V = 10 Jacobi 0.03 0.28 1.15 1.15 V = 10 Jacobi 0.23 2.79 2.8 2.7 V = 1000
5 10 50 1000 -0.27 -0.5 -0.75 -0.8 -0.09 -0.25 -0.4 -0.41 2.8
2.6
2.4
2.4
1.15
1.15
1.14
1.14
2.75
2.8
-
-
Table 1: Contraction constant (in fact minus its logarithm) in function of t for the explicit (Jacobi) and semi-explicit (Gauss-Seidel) version of our coupling algorithm. We observe divergence for V = 0 and t > 2 and convergence otherwise.
Jacobi (parallel) V Kmin 0 0:85 10 0:50 103 0:14
Gauss-Seidel (sequential) V Kmin 0 0:68 10 0:11
Table 2: Minimal contraction constant versus the Reynolds V for both sequential and parallel versions of the algorithm.
Inria
45
Convergence of Domain Decomposition Algorithms
Figure 5: Convergence of the Time Marching Algorithm: k tk ?0 k 1 k are plotted versus the iteration count n for t = 10?6 , v = 0 (Jacobi). Observe the very slow convergence. n+1
RR n 2435
n
46
Patrick Le Tallec , Moulay D. Tidriri
Figure 6: Convergence of the Time Marching Algorithm: k tk ?0 k 1 k are plotted versus the iteration count n for t = 10?1 , v = 0 (Jacobi). n+1
n
Inria
47
Convergence of Domain Decomposition Algorithms
Figure 7: Convergence of the Time Marching Algorithm: k tk ?0 k 1 k are plotted versus the iteration count n for t = 1, v = 0 (Jacobi). n+1
RR n 2435
n
48
Patrick Le Tallec , Moulay D. Tidriri
Figure 8: Divergence of the Time Marching Algorithm: k tk ?0 k 1 k are plotted versus the iteration count n for t = 10, v = 0 (Jacobi). n+1
n
Inria
49
Convergence of Domain Decomposition Algorithms
Figure 9: Convergence of the Time Marching Algorithm: k tk ?0 k 1 k are plotted versus the iteration count n for t = 10 and the ow velocity is equal to 1 (Jacobi). n+1
RR n 2435
n
50
Patrick Le Tallec , Moulay D. Tidriri
8 Conclusion We have analysed the convergence properties of a standard time marching algorithm for solving a domain decomposed advection-di usion problem. We were able to prove theoretically the unconditional stability and linear convergence of the fully implicit algorithm (x5). Using the maximum principle, we were also able to prove in x6 the linear convergence of the semi explicit algorithm used with t = +1 on the operator Lu = ?u + div(vu) + u if was suciently large, and the linear convergence of the xed point algorithm applied at each time step to the solution of the implicit problem, independently of any discretization strategy. When using the uncoupled semi-explicit algorithm in the general case, numerical evidence indicate that this algorithm is unstable for large values of t and small overlapping, and that it becomes linearly convergent when t is below a Reynolds dependent threshold (x7). This conditional stability is not a real issue for practical CFD problems because most solvers already require to use small time steps inside each domain. Nevertheless, it would be nicer to derive an uncoupled unconditionnally stable version of the present time marching algorithm.
Inria
Convergence of Domain Decomposition Algorithms
51
References [1] R. Glowinski, G. Golub and J. Periaux (eds), Proceedings of the First International
Symposium on Domain Decomposition Methods for Partial Di erential Equations, Paris, France, January 7-9, 1987 , (SIAM, Philadelphia, 1988). [2] T. Chan, R. Glowinski, J. Periaux and O. Widlund (eds), Proceedings of the Second International Symposium on Domain Decomposition Methods for Partial Di erential Equations, Los Angeles, California, January 1988 , (SIAM, Philadelphia, 1989). [3] T. Chan, R. Glowinski (eds), Proceedings of the Third International Symposium on Domain Decomposition Methods for Partial Di erential Equations, Houston, Texas, March 20-22, 1989 , (SIAM, Philadelphia, 1990).
[4] R. Glowinski, Y. Kuznetsov, G. Meurant, J. Periaux and O. Widlund (eds),Proceedings of the fourth international symposium on Domain Decomposition Methods for Partial Di erential Equations, Moscou, June 1990 , (SIAM, Philadelphia, 1991). [5] T. Chan, D. Keyes, G. Meurant, S. Scroggs and R. Voigt (eds), Proceedings of
the fth international symposium on Domain Decomposition Methods for Partial Di erential Equations, Norfolk, May 1991 , (SIAM, Philadelphia, 1992). [6] A. Quarteroni (ed), Proceedings of the sixth international symposium on Domain Decomposition Methods for Partial Di erential Equations, Como, June 1992 ,
(AMS, Providence, 1994). [7] Y. Achdou and O. Pironneau, A fast solver for Navier-Stokes Equations in the laminar regime using mortar nite element and boundary element methods, Technical Report 93-277 (Centre de Mathematiques Appliquees, Ecole Polytechnique, Paris, 1993) [8] A. D. Aleksandrov, Majoration of solutions of second order linear equations, Vestnik Leningrad univ. 21, 5-25(1966) English translation in Amer. Math. Soc. Transl. (2) 68, 120-143(1968). [9] I. Ya. Bakel'man, Theory of quasilinear elliptic equations. Siberian Math. J. 2, 179-186(1961).
[10] J. Bramble, R. Ewing, R. Parashkevov and J. Pasciak, Domain decomposition methods for problems with partial re nement SIAM J. Sci. Stat. Comp. 13 (1992) 397{410. [11] C. Canuto and A. Russo, On the Elliptic-Hyperbolic Coupling. I: the Advection Di usion Equation via the -formulation, Math. Models and Meth. Appl. Sciences (to appear) (1993). [12] D. Gilbarg, N. S. Trudinger, Elliptic partial di erential equations of second order.Berlin-Heidelberg-New York, Springer Verlag 1977.
RR n 2435
52
Patrick Le Tallec , Moulay D. Tidriri [13] W. Gropp and D. Keyes, Domain Decomposition Methods in Computational Fluid Dynamics, Int. J. Num. Meth. Fluids 14 (1992) 147-165. [14] P. Le Tallec and M. D. Tidriri, Coupling Navier-Stokes and Boltzmann. Submitted to J. Comp. Phy. [15] L. Marini and A. Quarteroni, An iterative procedure for domain decomposition methods: a Finite Element approach, in [1]. [16] L.D. Marini, A. Quarteroni, A relaxation procedure for domain decomposition methods using Finite Elements, Numer. Math. 55, (1989) 575{598. [17] A. Quarteroni, G. Sacchi Landriani and A Valli, Coupling of Viscous and Inviscid Stokes Equations via a Domain Decomposition Method for Finite Elements, Technical report UTM89-287 (Dipartimento di Mathematica, Universita degli Studi di Trento, 1989). [18] M. D. Tidriri, Couplage d'approximations et de modeles de types di erents dans le calcul d'ecoulements externes,these, Universite de Paris IX, May 1992. [19] M. D. Tidriri, Domain Decomposition for Incompatible nonlinear models. INRIA Research Report. To appear.
[20] M. D. Tidriri, Domain decomposition for compressible Navier-Stokes equations with di erent discretizations and formulations. To appear in J. Comp. Phy.
Inria
Convergence of Domain Decomposition Algorithms
53
Appendix The main result of this section relies on the notion of a contact set. If u is a continuous arbitrary function on , the upper contact set, denoted ?+ or ?+u , is the sub-set of , de ned by ?+ = fy 2 ; 9p(y) 2 IRn such that u(x) u(y) + p (x ? y) 8x 2 g: (126) We see that u is a concave function on i ?+ = . When u 2 C 1( ) we must have p = Du(y) in the relation (126). In addition, when u 2 C 2( ), the Hessian matrix D2 u = [Dij u] is negative on ?+ . In general, ?+ is closed in . If u is a continuous arbitrary function on , we de ne the "normal mapping" (y) = u (y) at point y 2 by (y) = fp 2 IRn; u(x) u(y) + p:(x ? y) 8x 2 g: (127) We can see that (y) is non empty i y 2 ?+ . In addition when u 2 C 1( ), we have (y) = Du(y) on ?+ ; in other words is the gradient eld of u on ?+ .
As a particular case of the Bakelman-Alexandrov ([8] and [9]) maximum principle, we have under the above notation. Lemma .1 For u 2 C 2 ( ) \ C o( ), we have :
d kaij Dij u=D kn;?+
nwn1=n with d the diameter of and wn the volume of a unit sphere in Rn .
sup u sup@ u +
For further details see [12]. We now proceed to the proof of Theorem 6.1, by following the steps of [12]. We take B^ = B1 (0) and the general case will be deduced by considering the coordinate transform, x ! x^ = (x ? y)=2R. We will begin, in rst step, by showing this result for u 2 C 2( ) \ W 2;n( ) and then in a second step we will deduce the result for u 2 W 2;n( ).
Step 1:
We suppose that u 2 C 2 ( ) \ W 2;n( ). For 1, we consider the cut o function de ned by (^x) = (1 ? jx^j2) : By di erentiation, we obtain D^ i = ?2 x^i (1 ? jx^j2) ?1 ;
RR n 2435
54
Patrick Le Tallec , Moulay D. Tidriri
By setting
D^ ij = ?2 ij (1 ? jx^j2) ?1 + 4 ( ? 1)^xi x^j (1 ? jx^j2) ?2 : v = u;
we then obtain a^ij D^ ij v = a^ij D^ ij u + 2^aij D^ i D^ j u + ua^ij D^ ij (f^ ? ^biD^ i u ? c^u) + 2^aij D^ i D^ j u + ua^ij D^ ij : Let ?+ = ?+v be the upper contact set v, in the sphere B^ ; we have : u > 0 on ?+ :
If x 2 @ B^ such that p:(x ? y) < 0 we indeed have v(x) = 0. Consequently v(y) + p:(x ? y) v(x) = 0:
Moreover, using the concavity of v on ?+ , we can estimate the following quantity : ^ j = (1=)jDv ^ ? uD ^ j: jDu Indeed,
^ (1=)( Dv ^ + u D ^ ) Du
^ j) (1=)( 1 ?v jx^j + ujD 2(1 + )?1= u:
In that way, we have on ?+ the following inequality : ?a^ij D^ ij v f(16 2 + 2 )^ ?2= + 2 j^bj?1= + c^gv + jf^j: Since c^ 0, we deduce the inequality ?a^ij D^ ij v f(16 2 + 2 )^ ?2= +
(128)
2 j^bj?1= gv + jf^j
c1 ^?2= v + jf^j; with c1 = c(n; ; ; ^) independent of c^.
Inria
Convergence of Domain Decomposition Algorithms
55
Consequently, by applying Lemma 6.1 on B^ , we obtain, for 2 : ^ sup v ( d1=n )( ^1 )kc1 ^ ?2= v + jf^jkn;B^ : nwn D B^ By using the relation (6.1), it comes ^ ^ sup v ( d1=n )c1 k?2= vkn;B^ + ( d1=n )( ^1 )kf^kn;B^ nwn nwn B^ c1 d^(k?2= vkn;B^ + (1=^)kf^kn;B^ )
c1 d^(k?2= v+ kn;B^ + (1=^)kf^kn;B^ ) c1 d^((supv+ )1?2= k(u+ )2= kn;B^ + (1=^)kf^kn;B^ );
where c1 is a constant depending only on n; ; and ^. Here, d^ is the diameter of B^ (d^ = 2). By using the Young inequality under the form ab "aq + "?r=q br
for q = (1 ? 2= )?1 and r = =2, we have
2 (supv+ )1?2= k(u+ )2= kn;B^ "supv+ + "1? =2 k(u+ )2= k = n;B^ ; 8" > 0:
By taking " = 1 ^ and plugging in our inequality on v, we obtain : 2c1 d 2 sup v (1=2) sup v+ + (1=2)1? =2 (c1 d^) =2 k(u+ )2= k = n;B^
B^
(129)
^ ^)kf^kn;B^ : +(c1 d= We want to prove the theorem for all p > 0. We will treat separately the cases p n and p > n. If p n, we set = 2n=p. In this case we have 2 + k(u+ )2= k = n;B^ = k(u )kp;B^ :
Plugging this in our inequality on v, we obtain : ^ ^)kf^kn;B^ : (1=2) sup v (1=2)1? =2 (c1 d^) =2 k(u+ )kp;B^ + (c1 d= B^
Consequently, we obtain the following inequality ;
RR n 2435
56
Patrick Le Tallec , Moulay D. Tidriri
Z
^ 2^)kf^kn;B^ g: sup v c2 f( ^ (u+ )p )1=p + (d= B^
B
On the sphere B1=2 (0), the cut o function satis es 1= (1=2) : It follows, then sup u
B1=2 (0)
sup (v= )
B1=2 (0)
2 sup v: B^
Finally we end up at the desired estimate
Z
^ 2^)kf^kn;B^ g: sup u c3 f( (u+ )p )1=p + (d= ^
B1=2 (0)
B
for u in W 2;n( ) \ C 2( ). The constant c3 above depend only on n; ; and ^, but is independent of c^. On the other hand if p > n, we have : 2n= < p; 8 2: Then, it follows (by assuming 2) jB^ j?1=(2n= )k(u+ )k2n= ;B^ jB^ j?1=pku+ kp;B^ : But 2 ku+ k2n= = k(u+ )2= k = n;B^ ;
and therefore, by processing as before, we obtain the desired estimate
Z
^ 2^)kf^kn;B^ g sup u c4 f( (u+ )p )1=p + (d=
B1=2 (0)
B^
for u in W 2;n( ) \ C 2( ). The constant c4 above depends only on n; ; and ^, but is independent of c^. Transformation x^ ! x. By construction, D^ ij = R?2Dij , thus ^ = R?2 and ^ = R2. In addition, we have jB j = wn (2R)n and jgjp;B^ = R?n=pjgjp;B .
Inria
Convergence of Domain Decomposition Algorithms
57
Written in term of x, the last inequality becomes n
sup u c4 f( 2jBwjn
BR(y)
Z
B
1=n (u+ )p dx)1=p + ( 2wn R )kf kn;B g;
with c4 a function of n; ; ^ = R2 and p. This is the desired estimate for u 2 W 2;n( ) \ C o ( ).
Step 2:
Now, let u 2 W 2;n( ). By density, let (um ) be a sequence of functions of C 2(B ), converging towards u in W 2;n(B ). The injection of W 2;n(B ) in C o (B ) is continuous, consequently (um ) converges uniformly towards u in B . We have Lum = L(um ? u) + Lu
f + L(um ? u): By setting, fm = L(um ? u), we observe by construction that fm converges towards 0 in Ln ( ). As um 2 W 2;n( ) \ C 2( ) and f~m = f + fm is in Ln ( ), the estimate (59) is valid also for um , so that we have
Z sup um ctef( jB1 j (u+m )p )1=p + R kf~kn;B g: B (y) B R
(130)
Using previous results and taking the limit, we have : Z sup u ctef( jB1 j (u+ )p )1=p + R kf kn;B g: BR (y)
B
Observe also that by replacing u by ?u, the theorem can be extended easily to the case of supersolutions and solutions of the equation : Lu = f:
RR n 2435
Unité de recherche Inria Lorraine, Technopôle de Nancy-Brabois, Campus scientifique, 615 rue du Jardin Botanique, BP 101, 54600 Villers Lès Nancy Unité de recherche Inria Rennes, Irisa, Campus universitaire de Beaulieu, 35042 Rennes Cedex Unité de recherche Inria Rhône-Alpes, 46 avenue Félix Viallet, 38031 Grenoble Cedex 1 Unité de recherche Inria Rocquencourt, Domaine de Voluceau, Rocquencourt, BP 105, 78153 Le Chesnay Cedex Unité de recherche Inria Sophia-Antipolis, 2004 route des Lucioles, BP 93, 06902 Sophia-Antipolis Cedex
Éditeur Inria, Domaine de Voluceau, Rocquencourt, BP 105, 78153 Le Chesnay Cedex (France) ISSN 0249-6399
Convergence Analysis of Domain Decomposition algorithms with full overlapping for the advection-diffusion problems. Patrick Le Tallec , Moulay D. Tidriri
N˚ 2435 Octobre 1994
PROGRAMME 6 Calcul scientifique, modélisation et logiciel numérique
ISSN 0249-6399
apport de recherche
1994
Convergence Analysis of Domain Decomposition algorithms with full overlapping for the advection-di usion problems. Patrick Le Tallec , Moulay D. Tidriri Programme 6 | Calcul scienti que, modelisation et logiciel numerique Projet MENUSIN Rapport de recherche n 2435 | Octobre 1994 | 57 pages
Abstract: The aim of this paper is to study the convergence properties of a Time
Marching Algorithm solving Advection-Di usion problems on two domains using incompatible discretizations. The basic algorithm is rst presented, and theoretical or numerical results illustrate its convergence properties. This study is based on spectral theory, a priori estimates and a Di-Giorgi-Nash maximum principle . Key-words: non symmetric elliptic operators, domain decomposition, incompatible grids, advection-di usion, time marching algorithms, overlapping, maximun principle. (Resume : tsvp) Universite Paris-Dauphine and INRIA, Domaine de Voluceau- Rocquencourt- B.P. 105- Le
Chesnay Cedex (France) email: [email protected] Department of Mechanical Engineering, Yale University, New Haven, CT 06520 email: [email protected]. This work has been done in 1991 at INRIA and has been supported by the Hermes Research program under grant number RDAN 86.1/3. The author is presently supported by the National Science Foundation under contract number ECS-8957475 and by the United Technologies Research Center.
Unité de recherche INRIA Rocquencourt Domaine de Voluceau, Rocquencourt, BP 105, 78153 LE CHESNAY Cedex (France) Téléphone : (33 1) 39 63 55 11 – Télécopie : (33 1) 39 63 53 30
Analyse de Convergence des Methodes de Decomposition de Domaines avec Recouvrement total Resume : Le but de l'article propose est d'etudier les proprietes de convergence d'un algorithme de marche en temps applique a la solution d'un probleme d'advection di usion pose sur deux sous-domaines recouvrants. Les discretisations utilisees sur chaque sous-domaine peuvent ^etre totalement independantes. L'article rappelle les principes de l'algorithme utilise, en presente une analyse theorique, et une illustration numerique. L'etude theorique s'appuie sur des estimations apriori, sur le principe du maximun de Di Giorgi-Nash, et sur une analyse spectrale unidimensionnelle. Mots-cle : operateurs elliptiques nonsymetriques, decomposition de domaines, grilles incompatibles, advection-di usion, algorithmes de marche en temps, recouvrement, principe du maximum.
Convergence of Domain Decomposition Algorithms
3
1 Introduction Domain decomposition methods have recently appear as an ecient strategy for solving large scale problems on parallel computers ([1], [2], [3], [4], [5], [6]). Nevertheless, they can also be used in order to couple di erent models [11], [17]. In this paper we will examine a domain decomposition strategy which can be applied to such situations. This approach was introduced in order to solve several diculties that occur in
uid mechanics. In particular, our aim is to introduce several sub-domains in order to do one of the following :
Solve di erent problems on each subdomain. Use di erent kinds of approximation methods on each subdomain [7]. Use \local re nement techniques" or \mesh adaptive techniques", locally, per subdomain ([10]).
The subdomains fully overlap and the coupling is achieved through \friction" forces acting on the internal boundary of the domain, these friction forces being updated by an explicit time marching algorithm. The theoretical study of our method will be done on an Advection-Di usion problem, which will serve as our model problem from now on. The analysis will be made at the continuous level, independently of any discretization strategy, which means that the derived results will be mesh independent. In the next section we will describe this model problem. In the third section we will present our algorithm for some special cases. The fourth section will treat the one-dimensional stationary problem. We will show from this study that we can improve the convergence of this method by introducing a relaxation parameter (see [5]) and by replacing one of the Dirichlet condition by a Neumann type condition. The fth section will focus on the resulting multidimensional case. We rst study the proposed strategy for the time implicit problem. We then establish a priori estimates and show from these estimates exponential convergence of our algorithm. A detailed study of the semi-discrete algorithm will be done in section 6. In the last section we study the numerical stability of the explicit algorithm. Practical applications of the proposed algorithm to real life CFD problems can be found in [14], [18], [19], and [20].
RR n 2435
4
Patrick Le Tallec , Moulay D. Tidriri
Γe Ω
Γι Γb
ΩV
Figure 1: description of computational domain.
2 The model problem Consider a bounded domain, of IRn such that its boundary @ is lipschitzian and
loc = v a connected domain of IRn with loc ( g. 1). Let : ?b = @ \ @ loc ; ( internal boundary) ?i = @ loc \ ; ( interface) ?1 = @ n?b :(far eld boundary) We denote by n the external unit normal vector to @ or @ loc . Let v a velocity eld given such that :
8 < divv = 0 in ; : v:n = 0 on ?b:
(1)
The model problem we want to solve is the following:
Inria
Convergence of Domain Decomposition Algorithms
5
Find ', a real valued function, de ned on and satisfying :
8 > div(v') ? ' = 0 in ; > < ' = '1 on ?1 ; > > :
(2)
' = 0 on ?b :
Above v is the ow velocity and is the di usion coecient. Most CFD algorithms will in fact consider the solution of this problem as the stationary solution of the evolution problem (3) described below : Find : (0; T ) ! IR such that,
8 @ > > @t + div(v ) ? > > < > > > > : (0)
= 0 in (0; T ); = 1 on ?1 (0; T ); = 0 on ?b (0; T );
(3)
= 0 in : The general CFD algorithm consists then in integrating (3) in time until reaching a stationary solution.
3 General Algorithm 3.1 Time-Continuous case
Let us introduce the local sub-domain loc (see g. 1) which has as external boundary ?i , and let us consider the trace loc of on the sub-domain loc , as an independent variable, to which we associate an arbitrary independent initial value ol 6= o j loc . We now replace the evolution problem (3) by the following evolution system : Find (resp. loc ) : ! IR (resp. loc ! IR) satisfying :
8 @ > + div(v ) ? = 0 in (0; T ); > > < @t = 1 on ?1 (0; T ); > > > : @ = @ loc on ?b (0; T ); @n
RR n 2435
@n
(4)
6
Patrick Le Tallec , Moulay D. Tidriri
8 @ > loc + div(v ) ? > loc loc = 0 in loc (0; T ); > < @t loc = 0 on ?b (0; T ); > > > : loc = on ?i (0; T );
(5)
(0) = o in ; loc(0) = ol in loc :
(6)
Remark 3.1 The global problem (4)-(6) has no no-slip boundary condition. This
suppresses the boundary layer which appears at low viscosity and facilitates the numerical solution of this problem. The boundary layers are modeled by the local problems (5) (6) which are only to be solved on a small domain loc , with a very ne discretisation if needed. The two problems are only coupled by their boundary conditions and not by volumic interpolation.
Remark 3.2 The \Neumann" version of this method consists in replacing in (5) the boundary condition on ?i by the corresponding Neumann condition :
@ loc = @ on ? : i @n @n
3.2 Time Discrete case
The general algorithm that we propose for the solution of our model problem (2) is as usual to integrate in time the evolution problem (4)-(5)-(6) until we reach a stationary solution. This integration in time is then achieved by the following uncoupled semi-explicit algorithm, where the operators are treated implicitely inside each subdomain and where the coupling boundary conditions are treated explicitely :
set 0loc = ol and 0 = 0; then, for n 0; nloc and n being known, solve successively 8 n+1 n loc ? loc + div(v n+1 ) ? n+1 = 0 in ; > > loc loc loc > t < nloc+1 = n on ?i ; > > : n+1 = 0 on ?b ;
(7)
loc
Inria
7
Convergence of Domain Decomposition Algorithms
8 n+1 n > ? + div(v n+1 ) ? n+1 = 0 in ; > > > < t n+1 = 1 on ?1 ; > > n+1 n+1 > > : @ = @ loc on ?b : @n
(8)
@n
Remark 3.3 We have a full decoupling between (7) and (8). They can be discretized and solved by two independent solution techniques.
Remark 3.4 The "Neumann" version of this method consists in replacing the Dirichlet boundary condition on ?i by the corresponding Neumann condition :
@ nloc+1 = @ n on ? : i @n @n
Remark 3.5 The fully implicit version of this method consists in replacing the condition :
nloc+1 = n on ?i
@ nloc+1 = @ n on ? ) by the condition : i @n @n nloc+1 = n+1 on ?i ; @ n+1 @ n+1 on ? ). (resp. loc = i @n @n (resp.
The two subproblems are then coupled at each time step.
Remark 3.6 If we replace in (8) by E de ned as follows :
E = n loc; and ?b by ?i , and if we set t = 1, we obtain a nonoverlapping version of our strategy, which is a standard Dirichlet-Neumann algorithm [15], [16].
Remark 3.7 The initial condition ol is not assumed to be equal to 0 on the local
subdomain loc because in most cases this condition is impossible to impose at the discrete level since the grid used on loc will be in general di erent from the grid used on . In addition, even if we assume ol = 0 , we will not have nloc = n on
loc unless we use the fully implicit algorithm on compatible grids.
RR n 2435
8
Patrick Le Tallec , Moulay D. Tidriri
4 Stationary one-dimensional case
For t = +1, the above algorithm can be written :
set 0loc = 0 and 0 = 0, then, for n 0; nloc and n being known, solve 8 > div(v nloc+1 ) ? nloc+1 = 0 in loc ; > < nloc+1 = n on ?i; > > : n+1 loc
(9)
= 0 on ?b ;
8 > div(v n+1 ) ? n+1 = 0 in ; > > < n+1 = 1 on ?1 ; > > n+1 n+1 > loc on ?b : : @ @n = @ @n
(10)
In one space dimension, if we take as domain the interval ]0; 1[ of IR decomposed into two sub-domains =]0; 1[ and loc =]h2; 1[ then we obtain
0 < h2 < 1;
8 (n)0 > v'2 ? '(2n)" > > < '(2n) (h2) > > > : '(2n) (1) 8 (n)0 > v'1 ? '(1n)" > > < '(1n) (0) > > 0 > : '1(n) (1)
(11)
= 0 on ]h2 ; 1[; = '(1n?1)(h2 );
(12)
= b; = 0 on ]0; 1[; (13)
= a; 0
= '2(n) (1): By introducing relaxation parameters, we can also introduce the following variants of the above algorithm :
Inria
9
Convergence of Domain Decomposition Algorithms
8 (n)0 > v'2 ? '(2n)" = 0 on ]h2 ; 1[; > > < '(2n)(1) = b; > > > : '(2n) (h2 ) = 2 '(1n?1)(h2) + (1 ? 2 )'(2n?1) (h2); 8 (n)0 > v'1 ? '(1n)" = 0 on ]0; 1[; > > < '(1n) (0) = a; > > 0 0 0 : '1(n) (1) = 1'2(n) (1) + (1 ? 1 )'1(n?1) (1):
(14)
(15)
We will now try to produce conditions under which algorithm (15)-(14) converges, and those for which this convergence is optimal. For this purpose, we write the interface solution under the form 0 '1(n) (1) = 0(1) + n; '(2n) (h2) = (h2) + n ;
(16) (17) where is the converged solution we are looking for. Using the analytical solutions of problems (14) (15), we obtain the following induction formula
with
n
n
!
= MIN
(n?1)
(n?1)
!
:
0 1 ?( v ) v h 1 ? 2 2 e (e 2 ? 1) B CC v ? v v MIN = B B @ 1v(1 ? 2)( v ) e( )v12(e( )h2 ? 1) + (1 ? ) CA 1 e( )(h2?1) ? 1 (e( )(h2 ?1) ? 1)
(18)
(19)
This iterative process converges if the spectral radius of matrix MIN is less than 1. A direct but tedious calculation then yields :
Lemma 4.1 The spectral radius of the transfer matrix of algorithm (15)-(14) is : p (MIN ) = max[ 12 jD D2 ? 4Rj] (20) with
RR n 2435
10
Patrick Le Tallec , Moulay D. Tidriri
D = 2 ? (1 + 2 ) + 12 e(?v=) (e(v=)h2 ? 1) e?v= ev(1h2 =) ? 1
(21)
R = (1 ? 1 )(1 ? 2 ):
(22)
From this calculation, we can get the following results :
i) When h2 goes to 1 (nonoverlapping), D goes to +1, and then, (MIN ) goes to +1. There is no-convergence at this limit. ii) The optimal convergence is obtained in the case where all the eigenvalues of the matrix MIN are zero, i.e., when : D = 0 and R = 0. The latter conditions imply in particular
1 = 1 or 2 = 1: If we choose, in addition, 1 = 2 , the condition D = 0 implies h2 = 0. In this case the sub-domain loc is equal to the whole domain, and the associate algorithm has no interest. iii) The convergence of the method depends symmetrically on both relaxation parameters. From ii), it is reasonable to take one of the i equal to 1 and call the other . By setting : (?v= ) (v= )h2 A = 1 ? e e(v=()(eh2 ?1) ??1 1) (23) we have then
(MIN ) = j1 ? Aj:
In this case, we get the following convergence results :
Theorem 4.1 1) There is optimal convergence (in 1 iteration) if (?v= ) (v= )h2 = opt = f1 ? e e(v=()(eh2 ?1) ??1 1) g?1 < 1:
(24)
(25)
2) The algorithm converges for all in ]0; A2 [.
Inria
11
Convergence of Domain Decomposition Algorithms
Corollary 4.1 1) The case without relaxation ( = 1) converges only if : 2 1; A i.e., by setting d = 1 ? h2 (overlapping length), only if : d v Log (1 + e2?v= ) (stability condition). 2) When v goes to zero, we must have d 21 . Remark 4.1 This theorem states that the application of algorithm (15)-(14) to the steady case converges only if the overlapping d is suciently large. In the same situation, we will see that the unsteady approach will converge to the same steady solution but with less restrictions on d. This motivates the introduction of the time marching algorithm of section 3. Moreover, this time marching technique is well adapted to nonlinear problems such as those encountered in uid mechanics.
5 Implicit time discretization
5.1 The general algorithm
This section deals with the convergence analysis of the proposed algorithm in multiple dimension when one uses the fully implicit version of our strategy (4)-(6) : Set 0loc = ol and 0 = 0; then, for n 0; nloc and n being known, solve
8 n+1 n ? + div(v n+1 ) ? n+1 = > > > > < t n+1 = > > n+1 > > : @ @n = 8 n+1 n loc ? loc + div(v n+1 ) ? n+1 > > loc loc > t < nloc+1 > > > : n+1 loc
RR n 2435
0 in ;
1 on ?1 ; n+1 loc on ?b ; @ @n = 0 in loc ; = n+1 on ?i ; = 0 on ?b :
(26)
(27)
12
Patrick Le Tallec , Moulay D. Tidriri
5.2 Convergence analysis
We have the following convergence result : Theorem 5.1 The solution of algorithm (26)-(27) converges linearly in H 1( ) towards the solution of the stationary problem (2), for all values of t and all choices of loc .
Proof Step 1: Local L2 estimates of ? loc
Substracting (27) from (26), multiplying the result by n+1 ? nloc+1 and integrating by parts, we obtain the following relation : Z 1 Z 1 2 n +1 n +1 ( ? loc ) ? ( n ? nloc )( n+1 ? nloc+1 )
loc t
locZ t (28) n +1 2 n +1 + jr( ? loc )j = 0:
loc
By using Cauchy-Schwarz, we deduce 1 k n+1 ? n+1 k2 n+1 ? n+1 j2 loc 1;2; loc loc 0;2; loc + j 2t (29) 1 n n 2 2t k ? lock0;2; loc : From the Poincare inequality, it follows that (30) k n+1 ? nloc+1k20;2; loc 1 + 21 tc k n ? nlock20;2; loc where c is the Poincare constant. By setting co = k o ? oloc k20;2; loc , we get by induction our basic L2 estimate k n+1 ? nloc+1 k20;2; loc ( 1 + 21 tc )n+1co: (31) By summing (29), we also obtain
k n+1 ? nloc+1k20;2; loc + 2 t
n X i=p
j i+1 ? iloc+1j21;2; loc
k p ? ploc k20;2; loc :
(32)
Inria
13
Convergence of Domain Decomposition Algorithms
In particular, we obtain for p = 0
k n+1 ? nloc+1k20;2; loc + 2 t
n X i=0
j i+1 ? iloc+1j21;2; loc
k o ? oloc k20;2; loc :
(33)
Step 2: Local H 1 estimates of n ? nloc Substracting the two rst equations in (26) and (27), multiplying the result by
xn = (
n+1 ? n+1 ) ? ( n ? n ) loc ; loc
and integrating on loc , we obtain 0 =
Z Z
loc
+
?
Z
jxnj2 +
loc
Z
t
(34)
div(v ( n+1 ? nloc+1 ))xn
(35)
loc ( n+1 n+1 )
r
?
loc
rxn
@ ( n+1 ? n+1 )xn: loc @ loc @n
The use of the third relation in (26) and the second relation in (27) then yields
Z
loc
jxnj2 +
Z
loc
div(v ( n+1 ? nloc+1 ))xn +
Z
loc
r( n+1 ? nloc+1 )rxn = 0:
By using now Cauchy-Schwarz, we obtain
Z
loc
jxnj2 +
Z
2 r( n+1 ? nloc+1)rxn kv2k1 j n+1 ? nloc+1 j21;2; loc
loc + 1 kxn k20;2; loc : 2
By construction of xn and Cauchy-Schwartz, the last inequality becomes
RR n 2435
14
Patrick Le Tallec , Moulay D. Tidriri
kvk21 j n+1 ? n+1j2 1 kxn k2 + j n ? nloc j21;2; loc 1 ; 2 ;
0 ; 2 ;
loc loc loc 2 2 2t
(36)
? 2 t j n+1 ? nloc+1j21;2; loc :
By using again relation (29), it follows that 2 kxnk20;2; loc k2vk1t (k n ? nlock20;2; loc ? k n+1 ? nloc+1 k20;2; loc ) + t (j n ? nloc j21;2; loc ? j n+1 ? nloc+1 j21;2; loc ) t (G(n) ? G(n + 1))
(37)
with G de ned by
2 G(n) = k2vk21 k n ? nloc k20;2; loc + j n ? nloc j21;2; loc : From (37), G is a decreasing function. This last property then yields
(n + 2 ? p)G(n + 1)
nX +1 i=p nX +1 i=p
(38)
G(i)
X
2 n+1
j i ? ilocj21;2; loc + k2vk21
i=p
k i ? ilock20;2; loc (39)
By using (29) again, we obtain
nX +1 +1 2 nX (n + 2 ? p)G(n + 1) 21t k i?1 ? iloc?1 k20;2; loc + k2vk21 k i ? iloc k20;2; loc i=p
i=p
But since k i ? iloc k0;2; loc is a decreasing function, we nally have (n + 2 ? p)G(n + 1) 2 (n + 2 ? p)( 21t + k2vk21 )k p?1 ? ploc?1 k20;2; loc ;
(40)
Inria
Convergence of Domain Decomposition Algorithms
15
that is 2 G(n + 1) ( 21t + k2vk21 )k p?1 ? ploc?1 k20;2; loc ; 8p n: Step 3: Global L2 estimate of n ? sta
(41)
We now de ne n+1 by
8 n+1 < ? sta in n loc; n+1 = : n+1 ? sta in loc; loc
(42)
with sta the solution of the stationary problem (2). By construction, n+1 satis es the following equations :
8 n+1 n ? + div(v n+1 ) ? n+1 = 0 in [ ( n ); > > loc loc > < t n+1 = 0 on @ ; > > > : n+1 continuous across ?i:
(43)
By multiplying this equation by n+1 and integrating by parts over loc and
n loc , we get Z n+1 ? n Z n+1 + jr n+1 j2
Z t1 @ n+1 ) n +1 2 (44) + ( 2 ( ) v:n ? n+1 @n Z@ loc 1 + ( ( n+1 )2(v:n) ? n+1 @ n+1 ) = 0: @n @ ( n loc ) 2 By taking into account the boundary conditions in (43) we obtain the following relation : Z n+1 ? n Z n+1 + jr n+1 j2
(45) Z @t n +1 n +1 n +1 ? @n ( loc ? ) = 0: ?i
RR n 2435
16
Patrick Le Tallec , Moulay D. Tidriri
But on loc ; nloc+1 ? n+1 satis es the equation ( nloc+1 ? n+1 ) ? ( nloc ? n ) + div[v ( nloc+1 ? n+1 )] ? ( nloc+1 ? n+1 ) = 0: (46) t Therefore, after multiplication by n+1 , integration by parts and from CauchySchwarz, we obtain
j
Z @ ( nloc+1 ? n+1 ) n+1 j 12 kxn k20;2; loc + 1 c21 k n+1 k20;2; loc @n 2c1 2 ?i 2 + 1 c2k n+1 k2 + kv k1 j n+1 ? n+1 j2 2c21
loc
1;2; loc
2
1
+ j nloc+1 ? n+1 j21;2; loc + j n+1 j21;2; loc ; 2 2 with c1 > 0 arbitrary. This latter inequality combined with (45) yields
0;2; loc
(47)
Z n+1 ? n n+1 + 2 j n+1 j21;2; ? c21k n+1 k20;2; loc t
21c2 kxnk20;2; loc + ( k2vck21 + 2 )j nloc+1 ? n+1 j21;2; loc : 2
1
1
By using Cauchy-Schwarz again we deduce : 1 n 2 1 n 2 1 ? c2 )k n+1 k2 + j n+1 j2 ( 2 0;2; 2 1;2; 2t k k0;2; + 2c2 kx k0;2; loc + t 1 1 2 k v k ( 2c21 + 2 )j nloc+1 ? n+1 j21;2; loc 1 (48) By using now Poincare inequality, we have (1 + (c ? 2c21)t)k n+1 k20;2; k n k20;2; + c2t kxn k20;2; loc 1
2 +t( kv k1 + )j n+1 ? n+1 j2
c21
loc
1;2; loc :
(49)
Inria
Convergence of Domain Decomposition Algorithms
17
This writes
k n+1k20;2; Ak nk20;2; + B1kxnk20;2; loc +B2 j nloc+1 ? n+1 j21;2; loc
with
(50)
A = 1 + (c ?1 2c2 )t ; B1 = c2t A;
1
1 v 21 c21
B2 = t( k k + )A:
Step 4: Final result
Let us choose c1 such that c ? 2c21 > 0. By induction, it follows ?1 Ai (B kxn?i k2 k n+1 k20;2; Apk n+1?p k20;2; + Ppi=0 1 0;2; loc +
B2j nloc+1?i ? n+1?i j21;2; loc ): Since A < 1 by assumption on c1, this implies
n X k n+1k20;2; Apk n+1?pk20;2; + A(B1 kxik20;2; loc + i=n+1?p n X j iloc+1 ? i+1j21;2; loc ): B2
(51)
(52)
i=n+1?p
Now, by using (37) and (32), we obtain
k n+1k20;2; Apk n+1?p k20;2; + A(B1 t (G(n + 1 ? p) ? G(n + 1))+ B2 21t k nloc+1?p ? n+1?p k20;2; loc ):
The same relation written between 0 and n + 1 ? p yields
RR n 2435
(53)
18
Patrick Le Tallec , Moulay D. Tidriri
k n+1?p k20;2; An+1?p k ok20;2; loc + A(B1 t (G(0) ? G(n ? p + 1)) +B2 21t k 0 ? 0loc k20;2; loc ):
By combining this relation with (53), we nally obtain k n+1 k20;2; An+1 k 0k20;2;
+ Ap+1 (B1 t G(0) + B2 21t k 0 ? 0loc k20;2; loc )
+ A(B1 t G(n + 1 ? p) + 2B2 t k nloc+1?p ? n+1?p k20;2; loc ):
By taking p such that n = 2p + q and by using (41) we conclude that
k n+1 k20;2; An+1C2 + Ap+1C3 + C4k ploc ? pk20;2; loc which implies from (31) the linear convergence of k n+1 k20;2; towards 0.
On the other hand by using (37) in (48) we obtain 1 ? c2 )k n+1 k2 + j n+1 j2 1 k n k2 ( 2 0;2; 2 1;2;
0;2;
t 1 2t 2 k v k + 2 (G(n) ? G(n + 1)) + ( 21 + 2 )j nloc+1 ? n+1 j21;2; loc : 2c1t 2c1 Therefore by using (29) we obtain 1 ? c2)k n+1 k2 + j n+1 j2 1 k n k2 ( 2 0;2; 2 1;2;
t 1 2t n 0;2; n 2 2 k ? k k v k + 2c2t (G(n) ? G(n + 1)) + ( 2c21 + 2 )( loc 2 t0;2; loc ): 1
(54)
(55)
(56)
1
Our result follows then from (41), (31) and the linear convergence of k n k0;2; .
6 Theoretical analysis of the explicit method
6.1 Local estimates
In this section, we will try to establish local estimates and a maximum principle for an arbitrary elliptic operator of second order. These tools are then needed in the convergence analysis of the explicit version of our coupling method.
Inria
19
Convergence of Domain Decomposition Algorithms
Let L be an operator written under the form
Lu = aij (x)Dij u + bi(x)Diu + c(x)u; for any u in W 2;n ( ), with a bounded domain of IRn . The coecients aij ; bi and c; i; j = 1; :::; n are de ned on . As usual, the repeated indices indicate a summation from 1 to n. We suppose that the operator L is strictly elliptic in in the sense that the matrix A of coecients [aij ] is strictly positive everywhere in . Let and denote respectively the smallest and the biggest eigenvalue of A. Let D denote the determinant of the matrix A and D = D1=n . We have 0 < D : We suppose in addition that the coecients aij ; bi and c are bounded in , and that there exists two positive real numbers and such that : = ;
(L is uniformly elliptic)
(57)
(jbj=)2 : (58) Now, we are in a position to state the principal result of this section, proved in annex. Theorem 6.1 Let u 2 W 2;n( ) and suppose that Lu f with f 2 Ln( ) and c 0. Then for all spheres B = B2R (y) of center y and radius 2R included in
and for all p > 0, we have : Z sup u ctef( 1 (u+)p)1=p + R kf k g; (59) B R (y )
jBj
B
n;B
where the constant cte depends on (n; ; R2; p), but is independent of c. Above u+ = max(u; 0).
Remark 6.1 The statement of the same theorem can be found in [12], under the assumption
jcj= :
(60) So, there the constant cte depends indirectly on c through . That is exactly what we want to avoid, since we need this constant to be independent of c, that is of t.
RR n 2435
20
Patrick Le Tallec , Moulay D. Tidriri
6.2 Technical results
In this section, we will try to apply the \maximum principle" of the previous section to study the model subproblems used in our algorithm (7)-(8).
6.2.1 Passing from the local to the global solution Let W be the sub-space of H 1( ) de ned by
W = fw 2 H 1( ); w = 0 on ?1 g: We then de ne the bilinear forms on W a(v; w) =
Z
rvrw +
(v; w) =
Z
Z
div(V v )w;
vw:
(61) (62) (63)
In the above equation, V is the velocity eld de ned by the relation (1), i.e., V satis es : 8 < divV = 0 in ; (64) : V:n = 0 on ?b : The rst basic problem associated to (8), can be written as follows : Find v 2 W satisfying :
a(v; w) + (1= )(v; w) =
Z
?b
gwd?; 8w 2 W:
(65)
Here the function g is given in H ?1=2(?b ) and the coecient is strictly positive. Moreover, from now on, we suppose that the coecients and satisfy the following relation : 1: (66) This hypothesis is not necessary but simpli es the proofs to come. Moreover, it is not restrictive, since we want the convergence for small corresponding to t small. Let d denotes the overlapping distance as in gure 2. Let then be a real number such that p 0 < < 3 =d;
Inria
21
Convergence of Domain Decomposition Algorithms
and set
p
k = =( ):
We then have the following result :
Theorem 6.2 Let v the solution of problem (65). If is suciently small, we have
q 1=2 1 C1 C2=d 1 + kV k1 d=C2 q
kvk1=2;?i
(1= )exp(?kd2=36)kg k?1=2;?b ;
(67)
where C1 and C2 are constants, with C1 depending only on d; ; and , but not on .
The proof of this theorem will be done in ve steps and is based on the local estimates described in Theorem 6.1 and on the maximum principle.
Step 1. Global estimates
By using relation (64), we have the equality :
Z
Z
vdiv(V v) = 1=2 div(V v 2)
Z
= 1=2 V:nv 2 ?
= 0; 8 v 2 W: By taking w = v in (65), we then obtain
Z
f jrvj
2 + (1= )v 2
g=
Z ?b
gv:
(68)
From this equality we deduce the following bound :
kvk21;2; kgk?1=2;?b kvk1=2;?b : Now by using the trace theorem, we obtain :
kvk1;2; (co= )kgk?1=2;?b; RR n 2435
(69)
22
Patrick Le Tallec , Moulay D. Tidriri
which implies in particular
kvk0;2; (co= )kgk?1=2;?b: Step 2. Nearby local estimates
(70)
Let L be the second order operator associate to problem (65) :
Lv = ? v + V:rv + (1= )v:
The operator
(71)
L0 = ?L
satis es the assumptions of Theorem 6.1, with c = ?1= ; f = 0. Applying then Theorem 6.1 for p = 2; y 2 i (the subdomain of width 3d externally delimited by ?i ) ( g. 2) and Ky = Bd=3 (y ), the sphere of center y and radius d=3 we obtain
kvk1;Ky c1kvk0;2;B2d=3(y): Therefore
kvk1;Ky c1kvk0;2; ;
where c1 is a constant depending only on ; ; d2 and (3=2d)n=2 . Moreover, there exist y1 ; : : :; yl belonging to i such that
(72)
2i = [y2 i B d6 (y ) [lj =1 Kyj = K with
Kyj = Bd=3 (yj ): By applying the relation (72) to each Kyj we obtain :
kvk1;K sup c1j kvk0;2; : j =1;:::;l
By setting c1 = sup c1j , we nally get j =1;:::;l
kvk1;K c1kvk0;2; : Step 3. Use of maximum principle
(73)
For any Mi in i ; we introduce (see g. (2)) :
Inria
Convergence of Domain Decomposition Algorithms
23
Bi = the ball centered on Mi of radius d=6, vi = exp[k(r2 ? d2=36)]kvk1;@Bi: The operator L applied to vi , can be written in polar coordinates (with r =
Mi M ) :
Therefore We set then : with
Lvi = 4(?k2 r2 ? k + k2 V:err + 41 )vi: Lvi 4(?k2r2 ? k2 jV:erjr + ( 41 ? k ))vi:
(74)
'(r; k) = a(k)r2 + b(k)r + c(k);
(75)
a(k) = ?k2 b(k) = ? k jV:e j
(76)
r
(77)
c(k) = 41 ? k:
(78)
2
We seek to satisfy the following relation :
0 inf '(r; k) for 0 r d6 :
As '(r; k) decreases on R+ , this will be satis ed if and only if '(d=6) 0; i.e. if and only if 2 2 kV k + 1 ? k 0: ? k 36d ? kd12 4
We replace k by its value. We therefore have to satisfy 2 d2 ? dkp V k 1 ? (36 ) 12 + 4 ? p 0: RR n 2435
24
Patrick Le Tallec , Moulay D. Tidriri
p
Multiplying by , it remains 1 (1 ? 2 d2 ) (1 + dkV k ): p 4 9 12 p From the constraint < 3 =d, we nally obtain after division
kV k ][1 ? 2d2 ]?1: '(r; k) 0 i 4p1 [1 + d12 9
(79)
p
From the relation (74) and the previous calculation, we deduce that for < 3 =d and satisfying the above inequality, we have In addition, by construction
Lvi 0 = Lv: vi v on @Bi :
Consequently, by using the maximum principle we end up at the following relation :
v vi on Bi :
In particular
v(Mi) exp(?kd2=36)kvk1;@Bi :
We do the same for ?v , and nally we have
jv(Mi)j exp(?kd2=36)kvk1;@Bi ; 8Mi 2 i: Step 4 : H 1 estimate Let 2 H 1( ) be such that : 8 < = 1 in 1 = ? loc ; : supp i [ 1 : We have from (65),
Z
(? v + div(V v ) + v= ) 2v = 0:
(80)
(81)
By using the Green's formula we deduce :
Inria
25
Convergence of Domain Decomposition Algorithms
Z
Z Z ? v 2 v = jr(v)j2 ? jrj2v 2:
On the other hand, we have :
Z
Z
div(V v ) 2v =
div(V 2 v 2=2) ?
Z
With the relations (82)-(83), (81) becomes
Z
0 =
Z
=
(83)
Z
Z
=
V:rv 2:
(82)
( jr(v )j2 + div(V 2v 2 =2) + 2 v 2= ) ? (v 2jr j2 + v 2V r )
Z
(jr(v)j2 + jvj2) + (1= ? ) 2 v 2 ? (v 2 jrj2 + v 2 V:r)
1
?
Z
(jrvj2 + jvj2) +
Z
i
Z
i
Z
(jr(v)j2 + jvj2) + (1= ? ) 2v2
(v 2jr j2 + v 2V:r ):
Hence, we obtain :
kvk
2 1;2; 1 +
Z
i
(jr
j
(v ) 2 +
jvj Z
i
2) +
Z
(1= ? ) 2v 2 =
(v 2jr j2 + v 2V r ):
With the relation (66), it follows
kvk21;2; 1[ i
Z
i
(v 2 jr j2 + v 2V:r )
kvk21; i
Z
i
( jr j2 + V:r )
(84) (85)
kvk21; i ( jj21;2; i + kk0;2; i kV kjj1;2; i ) kvk21; i jj21;2; i ( + kk0;2; i kV k=jj1;2; i ): RR n 2435
(86)
26
Patrick Le Tallec , Moulay D. Tidriri
If we take such that
jj0;2; i 1; jj21;2; i = c2=d;
where c2 is a constant, (86) then becomes
q q 1=2 kvk1;2; 1 kvk1; i c2=d 1 + kV k d=c2 : Step 5: Getting the nal estimate.
(87)
We have
kvk1;@Bi kvk1;K ;
(88)
kvk1;@Bi c1kvk0;2; :
(89)
since @Bi K . From (73), this implies
By using the relations (80) and (89) it follows :
jv(Mi)j exp(?kd2=36)c1kvk0;2; ; 8Mi 2 i:
Consequently we have
kvk1; i exp(?kd2=36)c1kvk0;2; :
(90)
kvk1; i c1co exp(?kd2=36)kgk?1=2;?b :
(91)
By using now the relation (70), we obtain :
Now, by using the relation (87) we obtain :
q 1=2 q kvk1;2; 1 coc1 c2=d 1 + kV k d=c2
(92)
(1= )exp(?kd2=36)kg k?1=2;?b:
Inria
27
Convergence of Domain Decomposition Algorithms
To conclude we use the trace theorem which yields
kvk1=2;?i c3kvk1;2; 1 : Consequently, we arrive at the following relation :
q q kvk1=2;?i coc1c3 c2=d(1 + kV k d=c2)1=2
(93)
(1= )exp(?kd2=36)kg k?1=2;?b : This is our theorem with C1 = coc1 c3 and C2 = c2.
6.2.2 Passing from the global to the local solution
Let W denote now the sub-space of H 1 ( loc ) de ned by W = fw 2 H 1( loc )=w = 0 on ?b g: We then set for any v and w in W
a(v; w) =
Z
loc
rv:rw + Z
(v; w) =
loc
Z
loc
div(V v )w;
vw;
(94) (95)
where V is the velocity eld de ned by the relation (64). The second basic problem associated to (7) can be written as : Find v 2 W such that
a(v; w) + (1= )(v; w) =
Z
?i
vj?i = h;
@v w; 8w 2 W; @n
where h is given in H 1=2(?i ). We rst have the following lemma ; Lemma 6.1 For suciently small, we have a(w; w) + (1= )(w; w) (=2)kwk21;2; loc ; 8w 2 W:
RR n 2435
(96) (97)
28
Patrick Le Tallec , Moulay D. Tidriri
Γi
Ωi
Γ1
y Ky
d/3
d/3
Γb
d/3
d/3
Ω 2ι
Mi Bi
δΒι
Figure 2: Description of the Domain loc and of the splitting used in the majoration of the local solution.
Inria
29
Convergence of Domain Decomposition Algorithms
Proof : Under the hypothesis 1= =2 + (1=2 )kV k21 , and using the Cauchy-Schwarz inequality, we obtain :
a(v; v) + (1= )(v; v) =
Z
loc
rv:rv +
Z
loc
V:rvv + (1= )
Z
loc
v2
krvk20;2 + (1= )kvk20;2 ? kV k1krvk0;2kvk0;2 krvk20;2 + (1= )kvk20;2 ? (=2)krvk20;2 ?(1=2)kV k21kvk20;2 (=2)(krvk20;2 + kvk20;2): We will also make the simplifying assumption (66). We have then the following result : Theorem 6.3 Under the notation of Theorem 6.2 and for suciently small, the solution of problem (96) is bounded by
k@v=@nk?1=2;?b
q
2 C1 C2=d 1 + 1 + kV2 k1
!
q 1=2 k V k 1 d=C 1+
2
(1 + 1= 2)exp(?kd2=36)khk1=2;?i : Here C1 and C2 are constants with C1 depending on d; v; and .
Proof
To prove this theorem, we will proceed as in the proof of Theorem 6.2.
Step 1 : Global estimates. RR n 2435
(98)
30
Patrick Le Tallec , Moulay D. Tidriri
By taking w = v in (96), we obtain :
Z
loc
jrvj2 +
Z
loc
(div(V v )v + (1= )v 2) =
Z ?i
Using the above Lemma and setting = =2, this implies
@v h: @n
(99)
(100) kvk21;2; loc k@v=@nk?1=2;?i khk1=2;?i : We will try now to get estimates on k@v=@nk?1=2;?i . From (96) and (64) we obtain : Z @v Z 1 vw): w = ( r v r w + (1 = ) V: r vw + ?i @n
loc Therefore, for any w in W , we have
j
Z @v @n wj krvk0;2; loc krwk0;2; loc + (1= )kV k1 krvk0;2; loc kwk0;2; loc ?i
+ 1 kv k0;2; loc kwk0;2; loc
(krvk20;2; loc + (1= 2)kV k21krvk20;2; loc + (1= 2)kvk20;2; loc)1=2 (krwk20;2; loc + kwk20;2; loc + (1= 2)kwk20;2; loc )1=2
[1 + (1 + kV2 k1) ]1=2kvk1;2; loc (1 + ?2 )1=2kwk1;2; loc : 2
From the trace theorem, this yields
k@v=@nk?1=2;?i (1 + ?2)1=2 1 + 1 + kV2 k1 2
!1=2
From the relations (100)-(101), we end up with
kvk1;2; loc (1 + ?2 )1=2 1 + 1 + kV2 k1 2
kvk1;2; loc: (101)
!1=2
khk1=2;?i
(102)
Inria
ΓV K
Convergence of Domain Decomposition Algorithms
31
Γi
Ωb
d/4
d/6
Bi
d/6
d/6
d/4
Mi d/6
Γb
Ω i
Figure 3: Description of the local domain loc and of the splitting used in the majoration of the global solution. and hence in particular
kvk0;2; loc (1 + ?2 )1=2 1 + 1 + kV2 k1 2
!1=2
khk1=2;?i :
(103)
Step 2 : Nearby local estimates
Let Ky = Bd=4 (y ) be the sphere centered on y and of radius d=4, with y belonging to ?V (see g. 3). By construction, ?V is the center surface of loc and i is subdomain of width 6d centered on ?V . Then by following the same argument as in step 2 of the proof of Theorem 6.2, we obtain : kvk1;Ky c1kvk0;2; loc ; (104) where c1 is a constant depending only on d; ; and . Let us consider now y1 ; : : :; yl in i such that
RR n 2435
32
Patrick Le Tallec , Moulay D. Tidriri
2i =
[ y2 i
B d6 (y)
[l j =1
Kyj = K:
By applying the relation (104) to each Kyj , we obtain
kvk1;K sup c1j kvk0;2; loc = c1kvk0;2; loc: j =1;:::;l
(105)
Step 3: Using the maximum principle For any Mi 2 i , we introduce (see g. 3) : a ball Bi centered on Mi and of radius d=6, the function vi = exp[k(r2 ? d2=36)]kvk1;@Bi. By construction of k (see previous x), '(r; k) is positive for all r 2 [0; d=6]. Then
by following the same argument as in step 3 of the proof of Theorem 6.2, we obtain the following inequality :
jv(Mi)j exp[?kd2=36]kvk1;@Bi: Step 4 : H 1 estimates Let us consider 2 H 1 ( loc), such that
(106)
= 1 in b ; (Figure 3) supp b [ i:
By choosing w = 2 v in (96) we obtain :
Z
loc
(? v + div(V v ) + (v= )) 2v = 0:
(107)
By using the Green's formula, as in step 4 of the previous theorem, we obtain :
kvk1;2; b[ i Then, if we take such that
Z
i
(v 2jr j2 + v 2V:r ):
kk0;2; i 1 Inria
33
Convergence of Domain Decomposition Algorithms
and
jj21;2; i = c2=d;
we obtain nally as in the proof of Theorem 6.2 :
q q 1=2 kvk1;2; b[ i kvk1; i c2=d 1 + kVk1 d=c2 : Step 5 : Getting the result Since @Bi K by construction, (105) and (106) imply kvk1; i exp(?kd2=36)c1kvk0;2; loc:
(108)
(109)
Furthermore by using the relation (103), it follows :
!
2 1=2 kvk1; i 1 + 1 + kV2 k1 c1(1 + 1= 2)1=2exp(?kd2 =36)khk1=2;?i:
By using the relation (108), we then obtain :
!
2 1=2 1 + k V k kvk1;2; b[ i 1 + 2 1 q q 1=2 c1 c2=d 1 + kVk1 d=c2 (1 + 1= 2)1=2exp(?kd2 =36)khk1=2;?i :
Before concluding we will establish an estimate of the quantity
k@v=@nk?1=2;?b : From (96), and by choosing
w 2 H 1 ( loc); with w = 0 on @ b \ @ i ; we obtain :
Z
b
(? v + div(V v ) + v= )w = 0:
After application of the Green's formula and from (64), we obtain :
RR n 2435
(110)
(111)
34
Patrick Le Tallec , Moulay D. Tidriri
Z @v Z w = (rv rw + (1= )V:rvw + 1 vw): @n ?b
b
By proceeding as in step 1 of the current proof, we obtain the following estimate :
! 2 1=2 1 + k V k k@v=@nk?1=2;?b (1 + 1= 2)1=2 1 + 2 1 kvk1;2; b:
(112)
The relations (111)-(112) nally yield the desired result.
6.3 Convergence of the explicit algorithm with t = 1 and small. Consider the following elliptic problem :
8 > > > < + V:r ? = 0 in ; = 1 on ?1 ; > > > : = 0 on ?b ;
(113)
that we want to solve by the fundamental algorithm (7)-(8) with "t = +1": The proposed algorithm can be written in this case as
set oloc = ol and o = o. then, for n 0; nloc and n being known, solve
8 n+1 > n+1 n+1 loc > > < + V:r loc ? loc = 0 in loc; nloc+1 = n on ?i ; > > > : nloc+1 = 0 on ?b ;
8 n+1 > + V:r n+1 ? n+1 = 0 in ; > > > < n+1 = 1 on ?1 ; > > n+1 > @ n+1 > : @ = loc on ?b : @n
(114)
(115)
@n
Inria
35
Convergence of Domain Decomposition Algorithms
We will show in this section that the algorithm described above converges. More precisely we have the following theorem.
Theorem 6.4 For suciently small, being the solution of the stationary pro-
blem (113), we have :
n+1 @ in H ?1=2(? ), i) @ @nloc converges towards @n b
ii) n+1 converges towards in H 1=2(?i ), iii) n+1 converges towards in H 1( ), iv) nloc+1 converges towards in H 1( loc). Proof:
In the rest of this section, we will take = . By the transformation n+1 ! reduce the problem to the case 1 = 0. After multiplication and integration by parts in (114), we obtain that nloc+1 = n on ?i and satis es :
n+1 ? with the solution of the stationary problem, we will
nloc+1 w + Z V:r n+1w + Z r n+1 rw loc loc
loc
loc
loc
Z
=
Z @ n+1 loc @n w; 8w 2 W: ?i
(116)
We can apply then Theorem 6.3 and obtain
q0 1 n+1 @ 0 2 loc k @n k?1=2;?b c1 c2=d 1 + 2 (1 + kV k1)
q 0 1=2 1 1 + kV k d=c
1
2
(1 + 1= 2)exp(?kd2 =36)k n k1=2;?i :
RR n 2435
(117)
36
Patrick Le Tallec , Moulay D. Tidriri
On the other hand, after multiplication by w and integration by parts in (115), we obtain the equality
Z n+1 Z Z Z @ n+1 n+1w + r n+1 rw = loc w + V: r @n w
?b
(118)
with w 2 H 1( ) and w = 0 on ?1 . By applying Theorem 6.2, we obtain :
k n+1 k
1=2;?i
q 1=2 1 c1 c2=d 1 + kV k1 d=c2 q
exp(?kd2=36)k@ nloc+1 =@nk?1=2;?b :
(119)
By using (117) and (119), we then get :
q k@ nloc+1=@nk?1=2;?b c01 c02=d(1 + 12 1 + kV k21) 1 q 1=2 1 + kV k1 d=c02 q 1=2 q 1 c1 c2 =d 1 + kV k1 d=c2 ! 2 d (1 + 1= 2 )exp ?k 18 k@ nloc =@nk?1=2;?b ; . Therefore for suciently small, the coecient of reduction will be with k = p dominated by the exponential term and will be then strictly less than 1, implying the linear convergence towards zero of
k@ nloc+1=@nk?1=2;?b : This is exactly assertion (i). By using (119), it also follows that n+1 tends to 0 in H 1=2(?i ). From (69) applied with g = @ nloc+1 =@n, we have in addition
k n+1 k1;2; co k@ nloc+1=@nk?1=2;?b ;
Inria
Convergence of Domain Decomposition Algorithms
37
and therefore k n+1 k1;2; converges to zero at the speed of k@ nloc+1 =@nk?1=2;?b . From (102) we also have
k n+1 k
(1 + 1= 2)1=2
1=2 1 2 1 + 2 (1 + kV k1 ) k n k1=2;?i
2 loc 1;2; loc and then k n+1 k1;2; also converges to zero at the speed of k n k1=2;?i .
6.4 Convergence of a xed point method for the implicit scheme.
The implicit scheme of section (5) couples the global and the local problem. To uncouple them, it is advisable to use the xed point algorithm below :
set oloc;0 = ol and o = o, being known, then for k 0; nkj+1 ?i solve
8 n+1 ? n > loc + div(v n+1 ) ? n+1 loc;k+1 > loc;k+1 = 0 in loc ; loc;k+1 > < t n+1 on ? ; (120) +1 nloc;k i > +1 = k > > : +1 nloc;k +1 = 0 on ?b ; 8 n+1 ? n > k+1 n+1 ) ? n+1 = 0 in ; > + div( v k+1 k+1 > < t +1 = (121) nk+1 1 on ?1 ; > > > : +1 =@n = @ n+1 =@n on ? : @ nk+1 b loc;k+1 We will study now the algorithm (120)-(121). By setting
RR n 2435
n+1 n+1 loc;k;q = loc;k+1 ? loc;q+1 ;
(122)
n+1 k;q = ( k ? nq +1 );
(123)
38
Patrick Le Tallec , Moulay D. Tidriri
we have that loc;k;q and k;q verify the following equations :
8 > > < > > :
loc;k;q =t + div (v loc;k;q) ? loc;k;q
8 > > > < > > > :
= 0 in loc ; loc;k;q = k?1;q?1 on ?i ; loc;k;q = 0 on ?b ;
(124)
k;q =t + div (v k;q) ? k;q
= 0 in ; k;q = 0 on ?1 ;
loc;k;q on ? : @@nk;q = @ @n b
(125)
If t is suciently small, we have from the previous section that k;q and loc;k;q +1 are Cauchy sequences converge linearly to zero. Hence the sequences nk +1 and nloc;k which converge linearly towards the unique solutions n+1 and nloc+1 of the implicit scheme. This guarantees the convergence of the above xed point algorithm.
Inria
Convergence of Domain Decomposition Algorithms
39
7 Numerical analysis of the stability of the algorithm (7)-(8) In this section we focus our attention on the numerical solution of the steady problem (2), using now the explicit time marching algorithm (7)-(8) studied in the previous paragraphs. We rst assume here that the boundary condition on ?b in (8) is explicit n+1
n
( @ @n = @ @nloc ) so that the resulting algorithm is parallel (Jacobi type). Here, denotes the domain surrounding the obstacle (an ellipse in our numerical study) as described in Figure 1. The global and local domains are discretized by fully overlapping compatible nite element grids. The global mesh contains 1378 nodes and 2662 elements (see gure 4). Further the time marching algorithm is being initialized by setting 0 to zero.
RR n 2435
40
Patrick Le Tallec , Moulay D. Tidriri
Figure 4: Description of the nite element mesh and of the local subdomain.
Inria
Convergence of Domain Decomposition Algorithms
41
In a rst step, the velocity eld is obtained by solving the following problem : div v = curl v = 0 (inviscid incompressible ow), v1 = (1; 0); v:n = 0 on the body?b; with a rst order nite element method using the same global mesh.
RR n 2435
42
Patrick Le Tallec , Moulay D. Tidriri
If we set v = 0, the algorithm may or may not converge depending on the values of t. More precisely, we observe that the algorithm converges linearly when t < 0 and is ndivergent otherwise. This is graphically shown in gures (4-7) where +1 ? n k k the values of tk 0 k1 are plotted versus the iteration count n for t taken on the values 10?6, 10?1, 1, and 10. Further, when the velocity is taken suciently large, the algorithm becomes unconditionnally stable. In particular, initializing our algorithm by 0 = 0 with kv1 k = 1, = 0:1 and t = 100 results in a converging scheme ( g. 8). By intuition such a behavior seems natural. An overestimation of the solution n at the interface ?i implies an overestimation of the friction forces on ?b . For suciently small time steps, this overestimation will not a ect the value of n+1 on ?i and can therefore be ignored at the next time step. If the Reynolds is suciently large, this error will only a ect the wake region but will not have any in uence at the interface ?i . To the contrary, for large t and small , the introduced error does a ect the value of n+1 on ?i . The in uence of the error on n+1 may be ampli ed throughout the iteration process.
Inria
Convergence of Domain Decomposition Algorithms
43
Another variant of the algorithm consists of replacing the explicit Dirichlet condition
nloc+1 = n on ?i in the algorithm (7)-(8)
by the following semi-implicit condition
nloc+1 = n+1 on ?i:
In fact, this implies replacing the previously parallel algorithm (Jacobi like ) by the sequential algorithm (Gauss-Seidel like ). When we solve the pure di usion problem (i.e. with ow velocity v = 0) with = 1 and t = 1 (respectively t = 2) we obtain a better convergence history : the speed of the new algorithm is linear and clearly faster than the parallel algorithm. the domain of convergence is moderately larger (see table 1). To study experimentally in more details the convergence behavior of both algorithms we assume that we have a linear behavior of our algorithm, and hence that the error at the iteration n will satisfy the following inequality
k n+1 ? nk1 K nk 1 ? 0k1:
The algorithm converges if K < 1. An estimate for K can be found by considering as in table 1 the ratio n+1 n ? n1 log k k 1 ?? 0k k1 = ?log K
1
which is displayed as a function of ( t) for n = 14 and di erent values of V = v . A negative value of this ratio means divergence of the algorithm. As expected, this ratio is positive for suciently small values of t and converges towards zero as t goes to zero. In this table, we observe that for V = 0; t < 0 2, the algorithm converges. However the convergence is slow since the minimal contraction constant Kmin (for the optimal value of t) is close to one (see table 2). For V = 10, then the algorithm appears to converge for a much larger range of values of t and the optimal contraction constant is much smaller. This is summarized on table 2 where we have displayed the best possible contraction constants for each of the coupling algorithms and di erent values of the Reynolds V = v .
RR n 2435
44
Patrick Le Tallec , Moulay D. Tidriri
t 1/1000 1/10 1/2 2 Gauss-Seidel 0.06 0.1 0.22 0.5 V =0 Jacobi 0.1 0.22 0 V =0 Gauss-Seidel 0.03 0.25 1.46 2.12 V = 10 Jacobi 0.03 0.28 1.15 1.15 V = 10 Jacobi 0.23 2.79 2.8 2.7 V = 1000
5 10 50 1000 -0.27 -0.5 -0.75 -0.8 -0.09 -0.25 -0.4 -0.41 2.8
2.6
2.4
2.4
1.15
1.15
1.14
1.14
2.75
2.8
-
-
Table 1: Contraction constant (in fact minus its logarithm) in function of t for the explicit (Jacobi) and semi-explicit (Gauss-Seidel) version of our coupling algorithm. We observe divergence for V = 0 and t > 2 and convergence otherwise.
Jacobi (parallel) V Kmin 0 0:85 10 0:50 103 0:14
Gauss-Seidel (sequential) V Kmin 0 0:68 10 0:11
Table 2: Minimal contraction constant versus the Reynolds V for both sequential and parallel versions of the algorithm.
Inria
45
Convergence of Domain Decomposition Algorithms
Figure 5: Convergence of the Time Marching Algorithm: k tk ?0 k 1 k are plotted versus the iteration count n for t = 10?6 , v = 0 (Jacobi). Observe the very slow convergence. n+1
RR n 2435
n
46
Patrick Le Tallec , Moulay D. Tidriri
Figure 6: Convergence of the Time Marching Algorithm: k tk ?0 k 1 k are plotted versus the iteration count n for t = 10?1 , v = 0 (Jacobi). n+1
n
Inria
47
Convergence of Domain Decomposition Algorithms
Figure 7: Convergence of the Time Marching Algorithm: k tk ?0 k 1 k are plotted versus the iteration count n for t = 1, v = 0 (Jacobi). n+1
RR n 2435
n
48
Patrick Le Tallec , Moulay D. Tidriri
Figure 8: Divergence of the Time Marching Algorithm: k tk ?0 k 1 k are plotted versus the iteration count n for t = 10, v = 0 (Jacobi). n+1
n
Inria
49
Convergence of Domain Decomposition Algorithms
Figure 9: Convergence of the Time Marching Algorithm: k tk ?0 k 1 k are plotted versus the iteration count n for t = 10 and the ow velocity is equal to 1 (Jacobi). n+1
RR n 2435
n
50
Patrick Le Tallec , Moulay D. Tidriri
8 Conclusion We have analysed the convergence properties of a standard time marching algorithm for solving a domain decomposed advection-di usion problem. We were able to prove theoretically the unconditional stability and linear convergence of the fully implicit algorithm (x5). Using the maximum principle, we were also able to prove in x6 the linear convergence of the semi explicit algorithm used with t = +1 on the operator Lu = ?u + div(vu) + u if was suciently large, and the linear convergence of the xed point algorithm applied at each time step to the solution of the implicit problem, independently of any discretization strategy. When using the uncoupled semi-explicit algorithm in the general case, numerical evidence indicate that this algorithm is unstable for large values of t and small overlapping, and that it becomes linearly convergent when t is below a Reynolds dependent threshold (x7). This conditional stability is not a real issue for practical CFD problems because most solvers already require to use small time steps inside each domain. Nevertheless, it would be nicer to derive an uncoupled unconditionnally stable version of the present time marching algorithm.
Inria
Convergence of Domain Decomposition Algorithms
51
References [1] R. Glowinski, G. Golub and J. Periaux (eds), Proceedings of the First International
Symposium on Domain Decomposition Methods for Partial Di erential Equations, Paris, France, January 7-9, 1987 , (SIAM, Philadelphia, 1988). [2] T. Chan, R. Glowinski, J. Periaux and O. Widlund (eds), Proceedings of the Second International Symposium on Domain Decomposition Methods for Partial Di erential Equations, Los Angeles, California, January 1988 , (SIAM, Philadelphia, 1989). [3] T. Chan, R. Glowinski (eds), Proceedings of the Third International Symposium on Domain Decomposition Methods for Partial Di erential Equations, Houston, Texas, March 20-22, 1989 , (SIAM, Philadelphia, 1990).
[4] R. Glowinski, Y. Kuznetsov, G. Meurant, J. Periaux and O. Widlund (eds),Proceedings of the fourth international symposium on Domain Decomposition Methods for Partial Di erential Equations, Moscou, June 1990 , (SIAM, Philadelphia, 1991). [5] T. Chan, D. Keyes, G. Meurant, S. Scroggs and R. Voigt (eds), Proceedings of
the fth international symposium on Domain Decomposition Methods for Partial Di erential Equations, Norfolk, May 1991 , (SIAM, Philadelphia, 1992). [6] A. Quarteroni (ed), Proceedings of the sixth international symposium on Domain Decomposition Methods for Partial Di erential Equations, Como, June 1992 ,
(AMS, Providence, 1994). [7] Y. Achdou and O. Pironneau, A fast solver for Navier-Stokes Equations in the laminar regime using mortar nite element and boundary element methods, Technical Report 93-277 (Centre de Mathematiques Appliquees, Ecole Polytechnique, Paris, 1993) [8] A. D. Aleksandrov, Majoration of solutions of second order linear equations, Vestnik Leningrad univ. 21, 5-25(1966) English translation in Amer. Math. Soc. Transl. (2) 68, 120-143(1968). [9] I. Ya. Bakel'man, Theory of quasilinear elliptic equations. Siberian Math. J. 2, 179-186(1961).
[10] J. Bramble, R. Ewing, R. Parashkevov and J. Pasciak, Domain decomposition methods for problems with partial re nement SIAM J. Sci. Stat. Comp. 13 (1992) 397{410. [11] C. Canuto and A. Russo, On the Elliptic-Hyperbolic Coupling. I: the Advection Di usion Equation via the -formulation, Math. Models and Meth. Appl. Sciences (to appear) (1993). [12] D. Gilbarg, N. S. Trudinger, Elliptic partial di erential equations of second order.Berlin-Heidelberg-New York, Springer Verlag 1977.
RR n 2435
52
Patrick Le Tallec , Moulay D. Tidriri [13] W. Gropp and D. Keyes, Domain Decomposition Methods in Computational Fluid Dynamics, Int. J. Num. Meth. Fluids 14 (1992) 147-165. [14] P. Le Tallec and M. D. Tidriri, Coupling Navier-Stokes and Boltzmann. Submitted to J. Comp. Phy. [15] L. Marini and A. Quarteroni, An iterative procedure for domain decomposition methods: a Finite Element approach, in [1]. [16] L.D. Marini, A. Quarteroni, A relaxation procedure for domain decomposition methods using Finite Elements, Numer. Math. 55, (1989) 575{598. [17] A. Quarteroni, G. Sacchi Landriani and A Valli, Coupling of Viscous and Inviscid Stokes Equations via a Domain Decomposition Method for Finite Elements, Technical report UTM89-287 (Dipartimento di Mathematica, Universita degli Studi di Trento, 1989). [18] M. D. Tidriri, Couplage d'approximations et de modeles de types di erents dans le calcul d'ecoulements externes,these, Universite de Paris IX, May 1992. [19] M. D. Tidriri, Domain Decomposition for Incompatible nonlinear models. INRIA Research Report. To appear.
[20] M. D. Tidriri, Domain decomposition for compressible Navier-Stokes equations with di erent discretizations and formulations. To appear in J. Comp. Phy.
Inria
Convergence of Domain Decomposition Algorithms
53
Appendix The main result of this section relies on the notion of a contact set. If u is a continuous arbitrary function on , the upper contact set, denoted ?+ or ?+u , is the sub-set of , de ned by ?+ = fy 2 ; 9p(y) 2 IRn such that u(x) u(y) + p (x ? y) 8x 2 g: (126) We see that u is a concave function on i ?+ = . When u 2 C 1( ) we must have p = Du(y) in the relation (126). In addition, when u 2 C 2( ), the Hessian matrix D2 u = [Dij u] is negative on ?+ . In general, ?+ is closed in . If u is a continuous arbitrary function on , we de ne the "normal mapping" (y) = u (y) at point y 2 by (y) = fp 2 IRn; u(x) u(y) + p:(x ? y) 8x 2 g: (127) We can see that (y) is non empty i y 2 ?+ . In addition when u 2 C 1( ), we have (y) = Du(y) on ?+ ; in other words is the gradient eld of u on ?+ .
As a particular case of the Bakelman-Alexandrov ([8] and [9]) maximum principle, we have under the above notation. Lemma .1 For u 2 C 2 ( ) \ C o( ), we have :
d kaij Dij u=D kn;?+
nwn1=n with d the diameter of and wn the volume of a unit sphere in Rn .
sup u sup@ u +
For further details see [12]. We now proceed to the proof of Theorem 6.1, by following the steps of [12]. We take B^ = B1 (0) and the general case will be deduced by considering the coordinate transform, x ! x^ = (x ? y)=2R. We will begin, in rst step, by showing this result for u 2 C 2( ) \ W 2;n( ) and then in a second step we will deduce the result for u 2 W 2;n( ).
Step 1:
We suppose that u 2 C 2 ( ) \ W 2;n( ). For 1, we consider the cut o function de ned by (^x) = (1 ? jx^j2) : By di erentiation, we obtain D^ i = ?2 x^i (1 ? jx^j2) ?1 ;
RR n 2435
54
Patrick Le Tallec , Moulay D. Tidriri
By setting
D^ ij = ?2 ij (1 ? jx^j2) ?1 + 4 ( ? 1)^xi x^j (1 ? jx^j2) ?2 : v = u;
we then obtain a^ij D^ ij v = a^ij D^ ij u + 2^aij D^ i D^ j u + ua^ij D^ ij (f^ ? ^biD^ i u ? c^u) + 2^aij D^ i D^ j u + ua^ij D^ ij : Let ?+ = ?+v be the upper contact set v, in the sphere B^ ; we have : u > 0 on ?+ :
If x 2 @ B^ such that p:(x ? y) < 0 we indeed have v(x) = 0. Consequently v(y) + p:(x ? y) v(x) = 0:
Moreover, using the concavity of v on ?+ , we can estimate the following quantity : ^ j = (1=)jDv ^ ? uD ^ j: jDu Indeed,
^ (1=)( Dv ^ + u D ^ ) Du
^ j) (1=)( 1 ?v jx^j + ujD 2(1 + )?1= u:
In that way, we have on ?+ the following inequality : ?a^ij D^ ij v f(16 2 + 2 )^ ?2= + 2 j^bj?1= + c^gv + jf^j: Since c^ 0, we deduce the inequality ?a^ij D^ ij v f(16 2 + 2 )^ ?2= +
(128)
2 j^bj?1= gv + jf^j
c1 ^?2= v + jf^j; with c1 = c(n; ; ; ^) independent of c^.
Inria
Convergence of Domain Decomposition Algorithms
55
Consequently, by applying Lemma 6.1 on B^ , we obtain, for 2 : ^ sup v ( d1=n )( ^1 )kc1 ^ ?2= v + jf^jkn;B^ : nwn D B^ By using the relation (6.1), it comes ^ ^ sup v ( d1=n )c1 k?2= vkn;B^ + ( d1=n )( ^1 )kf^kn;B^ nwn nwn B^ c1 d^(k?2= vkn;B^ + (1=^)kf^kn;B^ )
c1 d^(k?2= v+ kn;B^ + (1=^)kf^kn;B^ ) c1 d^((supv+ )1?2= k(u+ )2= kn;B^ + (1=^)kf^kn;B^ );
where c1 is a constant depending only on n; ; and ^. Here, d^ is the diameter of B^ (d^ = 2). By using the Young inequality under the form ab "aq + "?r=q br
for q = (1 ? 2= )?1 and r = =2, we have
2 (supv+ )1?2= k(u+ )2= kn;B^ "supv+ + "1? =2 k(u+ )2= k = n;B^ ; 8" > 0:
By taking " = 1 ^ and plugging in our inequality on v, we obtain : 2c1 d 2 sup v (1=2) sup v+ + (1=2)1? =2 (c1 d^) =2 k(u+ )2= k = n;B^
B^
(129)
^ ^)kf^kn;B^ : +(c1 d= We want to prove the theorem for all p > 0. We will treat separately the cases p n and p > n. If p n, we set = 2n=p. In this case we have 2 + k(u+ )2= k = n;B^ = k(u )kp;B^ :
Plugging this in our inequality on v, we obtain : ^ ^)kf^kn;B^ : (1=2) sup v (1=2)1? =2 (c1 d^) =2 k(u+ )kp;B^ + (c1 d= B^
Consequently, we obtain the following inequality ;
RR n 2435
56
Patrick Le Tallec , Moulay D. Tidriri
Z
^ 2^)kf^kn;B^ g: sup v c2 f( ^ (u+ )p )1=p + (d= B^
B
On the sphere B1=2 (0), the cut o function satis es 1= (1=2) : It follows, then sup u
B1=2 (0)
sup (v= )
B1=2 (0)
2 sup v: B^
Finally we end up at the desired estimate
Z
^ 2^)kf^kn;B^ g: sup u c3 f( (u+ )p )1=p + (d= ^
B1=2 (0)
B
for u in W 2;n( ) \ C 2( ). The constant c3 above depend only on n; ; and ^, but is independent of c^. On the other hand if p > n, we have : 2n= < p; 8 2: Then, it follows (by assuming 2) jB^ j?1=(2n= )k(u+ )k2n= ;B^ jB^ j?1=pku+ kp;B^ : But 2 ku+ k2n= = k(u+ )2= k = n;B^ ;
and therefore, by processing as before, we obtain the desired estimate
Z
^ 2^)kf^kn;B^ g sup u c4 f( (u+ )p )1=p + (d=
B1=2 (0)
B^
for u in W 2;n( ) \ C 2( ). The constant c4 above depends only on n; ; and ^, but is independent of c^. Transformation x^ ! x. By construction, D^ ij = R?2Dij , thus ^ = R?2 and ^ = R2. In addition, we have jB j = wn (2R)n and jgjp;B^ = R?n=pjgjp;B .
Inria
Convergence of Domain Decomposition Algorithms
57
Written in term of x, the last inequality becomes n
sup u c4 f( 2jBwjn
BR(y)
Z
B
1=n (u+ )p dx)1=p + ( 2wn R )kf kn;B g;
with c4 a function of n; ; ^ = R2 and p. This is the desired estimate for u 2 W 2;n( ) \ C o ( ).
Step 2:
Now, let u 2 W 2;n( ). By density, let (um ) be a sequence of functions of C 2(B ), converging towards u in W 2;n(B ). The injection of W 2;n(B ) in C o (B ) is continuous, consequently (um ) converges uniformly towards u in B . We have Lum = L(um ? u) + Lu
f + L(um ? u): By setting, fm = L(um ? u), we observe by construction that fm converges towards 0 in Ln ( ). As um 2 W 2;n( ) \ C 2( ) and f~m = f + fm is in Ln ( ), the estimate (59) is valid also for um , so that we have
Z sup um ctef( jB1 j (u+m )p )1=p + R kf~kn;B g: B (y) B R
(130)
Using previous results and taking the limit, we have : Z sup u ctef( jB1 j (u+ )p )1=p + R kf kn;B g: BR (y)
B
Observe also that by replacing u by ?u, the theorem can be extended easily to the case of supersolutions and solutions of the equation : Lu = f:
RR n 2435
Unité de recherche Inria Lorraine, Technopôle de Nancy-Brabois, Campus scientifique, 615 rue du Jardin Botanique, BP 101, 54600 Villers Lès Nancy Unité de recherche Inria Rennes, Irisa, Campus universitaire de Beaulieu, 35042 Rennes Cedex Unité de recherche Inria Rhône-Alpes, 46 avenue Félix Viallet, 38031 Grenoble Cedex 1 Unité de recherche Inria Rocquencourt, Domaine de Voluceau, Rocquencourt, BP 105, 78153 Le Chesnay Cedex Unité de recherche Inria Sophia-Antipolis, 2004 route des Lucioles, BP 93, 06902 Sophia-Antipolis Cedex
Éditeur Inria, Domaine de Voluceau, Rocquencourt, BP 105, 78153 Le Chesnay Cedex (France) ISSN 0249-6399
