PROCEEDINGS of the23rd International Congress on Acoustics
9 to 13 September 2019 in Aachen, Germany

A domain decomposition method with fast convergence for theHelmholtz equation
Denis DUHAMEL1
1

Ecole des Ponts ParisTech, FranceABSTRACT
Solving the Helmholtz equation by finite element methods is quite important in acoustics. When thefrequency or the size of the problem increase, large meshes are necessary and consequently heavy
computations are required. One possibility is to use domain decompositions for which the domain isdecomposed into subdomains on which the solutions can be computed more easily. This involves aniterative scheme where data are transmitted between subdomains from the precedent iteration. The main
problem is to have a low number of iterations so that the problem can be solved in a reasonable amount oftime.
In this work, we present a domain decomposition method based on two main features. The first one is touse extended domains with absorbing boundary conditions. The second feature is to decompose the wholedomain into one-dimensional or two-dimensional networks of subdomains so that double sweep
preconditioners can be used. Examples are shown where the number of iterations is usually low. This
number of iterations is also shown to depend slowly on the number of domains and the frequency.
Keywords: Domain decomposition, Finite Element, Helmholtz equation.

1. INTRODUCTION
We consider here the solution of the Helmholtz equation by finite element methods. When thefrequency increases many degrees of freedom are needed. In such cases, it can be interesting to use domain
decomposition methods for which the domain is decomposed into several subdomains over which thesolutions can be computed more easily by solving independent small size problems. Then an iterative
scheme needs to be developed to compute the solutions in the subdomains at a given step from thesolutions over the other subdomains at the precedent step. The problem is to get as few iterations as
possible and a number of iterations almost independent of the number of subdomains and the frequency.
Many works, such that the development of FETI methods or Schwartz methods, have been done in thepast on this subject, see (1-9) for instance. In overlapping Schwartz methods (3), the solutions arecomputed on overlapping domains and the boundary conditions for a given domain are defined at theboundary of this overlapping domain and so from the solutions at the precedent iteration obtained in theother domains that the present domain overlaps. In non-overlapping Schwartz methods (4, 5), the solutions
are computed in each domain and transmission operators are defined at the boundary of the domains to giveboundary conditions at a given iteration from the solutions of adjacent domains at the precedent iteration.
The number of iterations to achieve the convergence of the algorithm is deeply dependent on the quality ofthese transmission operators.
Different methods have been proposed to compute these transmission operators. The best one are
non-local such as the Dirichlet to Neumann operator but are too complex to be really useful (6). Different
local operators have been proposed, see for instance (7) and (8). A good one is the Perfectly Matched Layer(PML) as in (9) but it needs to define special layers around each domain with suitable absorption
properties.
Another possibility, which is developed here, is to define the transmission operator from the solutions ofthe Helmholtz equation in the subdomains completed with surrounding subdomains and with simple
absorbing conditions at the exterior boundary. Then the transmission operator is obtained from the solutionsin the surrounding subdomains. For being efficient, a preconditionner must be used with these methods
such as the double sweep preconditionner found in (10, 11,12). This method is developed for the cases
17496

denis.duhamel@enpc.fr

where the subdomain decomposition is made of one-dimensional or two-dimensional networks so that theglobal domain is first decomposed into a sequel of slices at the first level and then each slice can also bedecomposed into a second level of subslices.
The paper is organized as follows. In section 2, the method is presented for the Helmholtz equation
defined on a two-dimensional domain. Then in section 3, numerical results are given before a conclusion.

2. COMPUTATION BY DOMAIN DECOMPOSITION
2.1

General problem
The objective is to solve by domain decomposition, the Helmholtz equation given by(1)ο‫ ݌‬൅ ݇ ଶ ‫ = ݌‬0with p the pressure, k=ω/c the wavenumber, c the sound velocity, ω=2πf and f the frequency. This is definedడ௣
on a two-dimensional domain Ω with a Neumann boundary condition డ௡ ൌ ‫ ݍ‬on the exterior boundary Γof Ω.
Absorbing boundaries: i-1:i: i+ii-1

ii+1**i-2,i-1*i-1,i*i,i+1*i+1,i+2

Figure 1: Domain decomposed into one-dimensional slices
The classical way of solving the problem by domain decomposition is to divide the global domain into
subdomains, for instance into slices like in Figure 1, and to solve the sequel of following problems indexed
by n for each subdomain ȳ௜ ο‫݌‬௜௡ ൅ ݇ ଶ ‫݌‬௜௡ ൌ Ͳ݅݊ȳ௜߲‫݌‬௜௡௡ିଵ൅ ܵ‫݌‬௜௡ ൌ ݃௜௝‫݊݋‬Ȟ௜௝
(2)߲݊௜߲‫݌‬௜௡ൌ ‫ݍ‬௜ ‫݊݋‬Ȟ௜߲݊௜With Ȟ௜௝ the boundary between the domains i and j and S is an operator that try to simulate theimpedance of the domain that extends beyond the boundary. The values of g are updated by߲‫݌‬௝௡൅ ܵ‫݌‬௝௡ ‫݊݋‬Ȟ௜௝߲݊௜௡ିଵൌ െ݃௝௜൅ ʹܵ‫݌‬௝௡ ௡݃௜௝
ൌെ

(3)

7497

Decomposing ‫݌‬௜௡ ൌ  ‫ݒ‬௜௡ ൅ ‫ݓ‬௜ with ‫ݒ‬௜௡ the solution of problem (1) with ‫ݍ‬௜ ൌ Ͳ and ‫ݓ‬௜ the solution
௡ିଵwith ݃௜௝= 0, the problem can be put under the form݃௡ ൌ ‫݃ܣ‬௡ିଵ ൅ ܾ(4)௡With the vector ܾ obtained from ‫ݓ‬௜ , the matrix A from ‫ݒ‬௜௡ and ݃௡ is the vector of all the ݃௜௝fordifferent i and j. The problem (3) is then solved by the GMRES algorithm associated to the double sweep
preconditionner described in (10, 11, 12). The main problem consists in finding a good operator S such thatthe number of iterations to get the convergence is as low as possible and with a low dependency on thenumber of subdomains and the frequency.

2.2

Transmission operator
To get this operator S, one proposes to solve on extended domains by adding to the domain ȳ௜ thedomains ȳ௜ିଵ ȳ௜ାଵ associated to simple absorbing boundary conditions to the exterior of the global
domain, that is on the boundaries Ȟ௜ିଶǡ௜ିଵ and Ȟ௜ାଵǡ௜ାଶ, see Figure 1. So the problem (1) is transformedintoο‫݌‬௜௡ ൅ ݇ ଶ ‫݌‬௜௡ ൌ Ͳ݅݊ȳ௜ିଵ ‫ ׫‬ȳ௜ ‫ ׫‬ȳ௜ାଵ߲‫݌‬௜௡െ ݅݇‫݌‬௜௡ ൌ Ͳ‫݊݋‬Ȟ௜ିଶǡ௜ିଵ߲݊௜ିଵ߲‫݌‬௜௡െ ݅݇‫݌‬௜௡ ൌ Ͳ‫݊݋‬Ȟ௜ାଵǡ௜ାଶ߲݊௜ାଵ߲‫݌‬௜௡߲‫݌‬௜௡௡ିଵ൅ൌ ݃௜ǡ௜ିଵ‫݊݋‬Ȟ௜ିଵǡ௜߲݊௜ିଵ ߲݊௜߲‫݌‬௜௡߲‫݌‬௜௡௡ିଵ൅ൌ ݃௜ǡ௜ାଵ‫݊݋‬Ȟ௜ǡ௜ାଵ߲݊௜ ߲݊௜ାଵ߲‫݌‬௜௡ൌ Ͳ‫݊݋‬Ȟ௜ିଵ ‫ ׫‬Ȟ௜ ‫ ׫‬Ȟ௜ାଵ߲݊௜(5)௡ିଵ௡ିଵ
Where ݃௜ǡ௜ିଵand ݃௜ǡ௜ାଵare values obtained from the precedent iteration by relation (3) and the operator S isobtained as the normal derivative of the solution of the pressure on the extended domains ȳ௜ିଵ ȳ௜ାଵ.
More precisely on ȳ௜ିଵ one has ܵ‫ ݑ‬ൌడ௨డ௡೔షభ

on Ȟ௜ିଵǡ௜ with u the solution of the following problemο‫ ݑ‬൅ ݇ ଶ ‫ ݑ‬ൌ Ͳ݅݊ȳ௜ିଵ߲‫ݑ‬െ ݅݇‫ ݑ‬ൌ Ͳ‫݊݋‬Ȟ௜ିଶǡ௜ିଵ߲݊௜ିଵ߲‫ݑ‬ൌ Ͳ‫݊݋‬Ȟ௜ିଵ߲݊௜ିଵ‫ ݑ‬ൌ ‫݌‬௜௡ ‫݊݋‬Ȟ௜ିଵǡ௜(6)

And a similar definition for the domain ȳ௜ାଵ.
2.3

Two-dimensional network
The problem can be solved as in the precedent subsection or can be further decomposed into a

two-dimensional network of subdomains as in Figure 2. In this latter case, the sequel of following problems
indexed by m is solved.

7498

௡ǡ௠
௡ǡ௠
෩ ௜ǡ௝
ο‫݌‬௜௝൅ ݇ ଶ ‫݌‬௜௝ൌ Ͳ݅݊ȳ೙ǡ೘డ௣೔ೕడ௡೙ǡ೘డ௣೔ೕడ௡ೕషభ೙ǡ೘డ௣೔ೕడ௡ೕ௡ǡ௠െ ݅݇‫݌‬௜௝ൌ Ͳ‫݊݋‬Ȟ෨௜ǡ௝(7)డ௣೙ǡ೘௠ିଵ൅ డ௡೔ೕ ൌ ‫ݎ‬௝ǡ௝ିଵ‫݊݋‬Ȟ௝ିଵǡ௝ೕ൅೙ǡ೘డ௣೔ೕడ௡ೕశభ௠ିଵൌ ‫ݎ‬௝ǡ௝ାଵ‫݊݋‬Ȟ௝ǡ௝ାଵ

Absorbing boundaries

j+1: ij

j:i+1j-1

j-1

*i-1

i

i+1

Figure 2: Two-dimensional network of subdomains෩ ࢏ǡ࢐ is made of the domain ȳ௜ǡ௝ and its surrounding domains such thatWhereȳ෩ ௜ǡ௝ =ȳ௜ǡ௝ ‫ ׫‬ȳ௜ିଵǡ௝ ‫ ׫‬ȳ௜ାଵǡ௝ ‫ ׫‬ȳ௜ିଵǡ௝ିଵ ‫ ׫‬ȳ௜ǡ௝ିଵ ‫ ׫‬ȳ௜ାଵǡ௝ିଵ ‫ ׫‬ȳ௜ିଵǡ௝ାଵ ‫ ׫‬ȳ௜ǡ௝ାଵ ‫ ׫‬ȳ௜ାଵǡ௝ାଵȳ(8)෩ ௜ǡ௝ . The iteration continues on m until convergence in slice i toAnd Ȟ෨௜ǡ௝ is the exterior boundary of ȳ
get the solution of problem (5). Other points such as the definition of the transmission operator in slice i aredefined in a similar way to relations (2) and (6).

3. NUMERICAL EXAMPLES
3.1 One-dimensional decomposition
We consider the Helmholtz equation on a square domain of size 1m×1m. The sound velocity is
c=340m/s. The domain is broken down into vertical slices as in Figure 3. First the boundary condition
is defined by a plane wave such that the normal derivative of the pressure is ‫ ݍ‬ൌ ݅݇݊௫ ݁ ௜௞௫ . The numberof iterations is given in Table 1 and an example of solutions for the frequency 1000Hz is presented inFigure 4.
Table 1: Number of iterations for different frequencies f and number of subdomains n for a plane wavef݊ൌͷ݊ ൌ ʹͷ݊ ൌ ͷͲ݊ ൌ ͳͲͲ
500Hz
55551000Hz
55552000Hz
5955In this example, one can see that the number of iterations is low and almost independent of the7499

frequency and the number of subdomains. Figure 4 shows the plane wave on the whole domain withfull continuity at the boundaries of the subdomains.

Figure 3: Domain decomposition into 25 subdomains

Figure 4: Pressure for a plane wave at 1000Hz
Nextthecaseସగ௬‫ݍ‬ሺ‫ݔ‬ǡ ‫ݕ‬ሻ ൌ ݅݇௫ ݊௫  ቀofaସగସగ௬higherorderwavedefinedbyସగ ଶቁ ݁ ௜௞ೣ ௫ െ ௅ ݊௬  ቀ ௅ ቁ ݁ ௜௞ೣ ௫ with ݇ ଶ ൌ ݇௫ଶ ൅ ቀ ௅ ቁ and L=1m is presented.
௅The number of iterations is given in Table 2 and an example of results for the frequency 1000Hz ispresented in Figure 4.
Table 2: Number of iterations for different frequencies f and number of subdomains n for a higher modef݊ൌͷ݊ ൌ ʹͷ݊ ൌ ͷͲ݊ ൌ ͳͲͲ500Hz
41012131000Hz
71010107500

2000Hz

613713In this case, the number of iterations depends on the frequency and the number of subdomains.
However, it seems to reach rapidly a maximum value still independent of the frequency and the number ofsubdomains. Figure 5 presents the pressure for the frequency 1000Hz. A good continuous solution can beseen.

Figure 5: Pressure for the fourth mode at 1000Hz
3.2

Two-dimensional decomposition
We first consider the case of a rectangular domain of size 1m×1m with a boundary condition given by aplane wave eikx on the left and right boundaries as in the precedent subsections. The domain is now
divided into ݊ ൈ ݊ subdomains and each subdomain is mesh with triangular elements of degree 1. Anexample of decomposition is presented in Figure 6. Each subdomain has approximately 400 elements sothat the total number of elements will depend on the number of subdomains. Table 3 presents the number of
iterations for different frequencies and number of subdomains. The number inside the parenthesis is theaverage number of subiterations for solving problem (7) into a slice.
Table 3: Number of iterations for different frequencies and number of subdomainsfͷൈͷͳͲ ൈ ͳͲʹͲ ൈ ʹͲͶͲ ൈ ͶͲ500Hz
5 (8.3)4 (8.2)4 (9.9)4 (11.0)
1000Hz5 (8.5)4 (8.2)4 (8.9)5 (9.4)
2000Hz
9 (12.6)5 (8.4)6 (8.0)4 (8.2)
One can see that the numbers of iterations is low and with a low dependency on the frequency or thenumber of subdomains. Almost the same observation can be made for the average number of subiterations
in a slice. Figure 7 presents the equivalent of figure 5 for the same boundary conditions. The solution is thesame even if in the present case the computation is made on the two-dimensional network.

7501

yx
Figure 6: Domain decomposition into 10x10 subdomains

Figure 7: Pressure for the fourth mode at 1000Hz for 10x10 subdomains

4. CONCLUSION
We have presented transmission conditions which are intermediate between the local conditions on theboundary and domain conditions as provided by a PML layer. So it gives an equilibrium between the easeof implementation and the efficiency. It involves only modified matrices on the boundaries but no newmatrix in the domains. This was shown to converge in a low number of iterations both for one-dimensional7502

and two-dimensional networks of subdomains. Subsequent research should focus on the improvement of
the boundary condition and the consideration of domains with more complex shapes.

REFERENCES
1. Farhat C., Roux F.X., An unconventional domain decomposition method for an efficient parallel
solution of large-scale finite element systems, SIAM J. Sci. Stat. Comput., 13, 379-396, 1992.
2. Farhat C., Avery P., Tezaur R., Li J., FETI-DPH: A dual-primal domain decomposition method foracoustic scattering, Journal of Computational Acoustics, 13, 499-524, 2005
3. Gander M.J., Zhang H., Optimized Schwarz Methods with Overlap for the Helmholtz Equation, SIAM
Journal on Scientific Computing, 38, No. 5, A3195-A3219, 2016.
4. Gander M., Optimized Schwarz methods, SIAM J. Numer. Anal., 44, 699-731, 2006.
5. Gander M., Magoules F., Natar F., Optimized Schwarz methods without overlap for the Helmholtz
equation, SIAM J. Sci. Comput., 24 38-60, 2002.
6. Nataf F., Interface connections in domain decomposition methods, In Modern methods in scientificcomputing and applications (Montréal, QC, 2001), volume 75 of NATO Sci. Ser. II Math. Phys. Chem.,
pages 323-364. KluwerAcad. Publ., Dordrecht, 2002.
7. Thierry B., Vion A., Tournier S., El Bouajaji M., Colignon D., Marsic N., Antoine X., Geuzaine C.,
GetDDM: An open framework for testing optimized Schwarz methods for time-harmonic wave
problems, Computer Physics Communications, 203, 309-330, 2016
8. Boubendir Y., Antoine X., Geuzaine C., A quasi-optimal non-overlapping domain decompositionalgorithm for the Helmholtz equation, Journal of Computational Physics, 231, 262-280, 2012.
9. Stolk C.C., A rapidly converging domain decomposition method for the Helmholtz equation, Journalof Computational Physics, 241, 240-252, 2013.
10. Vion A., Geuzaine C., Double sweep preconditioner for optimized Schwarz methods applied to theHelmholtz problem, Journal of Computational Physics, 266, 171-190, 2014.
11. Eslaminia M., Guddati M.N., A double-sweeping preconditioner for the Helmholtz equation, Journalof Computational Physics, 314, 800-823, 2016.
12. Vion A., Bélanger-Rioux R., Demanet L., Geuzaine C., A DDM double sweep preconditioner for theHelmholtz equation with matrix probing of the DtN map, Proceedings of the 11th International
Conference on Mathematical and Numerical Aspects of Waves (WAVES 2013), Tunis, Tunisia, June 3-7,
2013.

7503