- Email: [email protected]
Artificial damping techniques for scalar waves in the frequency domain
Computers Math. Applic. Vol. 31, No. 8, pp. 1-12, 1996
Pergamon
Copyright(g)1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved S0898-1221(96)00025-9 0898-1221/96 $15.00 + 0.00
Artificial D a m p i n g Techniques for Scalar Waves in the Frequency D o m a i n S. KIM Department of Computational and Applied Mathematics, Rice University P.O. Box 1892, Houston, TX 77251-1892, U.S.A. skim~caam, rice. edu
M. LEE Center for Applied Mathematics, Purdue University W. Lafayette, IN 47907, U.S.A. miyoung~math, purdue, edu (Received October 1995; accepted November 1995)
A b s t r a c t - - A n artificial damping algorithm for solving the Helmholtz problem is considered. When the imaginary part of the wave number is small, the problem is known to be difficult to solve. In this paper, we propose an efficient artificial damping algorithm which can be viewed as a rational iteration. Each damped problem is solved incompletely by a nonoverlapping domain decomposition method. Convergence arguments are presented• Numerical results run on an nCUBE2 are presented to demonstrate the effectiveness of the method• K e y w o r d s - - A r t i f i c i a l damping, Wave propagation, Robin boundary condition, Domain decomposition method•
1. I N T R O D U C T I O N Let Ft = (0, 1) d, d = 2 or 3, with the b o u n d a r y F = 0Ft. Consider the following (complex-valued) scalar wave problem: 022 xEft, - A u - c~x)2 U -t- i q ( x ) 2 u = f ( x ) , cgu
+
w
c(x)
(1) u--0,
xEF,
where i is the imaginary unit, w > 0 is the angular frequency, c denotes the wave speed in the m e d i u m where the wave is propagating, q2 __ bw, where b > 0 is the friction coefficient, and u is the unit o u t w a r d normal to F. Here the coefficients are sufficiently regular so t h a t the existence and uniqueness of a solution of (1) lying in HS(ft) for some s > 1 for reasonable f are assured. T h e second equation of (1) represents a first-order outgoing radiation condition t h a t allows n o r m a l l y incident waves to pass out of Ft transparently. T h e problem (1) models the p r o p a g a t i o n of t i m e - h a r m o n i c waves such as electromagnetic waves, discretizations of the t i m e - d e p e n d e n t Schr5dinger equation by implicit difference schemes, inverse scattering problems, seismic waves, and underwater acoustics. In this paper, we are mainly interested in the case q -- 0,
or q is small.
(2)
I n this case, the linear system arising in the finite difference discretization of the problem (1) is difficult t o solve. In addition to having a complex-valued solution, the problem is neither coercive Typeset by .AAdS-TEX
2
S. KIM AND M. LEE
nor Hermitian symmetric. As a consequence, most standard iterative methods either converge slowly or fail to converge. (For q = 0, there is no convergent conjugate gradient (CG)-type algorithm for the problem [1-3].) For an iterative algorithm for solving (1), a coarse grid solution can be used to obtain a better initial guess a n d / o r to accelerate the convergence of the iteration. Such an idea has been widely studied in terms of multigrid (MG) methods where one uses a collection of successive coarser finite element grids; see [4-6] and references therein. There are two steps in MG methods, coarse grid correction and smoothing. The MG method is known to be effective for solving certain coercive elliptic problems. The reason why the MG method is so efficient is, roughly speaking, the following: in the coarse grid correction, the low and medium frequency components of the error (corresponding to small and medium eigenvalues) are significantly reduced, while the smoothing step reduces the high frequency components of the error. Since one needs to choose the mesh size h such that h. maxxef~(w/c(x)) is less than 2/3 to 1 for stability reasons [7], which is the same as taking at least 6 to 9 points per wavelength (27rc/w), it is hard to solve problem (1) even on the coarsest grid for large wave numbers. (So MG methods cannot be applied to (1) efficiently without designing a solver for the coarsest grid problem.) For such a difficult problem, we will propose a rational iteration instead of considering MG methods or CG-type accelerations. An outline of this paper is as follows. In Section 2, we propose an iterative algorithm that introduces artificial damping/attenuation and that can be viewed as a rational iterative algorithm. A convergence argument for the algorithm is presented. In Section 3, we consider a nonoverlapping domain decomposition iterative procedure to solve a fractional step of the rational algorithm presented in Section 2. So the domain decomposition iteration is considered as the inner iteration, while the rational iteration is the outer iteration. In Section 4, we present the main algorithm in which the inner iteration finishes incompletely and the coefficient of artificial damping is modified in early steps of the outer iteration. The method is implemented in an nCUBE2 parallel computer and some results are given in Section 5 to demonstrate the effectiveness of the method. In this section, we present heuristic strategies for finding efficient algorithm parameters introduced in Section 4. The last section includes conclusions and possible applications.
2. A R T I F I C I A L
DAMPING
When problem (1) is discretized by a finite difference/element method of a mesh size h, it can be represented as follows: find u satisfying Au = k.
(3)
Here A (= A h) is a symmetric complex matrix and k is the force vector. It is well known that the matrix A in (3) is poorly conditioned; as a consequence, most standard iterative methods either fail to converge or converge so slowly as to be impractical. Our next objective is to present an iterative algorithm for (3) which introduces artificial damping. First, we cite the following lemma [8]. LEMMA 2.1. Let A be an eigenvalue o[ A. Then, Im(~) > 0 for w2h sut~ciently small. Now, consider the following iterative procedure: for a given initial guess u °, solve the following, m_>O: (a) r m = k - Aura; (b)
if Ilrmll < e, stop;
(c) (d)
find pm such that (A + idI) p m = rra; um+l
(4)
= U m _1_ p r o .
Here 6 is the stopping criterion for the iteration, d is the coefficient of artificial damping, I is the identity matrix, and I1" II denotes the e2-norm. In Section 5, we will consider a heuristic strategy for finding an efficient coefficient of damping.
Artificial D a m p i n g T e c h n i q u e s
3
THEOREM 2.2. The iteration (4) converges for all d > O. PROOF. F r o m (4a), (4c), and (4d), U m + l 7__ U m + p m
= u m + (A + idI) -1 (k - A u m) =
(5)
+ (A + i d I ) - l A (u - urn).
U m
Let e m = u - um. Then, from (5), one has e m+l = e m - (A + i d X ) - l A e m
(6)
-- i d ( A + idI) -1 em. Let a(A) be the s p e c t r u m of the m a t r i x A. Then, Ilem+lH
<
max
- -
. Ilem]].
AEcr(A) A i/id
It follows from L e m m a 2.1 t h a t (4) converges for all d > 0.
|
REMARK. In algorithm (4), the p r o b l e m is how one can invert A + i d I efficiently. Since A + i d I is b e t t e r conditioned, one can find easily a convergent iterative algorithm for (4c). We do not have to solve (4c) directly/completely, in practical simulations. T h e p r o b l e m can be solved incompletely using an iterative algorithm. REMARK. T h e algorithm (4) can be rewritten as follows. Multiply (4d) by A + i d I and a p p l y (4a) and (4c) to the resulting equation. Then,
(A + idI) u "~+1 = k + i d u m.
(7)
Again, equation (7) can be rewritten as a discrete version of --Aum+l -
( w2 -~-iq
2) u m + l + i ~ 2 u m + l = f ( x ) + i ~ ? 2 u ra,
x•~,
(8)
Ou m+l
O~
+ iw-um+lc
= O,
x • F,
for some U > 0. For the 5-point finite difference discretization of a uniform grid size h, we can see 772 = d/h 2. In Section 3, we will consider a domain decomposition m e t h o d for the d a m p e d p r o b l e m (8). REMARK. Let Ud be the solution of
(A + idI) Ud = k.
(9)
T h e n , from (3) and (9), u
--
U d
----
id (A + i d I ) - 1 u,
and, therefore, Iku - udll <
max
.
(10)
T h e p e r t u r b e d solution u a can be a good initial guess for u; see Figure 2.1. T h e r e finite difference solutions of (9) on the line segment Y0.5 = {(Xl, x2) • gt : x2 = 0.5} are depicted for the point source f = 5(x0 - x), x0 = (0.5,0.5), and for the uniform grid size h = 1/100. We selected the p r o b l e m coefficients such t h a t q - 0, w = 70, and c = c4 in (26) below. F r o m Figure 2.1, one can see t h a t even t h o u g h the m a g n i t u d e of Ud is different from t h a t of u, the phase of Ud seems to be identical to t h a t of u.
S. K I M
AND
M.
LEE
Re(u}:
y=0.5
0.00004 - -
: d=O
..... : d = l
0.00003
0.00002
0.00001
-0.00001 i
.
.
.
.
0
i
.
.
.
.
0.2
i
i
,
,
0.4
i
!
J
.
•
•
'
0.6
.
.
.
.
'
0.8
,=0.5
Im{u)
0.00001
~ J
-0.00001
----: -0.00002
d=0
..... : d = l
0
0.2
0.4
0.6
0.8
X
Figure 2.1. Finite difference solutions of (9) on Y0.5 for the point source f = x o ---- ( 0 . 5 , 0 . 5 ) , w h e n q --- 0, w = 70, c = c4 in ( 2 6 ) below, and h = 1 / 1 0 0 .
3. S O L U T I O N
FOR
THE
DAMPED
6(xo -x),
PROBLEM
In this section, we will consider an iterative algorithm for inverting A + idI in (4c), or equivalently, for solving (8). We will employ the domain decomposition method in [9] for solving (8). Let {flj, j = 1 , . . . , M } be a partition of ~: M
j=l
Assume that f2j, j = 1 , . . . , M, are also rectangular/cubic regions. Let
r~ = r ~ O~j,
r~k = rk3 = o~j ~ O~k.
Let us consider the decomposition of the problem (8) over {~2j }. It can be easily checked that
Artificial Damping Techniques
5
the problem (8) is equivalent to the following. Find u~ +1, j = 1 , . . . , M, such that (a)
- - A U ? + ' --
K 2 u ? +1 + i772u? +1 = f + i~2u?,
x E fly,
OU? +1
(b)
O--~j + i~um+lc3
= O,
0 U ? +1
(c)
o.j
X E Fj,
(11)
C~U~n+l
+ iZ~'[+~ -
o.k
+ iZu~+l'
z e rjk,
where K is the wave number defined as K2-
022
c2
iq 2,
(12)
vj is the unit outward normal from flj, and fl is a complex function on the subdomain interfaces with Re(fl) > 0. In (11), the consistency conditions u~+l = u~+l'
Ovj
+
cgvk
O,
on Fjk
(13)
are replaced by the Robin interface boundary condition (11c). This replacement is often more convenient and it is not difficult to show that (11c) and (13) are equivalent; see [8-12]. The domain decomposition iterative algorithm for (11) can be defined as follows. Choose an initial guess u~n+l'°, j = 1 , . . . , M, and then build u~n+l'n, n > 1, recursively by solving (a) - - A u ? + l ' n - K 2 u ? +l'n + i~2u? +l'n = f + i~2u?, 0_ u ?_+ l ' n
(b)
~- $ 03 c U jra+l n --~ O,
x C F j,
Ovj 01U, n?_+_ 0/]j
(C)
-~- ifl U ? + l ' n __
x E ftj,
GUkO rn+l ,n-1 .~ rn+l,n-- 1 0/]k + zp u k
(14)
x E Fjk,
where u ~ is the limit of sequence {ujm n' : n _ 0}. The above algorithm was analyzed in [11] in differential, rather than discrete, level. Note that the Robin interface boundary condition (14c) imposes the continuity of the solution u and its flux. Since most finite element/difference methods admit discontinuities of the flux on the element interfaces, (14c) should be modified in discretized problems. We will approximate (14) by the five-point finite difference method of a uniform grid size h. Let A h u j be the centered five-point difference approximation of A u j , O¢,juj be the centered difference for ~ on the boundary Fj, and ~f, j k U j and Ob,jkuj denote the forward and the backward finite differences for ~ on the interface Fjk, respectively. (Here we assume an exterior bordering of the subdomains.) We can define a finite difference domain decomposition iterative algorithm as follows: for given u ~ +1'°, j = 1 , . . . , M, solve ujm+l,n , n ~ 1, satisfying ^ .~+1,,~ _ K 2 u ~ + l , n + i , 2 u ~.+ l , n = f + z~ 2uj.~ , (a) --,,hUj
(b) (C)
i W u m + l ''~ = O, O¢,j u ~ +1''~ + c 3 0
f, j k
um+l,n i~u?+l,n m + l , n - 1 -- ,~ m + l n - 1 j + = --~o,kj Uk ~ $~OUk ,
x C ~j, x E F j,
(15)
x E Fjk.
Algorithm (15) can be viewed as a generalized Schwarz splitting with minimum overlap [13]. Note that the modified Robin interface boundary condition (15c) imposes the continuity of the discrete solution only. One can prove the convergence of algorithm (15) for q2 + r]2 ;> O.
6
S. KIM AND M. LEE
A convergence argument and strategies for finding efficient parameter/9 for (15) can be found in [9,14]. It is numerically verified that the convergence rate of (15) is independent of the grid size h. Some applications of the algorithm to finite element methods can be found in [9]. In the following, we will illustrate the algebraic system of algorithm (15). First, let us consider the algebraic equations of (15) corresponding to interface nodal points; see Figure 3.1. Let F = f + i y 2 u ~ and drop superscript m and m + 1, for convenience. At the point C, (15a) reads (4 - h 2 K ~ + ih2zl~) uj~,c - u jn, E
- - U jn, W
- - U j,~ , S - - U jn, N
h2Fc,
(16)
where ujn,c = u~(C), the value of u~ at the point C, and the others are similarly defined. The term ujn,E in (16) evaluated at a point out of the subdomain f~j can be substituted by using (15c). Equation (15c) is written as n
n
n-1
uJ'r' - Uj'c + i/9 n ~tj, c h
n--1
Uk,E -- Uk,c + i / 3 h
--
,~-1
Uk,C ,
or, equivalently,
(17)
uj,Bn - (1 - i/~h) u~, c = Uk,En-1 _ (1 -- i/3h) u n-ik,c. Substituting (17) into (16), we have [4 - h 2 g 2 + ih2~l~ - (1 - i/3h)] uj~,c - ujn,w - uj~,s - Uj~,N n-1 = h 2 F c + Uk,E
_
(1 -
it3h) u ~ 1.
(18)
In the same manner, one can treat cross points arising in a box-type decomposition of the domain. Clearly, algorithm (15) is applicable to nonuniform meshes with minor changes. N w
c
E
S f/j
Fjk
~;,
Figure 3.1. Five point stencil at interface nodal points.
Now, we consider the matrix representation of algorithm (15). For simplicity, let the domain be decomposed into two subdomains; see Figure 3.2. Give a nodal point ordering as follows: we count nodes in f~l the first, those on F12 next, and those in ~2 the last. Then, by dropping the superscripts m and m + 1, one can rewrite (8) as
LorA11 A120 [Ul] If1] A22 A23
u2=
A32
u3
A33J
(19)
f3
°IUl [fl
Problem (11) can be represented by the following enhanced system: All A21 A21
A12 D1 C2
0 C1 D2
A23
u2
A23
u2
0
0
A32
A33
u3
D2 + C2,
and
=
f2 f2 f3
'
(20)
where A22 = D1 +
C1 =
(D1 - C2) -1 exists.
(21)
Artificial Damping Techniques Ill
IIit
llll II]1
IIII IIII| Illil[ •
IIIl!ll
~I
F12
|
f~2
Figure 3.2. Decomposition of the domain into two t~ subdomains.
The matrix A is called the enhanced matrix of A. It is not difficult d to check that (21) implies the existence o f . ~ - I (see [13]). Now, (15) can read
ifl Ul_X
- / u~/ f~ / u~-~ = M/u~x / e~ +-~ /-~-1 L u~ J where ~_
f3
[AI~
A12
0
0
[A21
Vl
0
0
0 D2
A23
0
A32
A33
[Oo
Lu~ -~ 0
_~=_
0 A21
'
0
(22)
000] 0
C1
C2 0
0 0
A23
Here we considered the following matrix splitting: A -- M - N.
(23)
In general, a domain decomposition method of the form (15) can be written, in its algebraic system, as (23). It is known that algorithm (22) converges if and only if 1
4. THE
ALGORITHM
In this section, we will consider a variant of algorithm (4). In it, we will employ an i n c o m p l e t e iteration of (15) and the coefficient of artificial damping ~ on the right-hand side of (1ha) will be set to be small for early outer iteration steps. T h a t is, algorithm (15) can be reformulated by replacing (1ha) with - A h u m+nl'j
+ i~2u~+1'~. = f +
-- K 2 u ?
i~?2mu ? ,
x • f~j,
where 1 < n < n . , for some n . > 1, ~/m -> 0, and ~rn --* 7" Now, we are ready to present the iterative algorithm. For a given initial guess uj0 ~_ ~ 0 [ ~ j , j = 1 , . . . , M , find u ~ +l'n, 1 < n < n., m > 0, such that (a) --
A h u jm + l ' n - K2u~n+l'n i~?2u~ +1''~ -. 2 ' -[" : f'~rnUj • U rn'bl'n
,
X E ~'~j,
i¢Ou m + l ' n -~ O,
X E rj,
(b)
v¢,3 ~
(c)
0f,jk u7 +1'~ + iZu7 +1'~ = -0b,kj ~ . ~ , ~ - 1 +,.u~-~ m.l,~-,,
(d)
+
m,n.
c J
(24)
u ~ +1'° = u ~ '~*,
x • f13,
where ~m >_ 0,
~m --* ~?,
and
uj0,n. = u j0,
x e r~k,
j = 1,..., M.
8
S. KIM AND M. LEE
It is not difficult to show that algorithm (24) is equivalent to the following iteration: for given ~0, find ~m+l, m > 0, by solving (a)
}m
(b) (c)
if II~mll < ~, stop; for ~m,0 = 0, find p-m,n , 1 < n < n., such that
(d)
M~"~m,n = ~ m , n - 1 + ym + i (din - d) ~m; ~m+l = ~m + ~)m,n.,
(25)
where
dm >_0 and
dm ~ d.
We are interested in determining the analytical convergence of algorithm (25). It is numerically verified that when the parameters d and dm are well chosen, algorithm (25) converges much faster (in computation time) than iterative algorithms which are not incorporating rational iterations. We will consider a heuristic method of finding efficient parameters in Section 5. One can apply algorithm (25) to real-valued symmetric positive definite systems of the form (3). (The term id should be replaced by d.) For such real-valued systems, one can prove the convergence of the corresponding algorithm. 5. N U M E R I C A L
RESULTS
In this section, we present numerical results to show the effectiveness of algorithm (24). The five-point finite difference method of uniform mesh size h is considered. The program is implemented in F O R T R A N with complex double precision, and is performed on an nCUBE2 parallel computer having 64 processors. We choose the following functions for the wave speed c(x, y): C 1 (X, y ) ---- 2,
c2(x, y) = eXY(2 + sin(~rx))(2 - sin(41ry)), c3(x, y)
f 3, I 1 + ex + sin(21rxy),
x < 0.5, elsewhere,
(26)
C4(X,y)={ 3,1.5, elsewhere.if x / ( x - 1)2 + ( Y - 1)2 < 2' The quantity q is chosen to be a nonnegative constant. The source function f is selected such that the true solution u(x, y) = ¢(x)"~2¢(y) , (2r) where ¢(x) = e i~(~-1) + e -i~x - 2. The right-hand side of (24b) is set correspondingly for a computational purpose. The error is estimated on the relative maximum norm r m, and the iteration is stopped when s~o _< 10 -4, where m
roo =
m
IIU IILo¢ (n)
,
=
ilUmll Lo¢ (f~)
,
where U m is the approximate solution of the m t h (outer) iteration. Zero initial values, U ° -- 0, are given. For a given artificial damping parameter 7, we choose rim----min (1, ~ ) . r ] , for some ~. (For computationally efficient ~'s, see (30) below.)
(28)
Artificial D a m p i n g Techniques
9
Each processor is assigned to solve a subproblem corresponding to one subdomain directly by using a Gauss elimination. For all results presented in this section, Mx [respectively, My] denotes the number of subdomains in x- [respectively, y-]direction, P is the number of processors used to perform the computation, and Tp (seconds) is the wall clock time to perform the given job on P processors. For the relaxation parameter ~ in (24), we applied the automatic procedure in [14], called ADOP. It is numerically tested that A D O P gives faster convergences t h a n any other (locally) constant complex-valued parameters and the convergence rate does not deteriorate as the mesh size h decreases. In Table 5.1, we compare the rate of convergence for various artificial damping ~?, when w = 50, q --- 0, 1/h == 128, and Mx x My -- 8 x 8. We choose ~ -- 1 and n. = 4. The number of outer iterations m, the m a x i m u m relative error r ~ , and the wall clock time T64 are presented for various z] and the wave speeds, ck, k -- 1, 2 and 3. Algorithm (24) did not converge when c = c2 and ~ is 0 or 4. Note that we do not know the convergence of (15) when q = ~ = 0 and t h a t we adopted an incomplete iteration. The superscript * indicates the total number of domain decomposition iterations when the artificial damping technique is not applied. From the table, one can see t h a t the artificial damping algorithm gives a fast convergence if the p a r a m e t e r is properly selected and the incomplete iteration is applied. We tested the classical iterative algorithms (relaxations and extrapolations) to solve the problems (TI = 0) presented in the table; they did not converge. Note that for such problems, conjugate gradient-type algorithms converge very slowly or m a y have possible breakdowns [1,2,15]. Table 5.1. w - - 5 0 , q - - 0 , C ~
7/
Ttl
1/h=
128, M~ x M u = 8 x S , ~ =
C1
t om o
C ~
T64
m
C2
1, a n d n , C ~
T om o
T64
m
=4.
C3
t om o
T64
0
313"
1.6e--2
90.4
206*
2.7e-2
60.5
4
39
1.6e-2
43.4
26
2.6e-2
30.0
6
23
1.6e--2
26.8
65
1.5e-2
70.5
24
2.6e--2
28.0
8
18
1.6e--2
21.6
38
1.5e-2
42.6
28
2.6e--2
32.1
i0
19
1.6e--2
22.6
53
1.5e-2
58.1
36
2.6e--2
40.4
12
22
1.6e--2
25.7
71
1.4e-2
76.7
46
2.7e-2
50.7
In numerical simulation of the wave problems of the form (1), one needs to choose the mesh size h such t h a t h(w/c) is not larger than 2/3 to 1 for a stability reason [7]. This choice of h corresponds to choosing at least 6 to 9 grid points per wavelength (21rc/w). In practice, 12 to 25 grid points per wavelength are often selected for accuracy reasons [7,16]. From practical computations, it is observed t h a t the rate of convergence of (24) is relatively sensitive to the p a r a m e t e r r/. (For the sensitivities to ~ and n,, see Table 5.3.) The quantity Q:=
Re ( K 2) I m ( K 2)
w2 c2
1 q2
is called the quality factor for the media, where K is the wave number defined in (12). For rocks in the earth, it is known t h a t seismic/acoustic waves propagate with the speed c being between 1 to 5 K m / s e c and Q between 100 to 1000. It is numerically verified t h a t one can select the p a r a m e t e r ~ for such problems as follows: Io:
~7
.
.
.
3c
.
lw .
.
6c
.
(29)
In Table 5.2, we check the convergence rate of algorithm (24) for various mesh sizes and wave speeds, when w = 75, q -- 3, and M~ × My = 32 × 2. The mesh sizes are selected as 1/96, 1/192, and 1/384. The parameters regarding the artificial damping technique are chosen as follows: CAHIa-A 31:8-B
10
S. KIM AND M. LEE
Table 5.2. w = 7 5 , q = 3 , Mx x M ~ = 3 2 x 2 , C=Cl
r/=10,~=1, andn. =5. c=e3
c=e2
1/h
m
m %0
T64
m
m roo
T64
m
m too
T64
96 192 384
30 28 28
6.3e-2 1.5e-2 3.9e-3
9.8 41.7 258.6
54 52 52
7.1e-2 1.7e-2 4.5e-3
16.9 76.5 472.3
37 36 36
7.7e-2 1.9e-2 5.0e-3
11.9 53.4 330.3
Table 5.3. w = 75, q ----3, c = c2, ~1= 10, 1/h = 192, and Mx x My = 16 x 4. n. = 1 N 308 308 295 274
n. : 2
~4
N
158.1 157.9 151.4 141.5
212 206 192 194
n. = 4 ~4
103.1 100.3 93.7 95.1
n, = 6
N
~4
N
~4
136 136 144 152
65.6 65.6 69.3 73.0
156 174 186 198
73.8 81.7 87.0 92.5
---- 10, ~ -- 1, a n d n . -- 5. I t seems t h a t t h e r a t e of convergence of t h e i t e r a t i v e a l g o r i t h m d o e s n o t d e t e r i o r a t e as t h e m e s h size decreases. F r o m t h e t a b l e we can see t h e s e c o n d - o r d e r convergence in L ° ° - n o r m for t h e a p p r o x i m a t e solutions. T a b l e 5.3 shows n u m e r i c a l results for various n . ' s a n d ~'s. T h e p r o b l e m coefficients are chosen as w = 75, q = 3, c = c2, a n d 77 = 10. W e set 1 / h = 192 a n d M x x M y = 16 x 4. T h e n u m b e r N d e n o t e s t h e t o t a l n u m b e r of d o m a i n d e c o m p o s i t i o n iterations. I n t h e table, one can see t h a t t h e a l g o r i t h m converges t h e b e s t w h e n n . -- 4 a n d ~ -- 1 or 2. I t is o b s e r v e d from v a r i o u s n u m e r i c a l s i m u l a t i o n s t h a t t h e a l g o r i t h m converges r a p i d l y w h e n ~] is chosen as in (29), 1
n.~-~.(M~-t-My),
and
~--1~2.
(30)
T h e a l g o r i t h m is t e s t e d t o solve t h e p r o b l e m s where p o i n t sources f are applied. W e choose f ( x ) -- ~ ( x 0 - x ) , where x = ( z , y ) a n d x0 -- (0.5,0.5), c ( x , y ) = c 4 ( x , y ) , w = 60, q -- 1, ~1 -- 8, = 1, a n d n . = 3. W h e n 1 / h = 80 a n d M s x M y = 8 x 8, t h e a l g o r i t h m (24) t a k e s 5.9 seconds (19 o u t e r i t e r a t i o n s ) on 64 processors. F i g u r e 5.1 p r e s e n t s t h e c o m p u t e d solution. I t seems t h a t m e s h r e f i n e m e n t s are r e q u i r e d n e a r t h e sources a n d t h e p o i n t s w h e r e t h e p r o b l e m coefficients are d i s c o n t i n u o u s . In view of saving c o m p u t e r s t o r a g e a n d avoiding c o m p l i c a t e d d a t a s t r u c t u r e s , it is i n t e r e s t i n g ( a n d r e a s o n a b l e ) to i n c o r p o r a t e local grid refinement over s u b d o m a i n s w h e r e t h e s t r o n g v a r i a t i o n s are e x p e c t e d a n d r e t a i n coarser grids elsewhere. D o m a i n d e c o m p o s i t i o n m e t h o d s (or t h e i r basic concepts) p r o v i d e a s y s t e m a t i c a n d e l e g a n t w a y t o i m p l e m e n t t h e a b o v e ideas. Now a v a r i a n t of t h e a l g o r i t h m (24) which i n c o r p o r a t e s finite e l e m e n t m e t h o d s a n d local grid r e f i n e m e n t s is b e i n g t e s t e d to solve some realistic s c a t t e r i n g p r o b l e m s .
6. C O N C L U S I O N S W e have c o n s i d e r e d an artificial d a m p i n g a l g o r i t h m for t h e H e l m h o l t z p r o b l e m . A convergence a r g u m e n t for t h e i t e r a t i v e a l g o r i t h m has b e e n presented. E a c h s t e p of t h e i t e r a t i o n was solved b y a d o m a i n d e c o m p o s i t i o n i t e r a t i v e p r o c e d u r e incompletely. So t h e a l g o r i t h m c a n b e viewed as a domain decomposition method accelerated by the rational iteration. T h e a l g o r i t h m was i m p l e m e n t e d on an n C U B E 2 a n d some n u m e r i c a l results were p r e s e n t e d . W e have discussed s t r a t e g i e s for finding c o m p u t a t i o n a l l y efficient a l g o r i t h m p a r a m e t e r s . I t has b e e n o b s e r v e d t h a t o u r a l g o r i t h m converges rapidly. (In fact, it showed faster convergences t h a n a n y o t h e r classical i t e r a t i v e m e t h o d s we tested: r e l a x a t i o n s , e x t r a p o l a t i o n s , C G - t y p e a l g o r i t h m s . ) F o r g e o p h y s i c a l a p p l i c a t i o n s , we should often solve t h e p r o b l e m (1) w i t h w b e i n g up t o a p p r o x i m a t e l y 600 a n d t h e wave s p e e d c h a v i n g big j u m p s . T h i s requires d i s c r e t e s o l u t i o n s in a
1l
Artificial Damping Techniques
Re(u)
0.00006 0.0000
o.oooo o /
-o.oooo2L/
/
1 Km"
ira(u)
O~ -0.0000
I Km
1 Km~ Figure 5.1. The solution for the point source f ( x ) = ~ (x0 - x), x0 = (0.5,0.5), w h e n c = c 4 , w = 6 0 , q = l , 7; = 8, ~ ---1, n* = 3 , 1 / h = 8 0 , a n d M x × M v = 8 × 8 .
piecewise linear approximation space of a mesh size h = 1/1500. (One can find many advantages for higher-order methods [16].) For such realistic problems, more efficient computational methods are to be found. REFERENCES 1. V. Faber and T. Manteuffel, Necessary and sufficient conditions for the existence of a conjugate gradient method, S I A M J. Numer. Anal. 21, 352-362 (1984). 2. W.D. Joubert and D.M. Young, Necessary and sufficient conditions for the simplification of the generalized conjugate-gradient algorithms, Linear Algebra AppL 8 8 / 8 9 , 449-485 (1987). 3. R.W. Freund, Conjugate gradient-type methods for linear systems with complex symmetric coefficient matrices, S I A M J. Sci. Stat. Comput. 13,425-448 (1992). 4. W. Hackbusch, Multigrid Methods and Applications, Springer-Verlag, Berlin, (1985). 5. J.H. Bramble, J.E. Pasciak and J. Xu, The analysis of multigrid algorithms for nonsymmetric and indefinite problems, Math. Comput. 51, 389-414 (1988). 6. C.C. Douglas and J. Douglas, Jr., A unified convergence theory for abstract multigrid or multilevel algorithms, serial and parallel, S I A M J. Numer. Anal. 30, 136-158 (1993).
12
S. KIM AND M. LEE
7. V.V. Shaidurov and E.I. Ogorodnikov, Some numerical methods of solving Helmholtz wave equation, In Mathematical and Numerical Aspects of Wave Propagation Phenomena, (Edited by G. Cohen, L. Halpern and P. Joly), pp. 73-79, SIAM, Philadelphia, (1991). 8. S. Kim, On the use of rational iterations and domain decomposition methods for solving the Helmholtz problem (submitted). 9. S. Kim, Domain decomposition iterative procedures for solving scalar waves in the frequency domain (submitted). 10. P.L. Lions, On the Schwarz alternating method III: A variant for nonoverlapping subdomains, In Third International Symposium on Domain Decomposition Method/or Partial Differential Equations, (Edited by T.F. Chan, R. Glowinski, J. Periaux and O.B. Widlund), pp. 202-223, SIAM, Philadelphia, (1990). 11. B. Despr~s, Domain decomposition method and the Helmholtz problem, In Mathematical and Numerical Aspects of Wave Propagation Phenomena, (Edited by G. Cohen, L. Halpern and P. Joly), pp. 44-51, SIAM, Philadelphia, (1991). 12. J. Douglas, Jr., P.S. Paes Leme, J.E. Roberts and J. Wang, A parallel iterative procedure applicable to the approximate solution of second order partial differential equations by mixed finite element methods, Numer. Math. 65, 95-108 (1993). 13. W.P. Tang, Generalized Schwarz splittings, SIAM J. Sci. Stat. Comput. 13, 573-595 (1992). 14. S. Kim, A parallelizable iterative procedure for the Helmholtz problem, Appl. Numer. Math. 14, 435-449 (1994). 15. A. Bayliss, C. Goldstein and E. Turkel, An iterative method for the Helmholtz equation, J. Comput. Phys. 49, 443-457 (1983). 16. A. Bayliss, C. Goldstein and E. Turkel, On accuracy conditions for the numerical computation of waves, J. Comput. Phys. 59, 396-404 (1985).
Pergamon
Copyright(g)1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved S0898-1221(96)00025-9 0898-1221/96 $15.00 + 0.00
Artificial D a m p i n g Techniques for Scalar Waves in the Frequency D o m a i n S. KIM Department of Computational and Applied Mathematics, Rice University P.O. Box 1892, Houston, TX 77251-1892, U.S.A. skim~caam, rice. edu
M. LEE Center for Applied Mathematics, Purdue University W. Lafayette, IN 47907, U.S.A. miyoung~math, purdue, edu (Received October 1995; accepted November 1995)
A b s t r a c t - - A n artificial damping algorithm for solving the Helmholtz problem is considered. When the imaginary part of the wave number is small, the problem is known to be difficult to solve. In this paper, we propose an efficient artificial damping algorithm which can be viewed as a rational iteration. Each damped problem is solved incompletely by a nonoverlapping domain decomposition method. Convergence arguments are presented• Numerical results run on an nCUBE2 are presented to demonstrate the effectiveness of the method• K e y w o r d s - - A r t i f i c i a l damping, Wave propagation, Robin boundary condition, Domain decomposition method•
1. I N T R O D U C T I O N Let Ft = (0, 1) d, d = 2 or 3, with the b o u n d a r y F = 0Ft. Consider the following (complex-valued) scalar wave problem: 022 xEft, - A u - c~x)2 U -t- i q ( x ) 2 u = f ( x ) , cgu
+
w
c(x)
(1) u--0,
xEF,
where i is the imaginary unit, w > 0 is the angular frequency, c denotes the wave speed in the m e d i u m where the wave is propagating, q2 __ bw, where b > 0 is the friction coefficient, and u is the unit o u t w a r d normal to F. Here the coefficients are sufficiently regular so t h a t the existence and uniqueness of a solution of (1) lying in HS(ft) for some s > 1 for reasonable f are assured. T h e second equation of (1) represents a first-order outgoing radiation condition t h a t allows n o r m a l l y incident waves to pass out of Ft transparently. T h e problem (1) models the p r o p a g a t i o n of t i m e - h a r m o n i c waves such as electromagnetic waves, discretizations of the t i m e - d e p e n d e n t Schr5dinger equation by implicit difference schemes, inverse scattering problems, seismic waves, and underwater acoustics. In this paper, we are mainly interested in the case q -- 0,
or q is small.
(2)
I n this case, the linear system arising in the finite difference discretization of the problem (1) is difficult t o solve. In addition to having a complex-valued solution, the problem is neither coercive Typeset by .AAdS-TEX
2
S. KIM AND M. LEE
nor Hermitian symmetric. As a consequence, most standard iterative methods either converge slowly or fail to converge. (For q = 0, there is no convergent conjugate gradient (CG)-type algorithm for the problem [1-3].) For an iterative algorithm for solving (1), a coarse grid solution can be used to obtain a better initial guess a n d / o r to accelerate the convergence of the iteration. Such an idea has been widely studied in terms of multigrid (MG) methods where one uses a collection of successive coarser finite element grids; see [4-6] and references therein. There are two steps in MG methods, coarse grid correction and smoothing. The MG method is known to be effective for solving certain coercive elliptic problems. The reason why the MG method is so efficient is, roughly speaking, the following: in the coarse grid correction, the low and medium frequency components of the error (corresponding to small and medium eigenvalues) are significantly reduced, while the smoothing step reduces the high frequency components of the error. Since one needs to choose the mesh size h such that h. maxxef~(w/c(x)) is less than 2/3 to 1 for stability reasons [7], which is the same as taking at least 6 to 9 points per wavelength (27rc/w), it is hard to solve problem (1) even on the coarsest grid for large wave numbers. (So MG methods cannot be applied to (1) efficiently without designing a solver for the coarsest grid problem.) For such a difficult problem, we will propose a rational iteration instead of considering MG methods or CG-type accelerations. An outline of this paper is as follows. In Section 2, we propose an iterative algorithm that introduces artificial damping/attenuation and that can be viewed as a rational iterative algorithm. A convergence argument for the algorithm is presented. In Section 3, we consider a nonoverlapping domain decomposition iterative procedure to solve a fractional step of the rational algorithm presented in Section 2. So the domain decomposition iteration is considered as the inner iteration, while the rational iteration is the outer iteration. In Section 4, we present the main algorithm in which the inner iteration finishes incompletely and the coefficient of artificial damping is modified in early steps of the outer iteration. The method is implemented in an nCUBE2 parallel computer and some results are given in Section 5 to demonstrate the effectiveness of the method. In this section, we present heuristic strategies for finding efficient algorithm parameters introduced in Section 4. The last section includes conclusions and possible applications.
2. A R T I F I C I A L
DAMPING
When problem (1) is discretized by a finite difference/element method of a mesh size h, it can be represented as follows: find u satisfying Au = k.
(3)
Here A (= A h) is a symmetric complex matrix and k is the force vector. It is well known that the matrix A in (3) is poorly conditioned; as a consequence, most standard iterative methods either fail to converge or converge so slowly as to be impractical. Our next objective is to present an iterative algorithm for (3) which introduces artificial damping. First, we cite the following lemma [8]. LEMMA 2.1. Let A be an eigenvalue o[ A. Then, Im(~) > 0 for w2h sut~ciently small. Now, consider the following iterative procedure: for a given initial guess u °, solve the following, m_>O: (a) r m = k - Aura; (b)
if Ilrmll < e, stop;
(c) (d)
find pm such that (A + idI) p m = rra; um+l
(4)
= U m _1_ p r o .
Here 6 is the stopping criterion for the iteration, d is the coefficient of artificial damping, I is the identity matrix, and I1" II denotes the e2-norm. In Section 5, we will consider a heuristic strategy for finding an efficient coefficient of damping.
Artificial D a m p i n g T e c h n i q u e s
3
THEOREM 2.2. The iteration (4) converges for all d > O. PROOF. F r o m (4a), (4c), and (4d), U m + l 7__ U m + p m
= u m + (A + idI) -1 (k - A u m) =
(5)
+ (A + i d I ) - l A (u - urn).
U m
Let e m = u - um. Then, from (5), one has e m+l = e m - (A + i d X ) - l A e m
(6)
-- i d ( A + idI) -1 em. Let a(A) be the s p e c t r u m of the m a t r i x A. Then, Ilem+lH
<
max
- -
. Ilem]].
AEcr(A) A i/id
It follows from L e m m a 2.1 t h a t (4) converges for all d > 0.
|
REMARK. In algorithm (4), the p r o b l e m is how one can invert A + i d I efficiently. Since A + i d I is b e t t e r conditioned, one can find easily a convergent iterative algorithm for (4c). We do not have to solve (4c) directly/completely, in practical simulations. T h e p r o b l e m can be solved incompletely using an iterative algorithm. REMARK. T h e algorithm (4) can be rewritten as follows. Multiply (4d) by A + i d I and a p p l y (4a) and (4c) to the resulting equation. Then,
(A + idI) u "~+1 = k + i d u m.
(7)
Again, equation (7) can be rewritten as a discrete version of --Aum+l -
( w2 -~-iq
2) u m + l + i ~ 2 u m + l = f ( x ) + i ~ ? 2 u ra,
x•~,
(8)
Ou m+l
O~
+ iw-um+lc
= O,
x • F,
for some U > 0. For the 5-point finite difference discretization of a uniform grid size h, we can see 772 = d/h 2. In Section 3, we will consider a domain decomposition m e t h o d for the d a m p e d p r o b l e m (8). REMARK. Let Ud be the solution of
(A + idI) Ud = k.
(9)
T h e n , from (3) and (9), u
--
U d
----
id (A + i d I ) - 1 u,
and, therefore, Iku - udll <
max
.
(10)
T h e p e r t u r b e d solution u a can be a good initial guess for u; see Figure 2.1. T h e r e finite difference solutions of (9) on the line segment Y0.5 = {(Xl, x2) • gt : x2 = 0.5} are depicted for the point source f = 5(x0 - x), x0 = (0.5,0.5), and for the uniform grid size h = 1/100. We selected the p r o b l e m coefficients such t h a t q - 0, w = 70, and c = c4 in (26) below. F r o m Figure 2.1, one can see t h a t even t h o u g h the m a g n i t u d e of Ud is different from t h a t of u, the phase of Ud seems to be identical to t h a t of u.
S. K I M
AND
M.
LEE
Re(u}:
y=0.5
0.00004 - -
: d=O
..... : d = l
0.00003
0.00002
0.00001
-0.00001 i
.
.
.
.
0
i
.
.
.
.
0.2
i
i
,
,
0.4
i
!
J
.
•
•
'
0.6
.
.
.
.
'
0.8
,=0.5
Im{u)
0.00001
~ J
-0.00001
----: -0.00002
d=0
..... : d = l
0
0.2
0.4
0.6
0.8
X
Figure 2.1. Finite difference solutions of (9) on Y0.5 for the point source f = x o ---- ( 0 . 5 , 0 . 5 ) , w h e n q --- 0, w = 70, c = c4 in ( 2 6 ) below, and h = 1 / 1 0 0 .
3. S O L U T I O N
FOR
THE
DAMPED
6(xo -x),
PROBLEM
In this section, we will consider an iterative algorithm for inverting A + idI in (4c), or equivalently, for solving (8). We will employ the domain decomposition method in [9] for solving (8). Let {flj, j = 1 , . . . , M } be a partition of ~: M
j=l
Assume that f2j, j = 1 , . . . , M, are also rectangular/cubic regions. Let
r~ = r ~ O~j,
r~k = rk3 = o~j ~ O~k.
Let us consider the decomposition of the problem (8) over {~2j }. It can be easily checked that
Artificial Damping Techniques
5
the problem (8) is equivalent to the following. Find u~ +1, j = 1 , . . . , M, such that (a)
- - A U ? + ' --
K 2 u ? +1 + i772u? +1 = f + i~2u?,
x E fly,
OU? +1
(b)
O--~j + i~um+lc3
= O,
0 U ? +1
(c)
o.j
X E Fj,
(11)
C~U~n+l
+ iZ~'[+~ -
o.k
+ iZu~+l'
z e rjk,
where K is the wave number defined as K2-
022
c2
iq 2,
(12)
vj is the unit outward normal from flj, and fl is a complex function on the subdomain interfaces with Re(fl) > 0. In (11), the consistency conditions u~+l = u~+l'
Ovj
+
cgvk
O,
on Fjk
(13)
are replaced by the Robin interface boundary condition (11c). This replacement is often more convenient and it is not difficult to show that (11c) and (13) are equivalent; see [8-12]. The domain decomposition iterative algorithm for (11) can be defined as follows. Choose an initial guess u~n+l'°, j = 1 , . . . , M, and then build u~n+l'n, n > 1, recursively by solving (a) - - A u ? + l ' n - K 2 u ? +l'n + i~2u? +l'n = f + i~2u?, 0_ u ?_+ l ' n
(b)
~- $ 03 c U jra+l n --~ O,
x C F j,
Ovj 01U, n?_+_ 0/]j
(C)
-~- ifl U ? + l ' n __
x E ftj,
GUkO rn+l ,n-1 .~ rn+l,n-- 1 0/]k + zp u k
(14)
x E Fjk,
where u ~ is the limit of sequence {ujm n' : n _ 0}. The above algorithm was analyzed in [11] in differential, rather than discrete, level. Note that the Robin interface boundary condition (14c) imposes the continuity of the solution u and its flux. Since most finite element/difference methods admit discontinuities of the flux on the element interfaces, (14c) should be modified in discretized problems. We will approximate (14) by the five-point finite difference method of a uniform grid size h. Let A h u j be the centered five-point difference approximation of A u j , O¢,juj be the centered difference for ~ on the boundary Fj, and ~f, j k U j and Ob,jkuj denote the forward and the backward finite differences for ~ on the interface Fjk, respectively. (Here we assume an exterior bordering of the subdomains.) We can define a finite difference domain decomposition iterative algorithm as follows: for given u ~ +1'°, j = 1 , . . . , M, solve ujm+l,n , n ~ 1, satisfying ^ .~+1,,~ _ K 2 u ~ + l , n + i , 2 u ~.+ l , n = f + z~ 2uj.~ , (a) --,,hUj
(b) (C)
i W u m + l ''~ = O, O¢,j u ~ +1''~ + c 3 0
f, j k
um+l,n i~u?+l,n m + l , n - 1 -- ,~ m + l n - 1 j + = --~o,kj Uk ~ $~OUk ,
x C ~j, x E F j,
(15)
x E Fjk.
Algorithm (15) can be viewed as a generalized Schwarz splitting with minimum overlap [13]. Note that the modified Robin interface boundary condition (15c) imposes the continuity of the discrete solution only. One can prove the convergence of algorithm (15) for q2 + r]2 ;> O.
6
S. KIM AND M. LEE
A convergence argument and strategies for finding efficient parameter/9 for (15) can be found in [9,14]. It is numerically verified that the convergence rate of (15) is independent of the grid size h. Some applications of the algorithm to finite element methods can be found in [9]. In the following, we will illustrate the algebraic system of algorithm (15). First, let us consider the algebraic equations of (15) corresponding to interface nodal points; see Figure 3.1. Let F = f + i y 2 u ~ and drop superscript m and m + 1, for convenience. At the point C, (15a) reads (4 - h 2 K ~ + ih2zl~) uj~,c - u jn, E
- - U jn, W
- - U j,~ , S - - U jn, N
h2Fc,
(16)
where ujn,c = u~(C), the value of u~ at the point C, and the others are similarly defined. The term ujn,E in (16) evaluated at a point out of the subdomain f~j can be substituted by using (15c). Equation (15c) is written as n
n
n-1
uJ'r' - Uj'c + i/9 n ~tj, c h
n--1
Uk,E -- Uk,c + i / 3 h
--
,~-1
Uk,C ,
or, equivalently,
(17)
uj,Bn - (1 - i/~h) u~, c = Uk,En-1 _ (1 -- i/3h) u n-ik,c. Substituting (17) into (16), we have [4 - h 2 g 2 + ih2~l~ - (1 - i/3h)] uj~,c - ujn,w - uj~,s - Uj~,N n-1 = h 2 F c + Uk,E
_
(1 -
it3h) u ~ 1.
(18)
In the same manner, one can treat cross points arising in a box-type decomposition of the domain. Clearly, algorithm (15) is applicable to nonuniform meshes with minor changes. N w
c
E
S f/j
Fjk
~;,
Figure 3.1. Five point stencil at interface nodal points.
Now, we consider the matrix representation of algorithm (15). For simplicity, let the domain be decomposed into two subdomains; see Figure 3.2. Give a nodal point ordering as follows: we count nodes in f~l the first, those on F12 next, and those in ~2 the last. Then, by dropping the superscripts m and m + 1, one can rewrite (8) as
LorA11 A120 [Ul] If1] A22 A23
u2=
A32
u3
A33J
(19)
f3
°IUl [fl
Problem (11) can be represented by the following enhanced system: All A21 A21
A12 D1 C2
0 C1 D2
A23
u2
A23
u2
0
0
A32
A33
u3
D2 + C2,
and
=
f2 f2 f3
'
(20)
where A22 = D1 +
C1 =
(D1 - C2) -1 exists.
(21)
Artificial Damping Techniques Ill
IIit
llll II]1
IIII IIII| Illil[ •
IIIl!ll
~I
F12
|
f~2
Figure 3.2. Decomposition of the domain into two t~ subdomains.
The matrix A is called the enhanced matrix of A. It is not difficult d to check that (21) implies the existence o f . ~ - I (see [13]). Now, (15) can read
ifl Ul_X
- / u~/ f~ / u~-~ = M/u~x / e~ +-~ /-~-1 L u~ J where ~_
f3
[AI~
A12
0
0
[A21
Vl
0
0
0 D2
A23
0
A32
A33
[Oo
Lu~ -~ 0
_~=_
0 A21
'
0
(22)
000] 0
C1
C2 0
0 0
A23
Here we considered the following matrix splitting: A -- M - N.
(23)
In general, a domain decomposition method of the form (15) can be written, in its algebraic system, as (23). It is known that algorithm (22) converges if and only if 1
4. THE
ALGORITHM
In this section, we will consider a variant of algorithm (4). In it, we will employ an i n c o m p l e t e iteration of (15) and the coefficient of artificial damping ~ on the right-hand side of (1ha) will be set to be small for early outer iteration steps. T h a t is, algorithm (15) can be reformulated by replacing (1ha) with - A h u m+nl'j
+ i~2u~+1'~. = f +
-- K 2 u ?
i~?2mu ? ,
x • f~j,
where 1 < n < n . , for some n . > 1, ~/m -> 0, and ~rn --* 7" Now, we are ready to present the iterative algorithm. For a given initial guess uj0 ~_ ~ 0 [ ~ j , j = 1 , . . . , M , find u ~ +l'n, 1 < n < n., m > 0, such that (a) --
A h u jm + l ' n - K2u~n+l'n i~?2u~ +1''~ -. 2 ' -[" : f'~rnUj • U rn'bl'n
,
X E ~'~j,
i¢Ou m + l ' n -~ O,
X E rj,
(b)
v¢,3 ~
(c)
0f,jk u7 +1'~ + iZu7 +1'~ = -0b,kj ~ . ~ , ~ - 1 +,.u~-~ m.l,~-,,
(d)
+
m,n.
c J
(24)
u ~ +1'° = u ~ '~*,
x • f13,
where ~m >_ 0,
~m --* ~?,
and
uj0,n. = u j0,
x e r~k,
j = 1,..., M.
8
S. KIM AND M. LEE
It is not difficult to show that algorithm (24) is equivalent to the following iteration: for given ~0, find ~m+l, m > 0, by solving (a)
}m
(b) (c)
if II~mll < ~, stop; for ~m,0 = 0, find p-m,n , 1 < n < n., such that
(d)
M~"~m,n = ~ m , n - 1 + ym + i (din - d) ~m; ~m+l = ~m + ~)m,n.,
(25)
where
dm >_0 and
dm ~ d.
We are interested in determining the analytical convergence of algorithm (25). It is numerically verified that when the parameters d and dm are well chosen, algorithm (25) converges much faster (in computation time) than iterative algorithms which are not incorporating rational iterations. We will consider a heuristic method of finding efficient parameters in Section 5. One can apply algorithm (25) to real-valued symmetric positive definite systems of the form (3). (The term id should be replaced by d.) For such real-valued systems, one can prove the convergence of the corresponding algorithm. 5. N U M E R I C A L
RESULTS
In this section, we present numerical results to show the effectiveness of algorithm (24). The five-point finite difference method of uniform mesh size h is considered. The program is implemented in F O R T R A N with complex double precision, and is performed on an nCUBE2 parallel computer having 64 processors. We choose the following functions for the wave speed c(x, y): C 1 (X, y ) ---- 2,
c2(x, y) = eXY(2 + sin(~rx))(2 - sin(41ry)), c3(x, y)
f 3, I 1 + ex + sin(21rxy),
x < 0.5, elsewhere,
(26)
C4(X,y)={ 3,1.5, elsewhere.if x / ( x - 1)2 + ( Y - 1)2 < 2' The quantity q is chosen to be a nonnegative constant. The source function f is selected such that the true solution u(x, y) = ¢(x)"~2¢(y) , (2r) where ¢(x) = e i~(~-1) + e -i~x - 2. The right-hand side of (24b) is set correspondingly for a computational purpose. The error is estimated on the relative maximum norm r m, and the iteration is stopped when s~o _< 10 -4, where m
roo =
m
IIU IILo¢ (n)
,
=
ilUmll Lo¢ (f~)
,
where U m is the approximate solution of the m t h (outer) iteration. Zero initial values, U ° -- 0, are given. For a given artificial damping parameter 7, we choose rim----min (1, ~ ) . r ] , for some ~. (For computationally efficient ~'s, see (30) below.)
(28)
Artificial D a m p i n g Techniques
9
Each processor is assigned to solve a subproblem corresponding to one subdomain directly by using a Gauss elimination. For all results presented in this section, Mx [respectively, My] denotes the number of subdomains in x- [respectively, y-]direction, P is the number of processors used to perform the computation, and Tp (seconds) is the wall clock time to perform the given job on P processors. For the relaxation parameter ~ in (24), we applied the automatic procedure in [14], called ADOP. It is numerically tested that A D O P gives faster convergences t h a n any other (locally) constant complex-valued parameters and the convergence rate does not deteriorate as the mesh size h decreases. In Table 5.1, we compare the rate of convergence for various artificial damping ~?, when w = 50, q --- 0, 1/h == 128, and Mx x My -- 8 x 8. We choose ~ -- 1 and n. = 4. The number of outer iterations m, the m a x i m u m relative error r ~ , and the wall clock time T64 are presented for various z] and the wave speeds, ck, k -- 1, 2 and 3. Algorithm (24) did not converge when c = c2 and ~ is 0 or 4. Note that we do not know the convergence of (15) when q = ~ = 0 and t h a t we adopted an incomplete iteration. The superscript * indicates the total number of domain decomposition iterations when the artificial damping technique is not applied. From the table, one can see t h a t the artificial damping algorithm gives a fast convergence if the p a r a m e t e r is properly selected and the incomplete iteration is applied. We tested the classical iterative algorithms (relaxations and extrapolations) to solve the problems (TI = 0) presented in the table; they did not converge. Note that for such problems, conjugate gradient-type algorithms converge very slowly or m a y have possible breakdowns [1,2,15]. Table 5.1. w - - 5 0 , q - - 0 , C ~
7/
Ttl
1/h=
128, M~ x M u = 8 x S , ~ =
C1
t om o
C ~
T64
m
C2
1, a n d n , C ~
T om o
T64
m
=4.
C3
t om o
T64
0
313"
1.6e--2
90.4
206*
2.7e-2
60.5
4
39
1.6e-2
43.4
26
2.6e-2
30.0
6
23
1.6e--2
26.8
65
1.5e-2
70.5
24
2.6e--2
28.0
8
18
1.6e--2
21.6
38
1.5e-2
42.6
28
2.6e--2
32.1
i0
19
1.6e--2
22.6
53
1.5e-2
58.1
36
2.6e--2
40.4
12
22
1.6e--2
25.7
71
1.4e-2
76.7
46
2.7e-2
50.7
In numerical simulation of the wave problems of the form (1), one needs to choose the mesh size h such t h a t h(w/c) is not larger than 2/3 to 1 for a stability reason [7]. This choice of h corresponds to choosing at least 6 to 9 grid points per wavelength (21rc/w). In practice, 12 to 25 grid points per wavelength are often selected for accuracy reasons [7,16]. From practical computations, it is observed t h a t the rate of convergence of (24) is relatively sensitive to the p a r a m e t e r r/. (For the sensitivities to ~ and n,, see Table 5.3.) The quantity Q:=
Re ( K 2) I m ( K 2)
w2 c2
1 q2
is called the quality factor for the media, where K is the wave number defined in (12). For rocks in the earth, it is known t h a t seismic/acoustic waves propagate with the speed c being between 1 to 5 K m / s e c and Q between 100 to 1000. It is numerically verified t h a t one can select the p a r a m e t e r ~ for such problems as follows: Io:
~7
.
.
.
3c
.
lw .
.
6c
.
(29)
In Table 5.2, we check the convergence rate of algorithm (24) for various mesh sizes and wave speeds, when w = 75, q -- 3, and M~ × My = 32 × 2. The mesh sizes are selected as 1/96, 1/192, and 1/384. The parameters regarding the artificial damping technique are chosen as follows: CAHIa-A 31:8-B
10
S. KIM AND M. LEE
Table 5.2. w = 7 5 , q = 3 , Mx x M ~ = 3 2 x 2 , C=Cl
r/=10,~=1, andn. =5. c=e3
c=e2
1/h
m
m %0
T64
m
m roo
T64
m
m too
T64
96 192 384
30 28 28
6.3e-2 1.5e-2 3.9e-3
9.8 41.7 258.6
54 52 52
7.1e-2 1.7e-2 4.5e-3
16.9 76.5 472.3
37 36 36
7.7e-2 1.9e-2 5.0e-3
11.9 53.4 330.3
Table 5.3. w = 75, q ----3, c = c2, ~1= 10, 1/h = 192, and Mx x My = 16 x 4. n. = 1 N 308 308 295 274
n. : 2
~4
N
158.1 157.9 151.4 141.5
212 206 192 194
n. = 4 ~4
103.1 100.3 93.7 95.1
n, = 6
N
~4
N
~4
136 136 144 152
65.6 65.6 69.3 73.0
156 174 186 198
73.8 81.7 87.0 92.5
---- 10, ~ -- 1, a n d n . -- 5. I t seems t h a t t h e r a t e of convergence of t h e i t e r a t i v e a l g o r i t h m d o e s n o t d e t e r i o r a t e as t h e m e s h size decreases. F r o m t h e t a b l e we can see t h e s e c o n d - o r d e r convergence in L ° ° - n o r m for t h e a p p r o x i m a t e solutions. T a b l e 5.3 shows n u m e r i c a l results for various n . ' s a n d ~'s. T h e p r o b l e m coefficients are chosen as w = 75, q = 3, c = c2, a n d 77 = 10. W e set 1 / h = 192 a n d M x x M y = 16 x 4. T h e n u m b e r N d e n o t e s t h e t o t a l n u m b e r of d o m a i n d e c o m p o s i t i o n iterations. I n t h e table, one can see t h a t t h e a l g o r i t h m converges t h e b e s t w h e n n . -- 4 a n d ~ -- 1 or 2. I t is o b s e r v e d from v a r i o u s n u m e r i c a l s i m u l a t i o n s t h a t t h e a l g o r i t h m converges r a p i d l y w h e n ~] is chosen as in (29), 1
n.~-~.(M~-t-My),
and
~--1~2.
(30)
T h e a l g o r i t h m is t e s t e d t o solve t h e p r o b l e m s where p o i n t sources f are applied. W e choose f ( x ) -- ~ ( x 0 - x ) , where x = ( z , y ) a n d x0 -- (0.5,0.5), c ( x , y ) = c 4 ( x , y ) , w = 60, q -- 1, ~1 -- 8, = 1, a n d n . = 3. W h e n 1 / h = 80 a n d M s x M y = 8 x 8, t h e a l g o r i t h m (24) t a k e s 5.9 seconds (19 o u t e r i t e r a t i o n s ) on 64 processors. F i g u r e 5.1 p r e s e n t s t h e c o m p u t e d solution. I t seems t h a t m e s h r e f i n e m e n t s are r e q u i r e d n e a r t h e sources a n d t h e p o i n t s w h e r e t h e p r o b l e m coefficients are d i s c o n t i n u o u s . In view of saving c o m p u t e r s t o r a g e a n d avoiding c o m p l i c a t e d d a t a s t r u c t u r e s , it is i n t e r e s t i n g ( a n d r e a s o n a b l e ) to i n c o r p o r a t e local grid refinement over s u b d o m a i n s w h e r e t h e s t r o n g v a r i a t i o n s are e x p e c t e d a n d r e t a i n coarser grids elsewhere. D o m a i n d e c o m p o s i t i o n m e t h o d s (or t h e i r basic concepts) p r o v i d e a s y s t e m a t i c a n d e l e g a n t w a y t o i m p l e m e n t t h e a b o v e ideas. Now a v a r i a n t of t h e a l g o r i t h m (24) which i n c o r p o r a t e s finite e l e m e n t m e t h o d s a n d local grid r e f i n e m e n t s is b e i n g t e s t e d to solve some realistic s c a t t e r i n g p r o b l e m s .
6. C O N C L U S I O N S W e have c o n s i d e r e d an artificial d a m p i n g a l g o r i t h m for t h e H e l m h o l t z p r o b l e m . A convergence a r g u m e n t for t h e i t e r a t i v e a l g o r i t h m has b e e n presented. E a c h s t e p of t h e i t e r a t i o n was solved b y a d o m a i n d e c o m p o s i t i o n i t e r a t i v e p r o c e d u r e incompletely. So t h e a l g o r i t h m c a n b e viewed as a domain decomposition method accelerated by the rational iteration. T h e a l g o r i t h m was i m p l e m e n t e d on an n C U B E 2 a n d some n u m e r i c a l results were p r e s e n t e d . W e have discussed s t r a t e g i e s for finding c o m p u t a t i o n a l l y efficient a l g o r i t h m p a r a m e t e r s . I t has b e e n o b s e r v e d t h a t o u r a l g o r i t h m converges rapidly. (In fact, it showed faster convergences t h a n a n y o t h e r classical i t e r a t i v e m e t h o d s we tested: r e l a x a t i o n s , e x t r a p o l a t i o n s , C G - t y p e a l g o r i t h m s . ) F o r g e o p h y s i c a l a p p l i c a t i o n s , we should often solve t h e p r o b l e m (1) w i t h w b e i n g up t o a p p r o x i m a t e l y 600 a n d t h e wave s p e e d c h a v i n g big j u m p s . T h i s requires d i s c r e t e s o l u t i o n s in a
1l
Artificial Damping Techniques
Re(u)
0.00006 0.0000
o.oooo o /
-o.oooo2L/
/
1 Km"
ira(u)
O~ -0.0000
I Km
1 Km~ Figure 5.1. The solution for the point source f ( x ) = ~ (x0 - x), x0 = (0.5,0.5), w h e n c = c 4 , w = 6 0 , q = l , 7; = 8, ~ ---1, n* = 3 , 1 / h = 8 0 , a n d M x × M v = 8 × 8 .
piecewise linear approximation space of a mesh size h = 1/1500. (One can find many advantages for higher-order methods [16].) For such realistic problems, more efficient computational methods are to be found. REFERENCES 1. V. Faber and T. Manteuffel, Necessary and sufficient conditions for the existence of a conjugate gradient method, S I A M J. Numer. Anal. 21, 352-362 (1984). 2. W.D. Joubert and D.M. Young, Necessary and sufficient conditions for the simplification of the generalized conjugate-gradient algorithms, Linear Algebra AppL 8 8 / 8 9 , 449-485 (1987). 3. R.W. Freund, Conjugate gradient-type methods for linear systems with complex symmetric coefficient matrices, S I A M J. Sci. Stat. Comput. 13,425-448 (1992). 4. W. Hackbusch, Multigrid Methods and Applications, Springer-Verlag, Berlin, (1985). 5. J.H. Bramble, J.E. Pasciak and J. Xu, The analysis of multigrid algorithms for nonsymmetric and indefinite problems, Math. Comput. 51, 389-414 (1988). 6. C.C. Douglas and J. Douglas, Jr., A unified convergence theory for abstract multigrid or multilevel algorithms, serial and parallel, S I A M J. Numer. Anal. 30, 136-158 (1993).
12
S. KIM AND M. LEE
7. V.V. Shaidurov and E.I. Ogorodnikov, Some numerical methods of solving Helmholtz wave equation, In Mathematical and Numerical Aspects of Wave Propagation Phenomena, (Edited by G. Cohen, L. Halpern and P. Joly), pp. 73-79, SIAM, Philadelphia, (1991). 8. S. Kim, On the use of rational iterations and domain decomposition methods for solving the Helmholtz problem (submitted). 9. S. Kim, Domain decomposition iterative procedures for solving scalar waves in the frequency domain (submitted). 10. P.L. Lions, On the Schwarz alternating method III: A variant for nonoverlapping subdomains, In Third International Symposium on Domain Decomposition Method/or Partial Differential Equations, (Edited by T.F. Chan, R. Glowinski, J. Periaux and O.B. Widlund), pp. 202-223, SIAM, Philadelphia, (1990). 11. B. Despr~s, Domain decomposition method and the Helmholtz problem, In Mathematical and Numerical Aspects of Wave Propagation Phenomena, (Edited by G. Cohen, L. Halpern and P. Joly), pp. 44-51, SIAM, Philadelphia, (1991). 12. J. Douglas, Jr., P.S. Paes Leme, J.E. Roberts and J. Wang, A parallel iterative procedure applicable to the approximate solution of second order partial differential equations by mixed finite element methods, Numer. Math. 65, 95-108 (1993). 13. W.P. Tang, Generalized Schwarz splittings, SIAM J. Sci. Stat. Comput. 13, 573-595 (1992). 14. S. Kim, A parallelizable iterative procedure for the Helmholtz problem, Appl. Numer. Math. 14, 435-449 (1994). 15. A. Bayliss, C. Goldstein and E. Turkel, An iterative method for the Helmholtz equation, J. Comput. Phys. 49, 443-457 (1983). 16. A. Bayliss, C. Goldstein and E. Turkel, On accuracy conditions for the numerical computation of waves, J. Comput. Phys. 59, 396-404 (1985).