An improved domain decomposition method for the 3D Helmholtz equation

An improved domain decomposition method for the 3D Helmholtz equation

Computer methods in applied mechanics and englneerlng Comput. ELSEVIER Methods Appl. Mech. Engrg. 162 (1998) 113- 124 An improved domain decompo...

728KB Sizes 19 Downloads 185 Views

Computer methods in applied mechanics and englneerlng Comput.

ELSEVIER

Methods

Appl. Mech. Engrg.

162 (1998)

113- 124

An improved domain decomposition method for the 3D Helniholtz equation A. Piacentini”, Department

N. Rosa

of Mathematics, University of Milan, via Saldini 50, 20133 Milano, Italy Received

5 May 1997; revised 28 October

1997

Abstract The efficiency of a non-overlapping domain decomposition method is strictly related to the transmission operator acting at the interfaces between subdomains. A new operator, arising from the non-reflecting boundary conditions theory, and a new iteration scheme lead to a dramatic improvement of the performances of the scheme which was originally proposed by Bnmo Desprks and which is recalled in the first

part of this paper. The convergence of the new method is proved and some numerical results show the comparison of the performances. 0 1998 Elsevier Science S.A. All rights reserved.

1. Introduction A great number of numerical simulations in acoustic analysis or in electromagnetism requires the solution of diffraction problems in the time-harmonic domain [ 15,17,18]. These are elliptic p.d.e. problems governed by the

Helmholtz equation and posed in an unbounded domain around an obstacle. In order to solve the Helmholtz equation by means of a numerical volume method (such as the finite element method or the finite difference method or a spectral method), it is necessary to truncate the unbounded domain up to an adequate artificial boundary. An appropriate condition on this boundary imposes, in an approximate way, that the outgoing waves are not reflected by the boundary. Therefore, these conditions are referred to as non-rejecting boundary conditions (NRBC) ([ 11,131 and the references herein). A high-order approximation of the exact, pseudodifferential NRBC can be very complicated (especially for the treatment of the bounding box corners) [9], whereas a simple, low-order approximation is effective only if the artificial boundary is at a quite large distance from the sources and the obstacles. As the oscillatory nature of the problems requires a great number of grid points per wavelength (typically 10 or 12 for the linear finite element method), the solution over a large 3D domain would demand the use of a prohibitive amount of computer memory. Furthermore, the stiffness matrix obtained after discretization is a sparse matrix, but it is complex, neither hermitian nor definite positive, because the Helmholtz operator is not definite positive or negative. Hence, the performances of the classical iterative methods would be, in this case, very unsatisfactory. A way to overcome these difficulties could be the decomposition of the domain into several non-overlapping subdomains, allowing the direct solution of a smaller problem in each subdomain [20]. This idea leads to an iterative method which converges to the solution of the problem if the solutions in the subdomains are suitably related by means of the boundary conditions at the interfaces. The domain decomposition method is, in its own nature, a parallel method.

* Corresponding

author. Present address:

CERFACS,

42 av. G. Coriolis,

0045-7825/98/$19.00 0 1998 Elsevier Science S.A. All rights reserved. PII: SOO4S-7825(97)00336-8

31057 Toulouse

Cedex, France.

114

A. Piacentini, N. Rosa I Comput. Methods Appl. Mech. Engrg. 162 (1998) 113-124

Recently, DesprCs [6-81 proposed a non-overlapping domain decomposition method in which the interface conditions are absorbing mixed conditions derived from the first-order NRBC. In this work the rate of convergence of the method will be improved by means of a new iteration scheme and of new interface conditions inspired by the second-order NRBC. The convergence of the new schemes is proved.

2. NRBCs

and domain decomposition

One way to obtain the exact NRBC over an infinite plane (say r which lets the variable orthogonal to the plane free to vary [9,15] ~(xl~~2,xg)-)~(5,~

:x3= 0)is the partial Fourier transform

5*,x,)

The NRBCs over finite surfaces are obtained The Helmholtz equation

as truncations

of the infinite one.

*u+
(1)

into an ordinary

differential

equation

w.r.t. x3

(2)

The solution in the half space x3 > 0 which propagates conditions on the plane x3 = 0

in the sense of x3 and which satisfies the continuity

is i( t,, t,, x3) = fi( t,, &,, 0+) e~(~Z~02’c*~“*X~ where the exponential is complex if 6’ < w2/c2. The derivative of the solution gives the operator A

-

E3(6,‘5*,0)=

(3)

linking

u and &l&c,

on the boundary

I’ : x3 = 0. From

(t*-$>“2P&

-

c

E*,O>

it follows g

r + T(+)

= 0

3

where

T is a non-local

pseudodifferential

operator

In the same way, it is possible to obtain the operator linking auldx, and u for the purely regressive wave. A NRBC can now be written imposing that the regressive wave vanishes on the artificial boundary. The same operators can be used in an interface condition for the domain decomposition method which imposes that the ingoing wave through an interface is equal to the outgoing wave from the subdomain which shares the common interface [4,19]. An elliptic problem of the form

A. Piacentini, N. Rosa I Comput. Methods Appl. Mech. Engrg. 162 (1998) 113-124

Au, + w2uk = -f,

in 0

a*, 8%

on the external

+ iwu, = 0

au, z + Tu, =

on

gk

k

115

boundary

an interface

has to be solved in each subdomain fik. It can be proved to be a well-posed problem as well as any approximation of it by means of a local variational condition. The proof makes use of classical functional analysis results and of the positiveness of T [10,22]. A local, variational condition can be obtained by means of a rational approximation of the symbol of the pseudodifferential operator T. Most of the energy is propagated by the low frequency (w.r.t. the space variables) components of the wave, so the approximation will be carried out for small values of [*, once the symbol of the operator has been rewritten in the form

(5) From now on, we assume, without loss of generality,

3. A first domain decomposition

c = 1.

scheme

In order to introduce the new domain decomposition method and to prove its convergence, we formerly recall in this section the method proposed by Despres [6] and those parts of the proof which are to be modified. Assuming a constant approximation of the square root 2

112

1-I

z1

(

cd2

)

Despres obtained the following scheme: Let {a,, . . . , a,} be a partition subdomains; then, given u”, for n = 1, 2, . . . solve for k = 1, . . . , K (Au;”

+ w’u;+’

= -f,

in LIk

+iwu;+’

=0

on c,,

of 0

into K non-overlapping

fl+1

auk

anext auk

+

iou;+’

= -an.

ank

+ iwuy

an

=

au,"

n+l

on Zik =

(6)

anjn aOk

I

The proof of the convergence involves the quantities e” = u - u” (where un is defined by uTn = ui) and e; = erflk. The linearity implies that e: satisfies the homogeneous Helmholtz problem, with the same boundary conditions. The key point of the demonstration is the following lemma (see also [6]). LEMMA 1. If aeilan, pseudo-energy

E L’(H&)

which satisfies the recursive ZTIa+’+a”*

i: /=I

PROOF.

The boundary

I

r,,,

Vk, then aeildn,

E L2(Ulk).

Therefore,

it is possible

to introduce

the

relation

le’l’ =

conditions

8’ of the iterative

(7) process satisfied by the error are

A. Piacentini, N. Rosa / Comput. Methods Appl. Mech. Engrg. 162 (1998) 113-124

116

-ae, = -iwe, It+1

?I+,

an,

ae,; - x

+ iwe]?

on Zkj

(8)

I

flfl de, II+1 = -iwe, 8%

on 4

(9)

Whence, aeLI&, E L*(dQ follows by induction, because ei E L*(&Q) V k, V n. We now split the boundary atik of each subdomain into the set of the interfaces Xkj and the external boundary 4, to get, by virtue of (9),

On the other hand, the interface

The last equality I,,

condition

relies on the technical

(8) on Zkj implies

lemma 4.1 in [6, p. 321, which states the isometry

l~-ioe12=Jlin~&+itie~*.

Since (9) yields

we get

The thesis readily follows after summation It follows from Lemma 1 that

over the iteration

index n.

lim n-F I r ex,lenI = 0 and

where we also used the continuous inclusion L2(dfilk) C Hp”*(dOk). Since any solution e E H’(O) of the Helmholtz equation Ae + w*e = 0 in 0 satisfies [3,5]

it follows that

(10)

A. Piacentini, N. Rosa I Comput. Methods Appl. Mech. Engrg. 162 (1998) 113-124

117

From now on we denote by C any constant independent of k, n. From ( 10) it is possible to prove the convergence of the domain decomposition scheme (6) starting from U: = 0 V k (see for all details [6, pp. 33-441).

4. A new iteration scheme We introduce two new iteration schemes that perform better than the original one. Following Despres [8], the interface condition in (6) can be relaxed: let y E (0, 1) be a relaxation parameter and let ?I+1

auk

+iwui+’ an,

= y(

-an+m. ““1

y)+(l-y)($$+iwU;)



(11)

The proof of the convergence remains unchanged. Otherwise, it is possible to use a Gauss-Seidel type scheme

au:*

n+ I

auk

-+iwu;+‘=--k+iou”’ an,

fork=l,...,K

I

I

(12)

withn*=n+l ifjk. In order to prove the convergence of the Gauss-Seidel type scheme we just need to show the analogous of Lemma 1. We start with a preliminary lemma. LEMMA

l,...,K

2. Vk=

(13) PROOF. Let us split the boundary a.C$ of each subdomain into the set of the interfaces skj and the external boundary 4, to get, by virtue of (10)

The conditions on 4 and on the interfaces ,Zjkimply (14) Forj
aeJ

I II zIL

%

I il 2

-

iwe,J+’

=

4k

de?+’

2

*

+

itie,?+’ - 2iwey+’

J

(15)

118

A. Piacentini, N. Rosa I Comput. Methods Appl. Mech. Engrg. 162 (1998) 113-124

The thesis easily follows by induction over k. The isometry lemma 4.1 in [6] is used in the following

relations:

(16) It is now possible

to prove the analogous

of Lemma

1.

LEMMA 3. The pseudo-energy

satisjes the recursive relation 8” + 4w2P$, 7

1 Ie,“j’ = 8’ J

(17)

It implies tle;+‘II~I,O~,6 C8” PROOF. Using

V k = 1,. . . , K

(18)

16) we get

whence, using (9) 8” + 4w2 C leJ’[’ = g’-’ i I$ Finally,

the proposition

(13) of Lemma 2 implies V k = 1, . . . , K

5. New interface operators It is possible to tune the operator introducing a balancing coefficient fi whose aim is to obtain the same rate of convergence to zero of the jumps across the interfaces [ulj, = ITr+ uj - TrZjk I+] and [au/&z] jk = lau,/an, + au,lan,(

A. Piacentini,

N. Rosa I Comput. Methods Appl. Mech. Engrg.

119

162 (1998) 113-124

1 au; i au;+’ -p ank +iwu[E+’ = ---+iwu,?

(19)

P anj

In the following, this method will be referred to as the ‘one parameter? method. A completely new operator can be obtained by means of the first-order approximation

of the symbol

i/2 +i” 2

(20)

cd2

By the inverse Fourier transform, P=Z--A,

this symbol corresponds

to the operator

1 202

in the space x,. x2, x3. where A, is the Laplace-Beltrami operator on I: The convergence of a scheme with an operator acting on u

au;

II+1

auk

+ iwFui+’

an,

= - an.

+ iwFuj” I

can still be proved if F : H ““(_Zj,) + H - “‘(_Zjk) is a definite positive operator. In this case F can be rewritten as S IL’+H-“~

F=SS*

and Lemma LEMMA

1 is modified

S*:H”‘+L2

as follows:

4. Let

(21) and

(22) Then, the pseudo-energy 2? = 8;. + q,, satisfies the recursive 8 “+‘+4W+ PROOF.

relation

j= I I rext

The interface

lej12 =

condition

+ iwS*ei+’

(23)

at the n + lth step of the process satisfied by the error is

= -S-l

2

+ ioS*ej” J

Whence

(24)

A. Piacentini, N. Rosu I Comput. Methods Appl. Mech. Engrg. I62 (1998) 113-124

The energy

gc,,

can be rewritten

The two identities

because

using the NRBC over 80

give the total energy

the quantity

Im >

vanishes

as the integral

has real value (see [6, Lemma 4.11). For the same reason, it is possible

to write

Therefore

The thesis follows by summation over the iteration index n. Therefore, the new interface condition will make use of a positive definite operator derived from P,with the introduction of a tuning coefficient a. A second tuning coefficient P, a relaxation parameter y and a Gauss-Seidel type iteration scheme will be employed as described in the previous sections. The new interface condition reads

+y In the following

(

u;*_K*fl;*+i’ cd2

au:*

po

anj

this method will be referred to as the ‘two parameters’

) method.

(25)

A. Piacentini, N. Rosa I Comput. Methods Appl. Mech. Engrg. 162 (1998) I13- 124

121

6. Variational formulation and discretization The three schemes can be written in the form Au + w2u = -f

in .Z2

Au+Bg=g

on r

(26)

where A = Z, B = -i/w Z on the external boundary and for the scheme proposed by Despres; A = Z, B = -i lpw Z for the second scheme (19) and A = Z - a/o* A,, B = -il@ Z for the third scheme (25). The right-hand side g is equal to 0 on the external boundary whereas it is calculated using the values of u and adan in the adjacent subdomains for the interface conditions. The variational formulation of the problem leads to a mixed finite element method [2 l] with the two unknown fields I( and r = au/an,,

in iI

(Ou.~~-W*UIF)dn-~~r~dT=In~~dn

tl+H’(0) (27)

Au$dT+

i

rBrt,bdT=

i

VlpsH-"2(r)

r g+dT

The operators A and B are variational

operators,

because

and the other operators are trivially variational. The integral over dZ can be neglected without loss of accuracy. The weak formulation was discretized using trilinear finite elements on a rectangular structured 3D grid for u and bilinear surface elements for r. No interpolation was needed for the discretization of the integrals over r, because the elements employed for r are the traces of the elements used for u. After discretization, a linear system is associated to each subdomain [14]. U

r

$2::

s:::lr:1

=[;I

If the operator A coincides with Z, then A,,, = AT,*. Anyway, the matrix A is always sparse and complex. The systems were solved by means of a direct Gaussian factorization with total pivoting implemented in the routine ME48 of the Harwell library. The matrices can be factorized only once because at any step of the iterative method only the right-hand side changes.

7. Numerical results The three different versions of the domain decomposition method were tested and compared over a large variety of cases. The curves shown here refer to the case of a cubic domain of side A with a pointwise source in the middle. The first two pictures refer to the effect of the new iteration scheme applied to the original method proposed by Despres. In Fig. 1 the convergence curves for the Jacobi scheme with or without relaxation are compared. In Fig. 2 the convergence curves for the Jacobi scheme and for the Gauss-Seidel scheme (both with optimal relaxation) are compared: the Gauss-Seidel scheme speeds up the convergence of a factor 2. It is interesting to notice that the optimal y is larger for the Gauss-Seidel scheme: this can be explained by the better accuracy of the right-hand side in this scheme. The next two pictures refer to the ‘one parameter’ method with A = Z and B = -iZl/30. In Fig. 3, the convergence of the original method with Gauss-Seidel scheme and relaxation and the new one are compared.

122

A. Piacentini,

N. Rosa

/ Comput. Methods Appl. Mech. Engrg. 162 (1998) 113-124

Fig. 1. Relaxation for the original method. Fig. 2. Gauss-Seidel vs. Jacobi for the original method.

Fig. 3. The ‘one parameter’ method. Fig. 4. ‘One parameter’ method: dependence on p.

The improvement is dramatic. Unfortunately, the efficiency of the new method is strictly related to the right choice of p and Fig. 4 shows the quite chaotic dependence of the number of iterations on the value of p. In the last pictures the ‘one parameter’ method is compared with the ‘two parameters’ one which has A = Z + (u/w* A,, B = --ill@. Fig. 5 shows a further improvement of the convergence. Moreover, the number of iterations depends smoothly on the parameters (Y and p and this is very important for applications. The optimal value of a confirms the prevision of l/2, while p seems to depend on the number of grid points per wavelength and on the distance between the source and the nearest interface, but these are only heuristic estimates (see Figs. 6 and 7). Finally, in Fig. 8, the method proposed by Despres (A = I, B = -ill@, Jacobi scheme, optimal relaxation) is compared with the ‘two parameters’ method. The original method achieves convergence in 10’02” on a Sparcstation 20 with 96Mb of RAM, while the new ‘two parameters’ operator leads the method to convergence in only 30”.

8. Conclusions

and perspectives

The dramatic improvement of the convergence speed due to the use of the new operator and scheme encourages the extension of the method to the cases of unstructured grids for complicated geometries and of

123

A. Piacentini, N. Rosa I Comput. Methods Appl. Mech. Engrg. 162 (1998) 113-124

%1

.

10

20

80

w

A

80

70

w

eo

WL * -1 -0.0

I

1w

4.8

4.7

-0.5

-0.8

-0.4

-0.3

-0.2

I

4.1

w

Fig. 5. The ‘two parameters’ Fig. 6. The ‘two parameters’

method.

method: dependence

on a.

--__

0.W

OS

I.

Et

0.4

0.45

a0

0.5

Fig. 7. The ‘two parameters’ Fig. 8. The ‘two parameters’

20

* 100

method: dependence

* 160

200

1

.

.

.

*

a2o200~40420500

on /3

method vs. the original

one.

heterogeneous media for which the speed of sound varies as C(X)=I/-. This last extension is already under study at the CERFACS, Toulouse, France, where the Despres method is discretized by means of a hybrid mixed finite element formulation [lo]. Furthermore, a theorethical study on how the convergence depends on the parameters (Y and p is an actual object of research.

Acknowledgments The authors would like to thank I? Joly from JNRIA-La CERFACS for their precious teaching and advice.

Rocquencourte

and F. Collino

from INRIA and

References [I] P.E. Bjorstad and O.B. Widlund, To overlap or not to overlap: a note on domain decomposition Sci. Stat. Comput. 10(.5) (1989) 1053-1061.

methods for elliptic problems,

SIAM J.

124

A. Piacentini,

N. Rosa I Comput. Methods Appl. Mech. Engrg. 162 (1998) 113-124

121 J. Bramble, J.E. Pas&k and H. Schatz, The construction of preconditioners for elliptic problems by substructuring, Math. Comput. 46 (1986) 361-369. [3] H. Brezis, Analyse Fonctionelle (Masson, Paris, 1983). 141 F. Collino and P. Joly, New absorbing boundary conditions for the finite element solution of 3D Maxwell’s equations, IEEE Trans. Magnetics 31(3) (Ptl) (1995) 1696-1701. [5] R. Dautray and J.L. Lions, Analyse Mathematique et Calcul Numerique (Masson, Paris, 1985). [6] B. Despres, Methodes de decomposition de domaine pour les probltmes de propagation d’ondes en regime harmonique, Ph.D. Thesis, Paris 9, 1991. [7] B. Despres, P. Joly and E. Roberts, A domain decomposition method for the harmonic Maxwell equation, in: Iterative Methods in Linear Algebra (Elsevier, Amsterdam, 1992) 475-484. [8] B. Despres, Domain decomposition method and the Helmholtz problem, in: Mathematical and Numerical Aspects of Wave Propagation Phenomena (SIAM, Philadelphia, 1991) 44-52. [9] B. Engquist and A. Majda, Radiation boundary conditions for acoustic and elastic wave propagation, Comm. Pure Appl. Math. 32 (1979) 313-357. [lo] S. Ghanemi, Methode de decomposition de domaine avec condition de transmission non locale, Ph.D. Thesis, Paris 9, 1995. [ 1 l] D. Givoli, Numerical Methods for Problems in Infinite Domain (Elsevier Science Publishers, Amsterdam, 1992). [12] R. Glowinski and M. Wheeler, Domain decomposition methods and mixed finite element methods, in: First International Symposium on Domain Decomposition Methods for Partial Differential Equations (SIAM, Philadelphia, 1988) 144-171. [l3] L. Halpem and J. Rauch, Absorbing boundary conditions for diffusion equations, Numer. Math. 71(2) (1995) 185-224. [14] C. Johnson, Numerical Solution of Partial Derivative Equations by the Finite Element Method (Cambridge University Press, Cambridge, 1987). [IS] P. Joly et al., Les bases de la theorie de la diffraction, in: Ecole d’Aquitaine: Ondes et matiire, 26-30 settembre 1994 Maubuisson, F; tea-cesta, Bordeaux, 1994. [ 161 A. Kufner, 0. John and S. FuEik, Function Spaces (Noordhoff International Publishing, Leyden, 1977). [17] L.D. Landau and E.M. Lifshits, Fluid Mechanics (Pergamon Press, Oxford, 1959). [ 181 Z. Maekawa and P. Lord, Environmental and Architectural Acoustics (Fn Spon, London, 1994). [l9] F. Nataf, F. Rogier and E. de Sturler, Optimal interface conditions for domain decomposition methods, Ecole pdytechnique Internal Report 301, April 1994. [20] A. Quarteroni et al. (eds.), Domain Decomposition Methods in Science and Engineering, June 15-19, Como, Italy (Rhode Island, Providence, 1994). (211 J.E. Roberts and J.M. Thomas, Mixed and hybrid methods, in: Handbook of Numerical Analysis,Vol. 3 (North Holland, Amsterdam). [22] L.N. Trefethen and L. Ha&m, Well-posedness of the one-way wave equations and absorbing boundary conditions, Math. Comput. 47 (1986) 421-435.