A Multigrid Method Enhanced by Krylov Subspace Iteration for Discrete Helmholtz Equations.
Standard multigrid algorithms have proven ineffective for the solution of discretizations of Helmholtz equations. In this work we modify the standard algorithm by adding GMRES iterations at coarse levels and as an outer iteration. We demonstrate the algorithm's effectiveness through theoretical analysis of a model problem and experimental results. In particular, we show that the combined use of GMRES as a smoother and outer iteration produces an algorithm whose performance depends relatively mildly on wave number and is robust for normalized wave numbers as large as 200. For fixed wave numbers, it displays grid-independent convergence rates and has costs proportional to the number of unknowns.
Introduction
. Multigrid algorithms are effective for the numerical solution
of many partial differential equations, providing a solution in time proportional to
the number of unknowns. For some important classes of problems, however, standard
multigrid algorithms have not been useful, and in this paper we focus on developing
effective multigrid algorithms for one such class, the discrete Helmholtz equation.
Our main interest lies in solving exterior boundary value problems of the form
on\Omega ae R d (1.1)
@\Omega (1.2)
such as arise in the modeling of time-harmonic acoustic or plane-polarized electromagnetic
scattering by an obstacle. The boundary \Gamma represents the scattering obstacle,
and the boundary operator B can be chosen so that a Dirichlet, Neumann or Robin
boundary condition is imposed. The original unbounded domain is truncated to the
finite
domain\Omega by introducing the artificial boundary \Gamma 1 on which the radiation
boundary condition (1.3) approximates the outgoing Sommerfeld radiation condition.
Depending on what type of radiation condition is chosen, M can be either a (local)
differential operator or a global integral operator coupling all points on \Gamma 1 (see [14]).
The data for the problem are given by the right hand side f and the boundary data
g. In the most common case, f j 0 and \Gammag is the boundary data of an incident
plane wave. The critical parameter is the wave number k, which is positive in the
case of unattenuated wave propagation. Due to the radiation boundary condition,
the solution of (1.1)-(1.3) is a complex-valued function u
Discretization of (1.1)-(1.3) by finite differences or finite elements leads to a linear
system of equations
in which the coefficient matrix A is complex-symmetric, i.e., not Hermitian. Moreover,
for large values of the wave number k, it becomes highly indefinite.
It is this indefiniteness that until recently has prevented multigrid methods from
being applied to the solution of the discrete equations with the same success as these
methods have enjoyed for symmetric positive definite problems. As will be illustrated
in Section 2, the difficulties with standard multigrid methods applied to Helmholtz
problems concern both of the main multigrid components: smoothing and coarse grid
correction. In particular, standard smoothers such as Jacobi or Gau-Seidel relaxation
become unstable for indefinite problems since there are always error components-
usually the smooth ones-which are amplified by these smoothers. The difficulties
with the coarse grid correction are usually attributed to the poor approximation of
the Helmholtz operator on very coarse meshes, since such meshes cannot adequately
resolve waves with wavelength of which the solution primarily consists. We
show, however, that although the coarse grid correction is inaccurate when coarse-grid
eigenvalues do not agree well with their fine-grid counterparts, coarse meshes can still
yield useful information in a multigrid cycle.
In this paper, we analyze and test techniques designed to address the difficulties
in both smoothing and coarse grid correction for the Helmholtz equation. For
smoothing, our approach is to use a standard, damped Jacobi, relaxation when it
works reasonably well (on fine enough grids), and then to replace this with a Krylov
subspace iteration when it fails as a smoother. Earlier work such as in Bank [2]
and Brandt and Ta'asan [11] have employed relaxation on the normal equations in
this context. Krylov subspace smoothing, principally using the conjugate gradient
method, has been considered by a variety of authors [3, 7, 8, 27, 29].
For coarse grid correction, we identify the type and number of eigenvalues that
are handled poorly during the correction, and remedy the difficulty by introducing an
outer acceleration for multigrid; that is, we use multigrid as a preconditioner for an
outer, Krylov subspace, iteration. This approach has been used by many authors, e.g.
[28, 31] but only for problems in which the coarse grid is restricted to be fairly fine.
It has also been used in other settings [20, 23]. Any Krylov subspace method is an
option for both the smoother and the outer iteration; we use GMRES [25]. In contrast
to many multilevel strategies [2, 6, 9, 31], the resulting algorithm has no requirements
that the coarse grid be sufficiently fine. For approaches based on preconditioning
indefinite problems by preconditioners for the leading term, see [4, 5, 15, 32].
In more recent work, Brandt and Livshits [10] have developed an effective multi-grid
approach for the Helmholtz equation based on representing oscillatory error components
on coarse grids as the product of an oscillatory Fourier mode and a smooth
amplitude-or ray-function. The standard V-cycle is augmented by so-called ray cy-
cles, in which the oscillatory error components are eliminated by approximating the
associated ray functions in a multigrid fashion. This wave-ray methodology has also
been combined by Lee et al. [21] with a first-order system least-squares formulation for
the Helmholtz equation. These approaches require construction of and bookkeeping
for extra grids associated with the ray functions.
An outline of the paper is as follows. In Section 2, we perform a model problem
analysis, using a one-dimensional problem to identify the difficulties encountered by
both smoothers and coarse grid correction, and supplementing these observations
with an analysis of how dimensionality of the problem affects the computations. In
Section 3, we present the refined multigrid algorithms and test their performance on
a set of two-dimensional benchmark problems on a square domain. In particular, we
demonstrate the effectiveness of an automated stopping criterion for use with GMRES
smoothing, and we show that the combined use of GMRES as a smoother and outer
iteration produces an algorithm whose performance depends relatively mildly on wave
number and is robust for wave numbers as large as two hundred. In Section 4, we show
the performance of the multigrid solver on an exterior scattering problem. Finally, in
Section 5, we draw some conclusions.
2. Model Problem Analysis. Most of the deficiencies of standard multigrid
methods for solving Helmholtz problems can be seen from a one-dimensional model
problem. Therefore, we consider the Helmholtz equation on the unit interval (0; 1)
with homogeneous Dirichlet boundary conditions
This problem is guaranteed to be nonsingular only if k 2 is not an eigenvalue of the
negative Laplacian, and we will assume here that this requirement holds. The problem
is indefinite for which is the smallest eigenvalue of the negative Laplacian.
Finite difference discretization of (2.1) on a uniform grid containing N interior
points leads to a linear system of equations (1.4) with the N \Theta N coefficient matrix
denotes the mesh
width and I denotes the identity matrix. Under the assumptions on k above, it is
well-known (see [26]) that for sufficiently fine discretizations, the discrete problems
are also nonsingular. We also assume that all coarse grid problems are nonsingular.
The eigenvalues of A are
and the eigenvectors are
2h [sin ij-h] N
The choice of Dirichlet boundary conditions in (2.1) allows us to perform Fourier
analysis using these analytic expressions for eigenvalues and eigenvectors. In experiments
described in Section 3, we will examine how our observations coincide with
performance on problems with radiation conditions, which are nonsingular for all k
[18]. Aspects of the algorithm that depend on the dimensionality of the problem will
be considered at the end of this section.
2.1. Smoothing. For the smoothing operator, we consider damped Jacobi re-
laxation, defined by the stationary iteration
em um denote the residual and error vectors at
step m, respectively. denotes the matrix consisting of the diagonal
of A, and ! is the damping parameter. The associated error propagation matrix
is and the eigenstructure of this matrix governs the behavior
of the error em . Since D is a multiple of the identity matrix, S! is a
polynomial in A and hence shares the same system of orthonormal eigenvectors (2.3).
The eigenvalues of S! are
Thus, the eigenvalue - j of S! is the damping factor for the error component corresponding
to the eigenvalue - j of A.
We now consider the effects of damped Jacobi smoothing on three levels of grids:
fine, coarse, and intermediate.
2.1.1. Fine Grids. The fine grid mesh size is determined by accuracy requirements
on the discretization, and this allows us to make certain assumptions on the
size of h versus k on the fine grid. Recall that the wavelength - associated with a
time-harmonic wave with wave number k ? 0 is given by 2-=k. The quantity
is the number of mesh points per wavelength, and it measures the approximability of
the solution on a given mesh. A commonly employed engineering rule of thumb [17]
states that, for a second-order finite difference or linear finite element discretization,
equivalently, kh -=5 (2.5)
is required, and we will enforce (2.5) in all experiments. We also note that, for reasons
of stability, a bound on the quantity h 2 k 3 is also required [18]; for high wave numbers
this bound is more restrictive than the bound on kh.
As a consequence of (2.5), the quantity multiplying the smoothing parameter !
in (2.4) will vary between about \Gamma1=4 and 9=4 for
smoothing results in a slight amplification of the most oscillatory modes as
well as of the smoothest modes. One can adjust ! so that the most oscillatory mode
is damped, and this is the case as long as For S! to
be an effective smoother, ! is usually chosen to maximize damping for the oscillatory
half of the spectrum. This leads to the choice
which is equal to the familiar optimal value of 2=3 for the Laplacian [22, p. 11] when
But the smoothest mode is amplified for any
positive choice of ! when the discrete problem is indefinite, and this is the case for
the discrete Helmholtz operator . As can be seen from (2.4), more
smooth-mode eigenvalues of S! become larger than one in magnitude as h is increased,
thus making damped Jacobi-as well as other standard smoothers-increasingly more
unstable as the mesh is coarsened.

Figure

2.1 shows the damping factors - j for each of the eigenvalues - j of A for
wave number on a grid with 31. The maximal amplification occurs for
the smoothest mode, corresponding to the leftmost eigenvalue of A. When
this amplification factor is approximately equal to

Figure

2.2 shows how ae varies with kh. Limiting this largest amplification factor, say
to ae - 1:1, would lead to the mesh size restriction kh - 0:52, somewhat stronger than
(2.5). One also observes that, for kh ?
6, this mode is once again damped.
In summary, the situation on the finest grids is similar to the positive definite
case, except for the small number of amplified smooth modes whose number and
amplification factors increase as the mesh is coarsened.
l
Fig. 2.1. The damping factors for the damped Jacobi relaxation plotted against the eigenvalues
of A (+) for
kh
Fig. 2.2. The variation of the damping/amplification factor of the smoothest mode as a function
of kh for
2.1.2. Very Coarse Grids. As the mesh is coarsened, the eigenvalues of A that
correspond to the larger eigenvalues of the underlying differential operator disappear
from the discrete problem, while the small ones-those with smooth eigenfunctions-
remain. This means that, for a fixed k large enough for the differential problem to
be indefinite, there is a coarsening level below which all eigenvalues are negative. For
the model problem (2.1), this occurs for kh ? 2 cos(-h=2) for any fixed k ? -. In this
(negative definite) case, the damped Jacobi iteration is convergent for
with
and the spectral radius of S! is minimized for This would permit the use of
(undamped) Jacobi as a smoother on very coarse grids, but we shall not make use of
this.
2.1.3. Intermediate Grids. What remains is the difficult case: values of kh
for which the problem is not yet negative definite but for which a large number of
smooth modes are amplified by damped Jacobi relaxation. Jacobi smoothing and
other standard smoothers are therefore no longer suited, and it becomes necessary to
use a different smoothing procedure. In [11] and [16] it was proposed to replace classical
smoothers with the Kaczmarz iteration, which is Gau-Seidel relaxation applied to
the symmetric positive-definite system AA   for the auxiliary variable v defined
by A   This method has the advantage of not amplifying any modes, but it
suffers from the drawback that the damping of the oscillatory modes is very weak. In
the following section we propose using Krylov subspace methods such as GMRES for
smoothing. These methods possess the advantage of reducing error components on
both sides of the imaginary axis without resorting to the normal equations.
2.2. Coarse Grid Correction. The rationale behind coarse-grid correction is
that smooth error components can be well represented on coarser grids, and hence
a sufficiently good approximation of the error can be obtained by approximating
the fine grid residual equation using the analogous system on a coarser mesh. This
assumes both that the error consists mainly of smooth modes and that the solution
of the coarse grid residual equation is close to its counterpart on the fine grid. In this
section, we present an analysis of what goes wrong for the Helmholtz problem.
2.2.1. Amplification of Certain Modes. Assume the number of interior grid
points on the fine grid is odd, and consider the next coarser mesh, with
interior points. We identify R N and R n , respectively, with the spaces of grid functions
on these two meshes that vanish at the endpoints, and we indicate the mesh such
vectors are associated with using the superscripts h and H. Let e
the fine grid error, let r denote the residual, and let denote the
coarse mesh size. Let the coarse-to-fine transformation be given by the interpolation
operator I h
\Theta I h

The following indication of what can go wrong with the (exact) coarse grid correction
was given in [11]: consider a fine-grid error e consisting of only the
smoothest eigenvector v h of A h with associated eigenvalue - h . The fine-grid residual
is thus given by r since we are assuming that v h is smooth,
its restriction - r H := I H
h to the coarse grid will again be close to an
eigenvector of the coarse-grid operator A H , but with respect to a slightly different
eigenvalue - H . The coarse grid version of the correction is
I H
Hence the error on the fine grid after the correction is
where we have assumed that the smooth mode v h is invariant under restriction followed
by interpolation. This tells us that, under the assumption that the restrictions
of smooth eigenvectors are again eigenvectors of A H , the quality of the correction
depends on the ratio - h =- H . If the two are equal, then the correction is perfect,
but if the relative error is large, the correction can be arbitrarily bad. This occurs
whenever one of - h , - H is close to the origin and the other is not. Moreover, if - h
and - H have opposite signs, then the correction is in the wrong direction.
We now go beyond existing analysis and examine which eigenvalues are problematic
in this sense for finite differences; a similar analysis can also be performed for
linear finite elements. Consider the coarse-grid eigenfunctions v H
. To
understand the effects of interpolation of these grid functions to the fine grid, we must
examine both the first n fine-grid eigenfunctions fv h
and their complementary
modes fv h
are related by
\Theta v h

\Theta v h

. As is easily
verified, there holds [12]
I h
with c j := cos j-h=2 and s j := sin j-h=2,
If full weighting is used for the restriction operator I H
componentwise
\Theta
I H


\Theta

\Theta

and the relation I h
. The following mapping properties are easily established

I H
with c j and s j as defined above.
If A H denotes the coarse-grid discretization matrix, then the corrected iterate
~
possesses the error propagation operator C := I \Gamma
I h
h A h . Denoting the eigenvalues of A h and A H by f- h
respectively, we may summarize the action of C on the eigenvectors using (2.9) and
as follows:
Theorem 2.1. The image of the fine-grid eigenfunctions fv h
under the
error propagation operator C of the exact coarse grid correction is given by
As a consequence, the two-dimensional spaces spanned by a smooth mode and its
complementary mode are invariant under
The following result shows the dependence of the matrices C j on
Theorem 2.2. Using the notation defined above, there holds
Moreover,
lim
Proof. Both (2.13) and (2.14) are simple consequences of (2.12) and the representation
(2.2) of the eigenvalues - h
.
Application of the error propagation operator to a smooth mode v h
If the entries of the first column of C j are small, then this mode is damped by
the coarse grid correction. However, if the (1; 1)-entry is large then this mode is
amplified, and if the (2; 1)-entry is large (somewhat less likely), then the smooth
mode is corrupted by its complementary mode. As seen from (2.13), these difficulties
occur whenever - H
is small in magnitude. From the limits (2.14), it is evident that no
such problems arise in the symmetric positive-definite case (a fact that is well-known),
but they also do not occur when kh is very large, i.e., when the coarse grid Helmholtz
operator is negative definite. These observations can be extended by returning to
(2.8) and using (2.2), wherein it holds that
That is, the coarse-grid correction strongly damps smooth error modes for either very
small or very large values of kh, but it may fail to do so in the intermediate range
associated with a smooth mode.
We also note that in the limit the eigenvalues of C j are 0 and 1, so that
C j is a projection, and in this case the projection is orthogonal with respect to the
inner product induced by the symmetric and positive definite operator A h . The
projection property is lost for k ? 0, since the coarse grid operator as we have defined
it fails to satisfy the Galerkin condition A
h A h I h
H . (The Galerkin condition is,
however, satisfied e.g. for finite element discretizations with interpolation by inclusion)
Moreover, regardless of the type of discretization, the term A h -orthogonality ceases
to makes sense once k is sufficiently large that A h is indefinite.
2.2.2. Number of Sign Changes. In this section, we discuss the number of
eigenvalues that undergo a sign change during the coarsening process, and thereby
inhibit the effectiveness of coarse grid correction. This is the only aspect of the
algorithm that significantly depends on the dimensionality of the problem. Thus, here
we are considering the Helmholtz equation (1.1) on the d-dimensional unit cube (0; 1) d ,
with homogeneous Dirichlet boundary conditions. We consider standard
finite differences (second order three-point, five-point or seven-point discretization of
the Laplacian in one, two or three dimensions, respectively), as well as the class of
low order finite elements consisting of linear, bilinear or trilinear elements.
We first state the issue more precisely using finite differences. In d dimensions,
the eigenvalues of the discrete operator on a grid with mesh size h and N grid points
in each direction are
d
sin
For any fixed multi-index I, this eigenvalue is a well-defined function of h that converges
to the corresponding eigenvalue of the differential operator as h ! 0. Our
concern is the indices for which this function changes sign, for these are the troublesome
eigenvalues that are not treated correctly by some coarse grid correction. As
the mesh is coarsened, the oscillatory modes (j i ? N=2 for some i) are not represented
on the next coarser mesh, but the smooth-mode eigenvalues f- H
I g are slightly
shifted to the left with respect to their fine-grid counterparts f- h
I g, and some of these
eigenvalues change sign at some point during the coarsening process.
The following theorem gives a bound, as a function of k, on the maximal number
eigenvalue sign changes occurring on all grids.
Theorem 2.3. For finite difference discretization of the Helmholtz equation with
Dirichlet boundary conditions on the unit cube in d dimensions 3), the
number of eigenvalues that undergo a change in sign during the multigrid coarsening
process is bounded above by
3:
For the finite element discretizations, the number of sign changes is bounded above by
pj
3:
Proof. For finite differences, let
fine denote the number of negative eigenvalues
on some given fine grid, and let
lim denote the number of negative eigenvalues of
the continuous Helmholtz operator. Because eigenvalues (2.16) with the same index
I shift from right to left with grid coarsening, it follows that
this is an equality for all fine enough grids, as the discrete eigenvalues tend to the
continuous ones. To identify
lim , consider the continuous eigenvalues
It is convenient to view the indices of these eigenvalues as lying in the positive orthant
of a d-dimensional coordinate system. The negative eigenvalues are contained in the
intersection of this orthant with a d-dimensional sphere of radius k=- centered at
the origin. Let N denote this intersection, and let -
N denote the d-dimensional cube
enclosing N . The number of indices in -
N is bk=-c d , and the number in N is aebk=-c d ,
where
is the ratio of the volume of N to that of -
N . It follows that
3:
Now consider the eigenvalues of discrete problems. Again, since sign changes
occur from right to left with coarsening, the mesh size that yields the maximum
number of negative eigenvalues is the smallest value h for which the discrete operator
is negative semidefinite. With N mesh points in each coordinate direction, this is
equivalent to
d sin 2 N-h
3:
Thus,
d=k, and
d
3:
Combining (2.19) with the fact that
fine
fine
The latter difference, shown in (2.17), is then a bound on number of sign changes.
For finite elements, we are concerned with the eigenvalues of the coefficient matrix
A h , but it is also convenient to consider the associated operator A h defined on the
finite element space V h . The eigenvalues of A h are those of the generalized matrix
eigenvalue problem
A h u
where M h is the mass matrix. These eigenvalues tend to those of the continuous
operator. Moreover, since V H is a subspace of V h , the Courant-Fischer min-max
theorem implies that eigenvalues oe h and oe H with a common index shift to the right
with coarsening (or to the left with refinement). In addition, since M h is symmetric
positive-definite, Sylvester's inertia theorem implies that the number of negative
eigenvalues of A h is is the same as that of (2.21). It follows from these observations
Sign changes
Sign changes
Fig. 2.3. Indices of eigenvalues undergoing a sign change during coarsening of an N \Theta N finite
element grid with during further coarsening of the next coarser
(n \Theta n) grid with
that the maximal number of negative eigenvalues of A h is bounded above by the fine
grid limit
lim .
This is also a bound on the number of sign changes. It can be improved by
examining the eigenvalues of A h more closely. Using the tensor product form of the
operators, we can express these eigenvalues as
1), the indices run from 1 to N and
Consider the requirement - h
so that A h is negative
semidefinite. This is equivalent to
Since the expression - j =- j is monotonically increasing with j, the largest eigenvalue
in d dimensions equals zero if
-N
12d=k. For this value of h, there are
12d) d negative eigen-
values, and on coarser meshes, the problem remains negative definite. Consequently,
none of these
quantities undergo a sign change, giving the bound j \Gamma
of
(2.18).

Figure

2.3 gives an idea of how sign changes are distributed for bilinear elements in
two dimensions. At levels where the changes take place, the indices of the eigenvalues
lie in a curved strip in the two-dimensional plane of indices. Typically, there is one
level where the majority of sign changes occur. As k is increased and h decreased
correspondingly via (2.5), the shape of these strips remains fixed but the number of
indices contained in them grows like O(h \Gammad however, that (2.5) is
not needed for the analysis.) The behavior for finite differences is similar.
The remedy suggested in [11] for these difficulties consists of maintaining an
approximation of the eigenspace V H of the troublesome eigenvalues. A projection
scheme is then used to orthogonalize the coarse grid correction against V H , and the
coefficients of the solution for this problematic space are obtained separately. Since
it involves an explicit separate treatment of the problematic modes, this approach is
restricted to cases where there are only very a small number of these.
3. Incorporation of Krylov Subspace Methods. In view of the observations
about smoothing in Section 2.1 and coarse grid correction in Section 2.2, we modify
the standard multigrid method in the following way to treat Helmholtz problems:
ffl To obtain smoothers that are stable and still provide a strong reduction of
oscillatory components, we use Krylov subspace iteration such as GMRES as
smoothers on intermediate grids.
ffl To handle modes with eigenvalues that are either close to the origin on all
grids-and hence belong to modes not sufficiently damped on any grid-or
that cross the imaginary axis and are thus treated incorrectly by some coarse
grid corrections, we add an outer iteration; that is, we use multigrid as a
preconditioner for a GMRES iteration for (1.4).
We will demonstrate the effectiveness of this approach with a series of numerical
experiments. In all tests the outer iteration is run until the stopping criterion
is satisfied, where Aum is the residual of the mth GMRES iterate and the
norm is the vector Euclidean norm. The multigrid algorithm is a V-cycle in all cases;
the smoothing schedules are specified below.
3.1. GMRES Accelerated Multigrid. We begin with an experiment for the
one-dimensional Helmholtz equation on the unit interval with forcing term
inhomogeneous Dirichlet boundary condition on the left and Sommerfeld
condition on the right. We discretize using linear finite elements on a uniform grid,
where the discrete right hand side f is determined by the boundary conditions. We
apply both a V-cycle multigrid algorithm and a GMRES iteration preconditioned by
the same V-cycle multigrid method. The smoother in these tests is one step of damped
Jacobi iteration for both presmoothing and postsmoothing, using in
(2.6). The initial guess was a vector with normally distributed entries of mean zero
and variance one, generated by the Matlab function randn.

Table

3.1 shows the iteration counts for increasing numbers of levels beginning
with fine grids containing elements and for wave numbers
which correspond to two and four wavelengths in the unit inter-
val, respectively. We observe first that both methods display typical h\Gammaindependent
multigrid behavior until the mesh size on the coarsest grid reaches kh -=2. (With
256 elements, this occurs for coarsest mesh 1=8, and for
coarsest 1=16). At this point both methods require noticeably
more iterations, the increase being much more pronounced in the stand-alone
multigrid case. When yet coarser levels are added, multigrid diverges, whereas the
# levels MG GMRES MG GMRES MG GMRES MG GMRES

Table
Iteration counts for multigrid V-cycle as a stand-alone iteration and as a preconditioner for
GMRES applied to the one-dimensional model Helmholtz problem, with damped Jacobi smoothing.
A dash denotes divergence of the iteration.
128 \Theta 128 elements 256 \Theta 256 elements
# levels MG GMRES MG GMRES MG GMRES MG GMRES

Table
Iteration counts for the two-dimensional problem for fine grids with
128 \Theta 128 and 256 \Theta 256 meshes. A dash denotes divergence of the iteration.
multigrid preconditioned GMRES method again settles down to an h-independent
iteration count, which does, however, increase with k.

Table

3.2 shows the same iteration counts for the two-dimensional Helmholtz
problem on the unit square with a second order absorbing boundary condition (see
[1, 13]) imposed on all four sides and discretized using bilinear quadrilateral finite
elements on a uniform mesh. Since the problem cannot be forced with a radiation
condition on the entire boundary, in this and the remaining examples of Section 3,
an inhomogeneity was imposed by choosing a discrete right hand side consisting of
a random vector with mean zero and variance one, generated by randn. The initial
guess was identically zero. (Trends for problems with smooth right hand sides were
the same.) In addition, for all two-dimensional problems, we use two Jacobi pre- and
postsmoothing steps whenever Jacobi smoothing is used. The damping parameter !
is chosen to maximize damping of the oscillatory modes. For the grids on which we
use damped Jacobi smoothing this optimum value was determined to be 8=9. The
results show the same qualitative behavior as for the one-dimensional problem in that
stand-alone multigrid begins to diverge as coarse levels are added while the GMRES-
accelerated iteration converges in an h-independent number of iterations growing with
# elements on coarsest grid 512 256 128 64
GMRES iterations 152 78 42 25

Table
As more coarse grid information is used, the number of iterations decreases, for the one-dimensional
problem with and a fine grid containing
k, although with a larger number of iterations than in the one-dimensional case.
A natural question is whether corrections computed on the very coarse grids, in
particular those associated with mesh widths larger than 1/10 times the wavelength
any contribution at all towards reducing the error. We investigate this
by repeating the GMRES accelerated multigrid calculations for the one-dimensional
problem with omitting all calculations-be they smoothing
or direct solves-on an increasing number of coarse levels. The results are shown
in

Table

3.3. The leftmost entry of the table shows the iteration counts when no
coarse grid information is used, i.e., for GMRES with preconditioning by two steps of
iteration. Reading from left to right, subsequent entries show the counts when
smoothings on a succession of coarser grids are included, but no computations are
done at grid levels below that of the coarsest grid. For the rightmost entry, a direct
solve was done on the coarsest mesh; this is a full V-cycle computation. The results
indicate that the computations on all grids down to that at level 2, which has eight
elements and only two points per wavelength, still accelerate the convergence of the
outer iteration.
These results show that, although multigrid by itself may diverge, it is nevertheless
a powerful enough preconditioner for GMRES to converge in an h-independent number
of steps. Two additional questions are whether replacing the unstable Jacobi smoother
with a Krylov subspace iteration leads to a convergent stand-alone multigrid method,
and how sensitive convergence behavior is as a function of the wave number k. We
address the former in the following section.
3.2. GMRES as a Smoother. In this section we replace the unstable Jacobi
smoother with GMRES smoothing. We use GMRES on all levels j where kh j -
1=2 and continue using damped Jacobi relaxation when kh choice
is motivated by the discussion at the end of Section 2.1.1, and it ensures that the
largest amplification factor for the Jacobi smoother does not become too large. The
results of Section 2.1.2 show that we could switch back to Jacobi smoothing for very
coarse meshes, but we have not explored this option.
3.2.1. Nonconstant Preconditioners. This introduces a slight complication
with regard to the outer GMRES iteration when multigrid is used as a preconditioner.
The inner GMRES smoothing steps are not linear iterations, and therefore a different
preconditioner is being applied at every step of the outer iteration. A variant of
GMRES able to accommodate a changing preconditioner (known as flexible GMRES
is due to Saad [24]. It requires the following minor modification of the
standard (right preconditioned) GMRES algorithm: if the orthonormal basis of the
(m+1)st Krylov space Km+1 (AM in the case of a constant preconditioner M is
denoted by then the Arnoldi relation AM
Hm
holds with an (m upper Hessenberg matrix ~
Hm . If the preconditioning and
matrix multiplication step
z
is performed with a changing preconditioner results in the modified
Arnoldi relation
The residual vector is now minimized over the space
need no longer be a Krylov space. This requires storing
the vectors fz j g in addition to the orthonormal vectors fv j g, which form a basis of
3.2.2. Hand-Tuned Smoothing Schedules. Numerical experiments with a
fixed number of GMRES smoothing steps at every level did not result in good perfor-
mance. To get an idea of an appropriate smoothing schedule, we proceed as follows.
For given k, we calculate the number o max of FGMRES iterations needed with j-
level multigrid preconditioning, where we use Jacobi smoothing on all grids for which
do a direct solve at the next coarser grid, making j grids in all. We
then replace the direct solve on the coarsest grid of the j-level scheme with GMRES
smoothing on this grid, coupled with a direct solve on the next coarser grid, and
determine the smallest number m j of GMRES smoothing steps required for the outer
iteration to converge in omax steps. For example, for the first line of Table 3.4, 6
outer FGMRES steps were needed for a 5-level scheme, and then m
smoothing steps were needed for the outer iteration of the new 6-level preconditioner
to converge in 6 steps. When the number m j has been determined, we could fix
the number of GMRES smoothing steps to m j on this grid, add one coarser level,
determine the optimal number of GMRES smoothing steps on the coarser grid and
continue in this fashion until the maximal number of levels is reached. This approach
is modified slightly by, whenever possible, trying to reduce the number of smoothings
on finer levels once coarser levels have been added. This is often possible, since replacing
the exact solve on the coarsest grid with several GMRES smoothing steps often
has a regularizing effect, avoiding some damage possibly done by an exact coarse
grid correction in modes whose eigenvalues are not well represented on the coarse
grid. This hand-tuning procedure gives insight into the best possible behavior of this
algorithm.
In contrast to classical linear smoothers, whose damping properties for different
modes is fixed, the damping properties of GMRES depend on the initial residual. In
particular, since GMRES is constructed to minimize the residual, it will most damp
those modes that lead to the largest residual norm reduction. For this reason, we
will favor post-smoothing over pre-smoothing to prevent the unnecessary damping of
smoother modes that should be handled by the coarse-grid correction. We do include
two GMRES pre-smoothing steps to avoid overly large oscillatory components in the
residual prior to restricting it to the next lower level, which could otherwise lead to
spurious oscillatory error components being introduced by the coarse grid correction.
The results are shown in Table 3.4. The entry 'D' denotes a direct solve on the
corresponding level and 'J' indicates that damped Jacobi smoothing was used on this
level. Looking at the smoothing schedules, we observe a 'hump' in the number of
GMRES smoothing steps on the first two levels on which GMRES smoothing is used.
Below this, the number decreases and is often zero for the coarsest levels. However,
256 \Theta 256,
# levels Smoothing schedule MG FGMRES
9 J J J J 13
128 \Theta 128,
# levels Smoothing schedule MG FGMRES
256 \Theta 256,
# levels Smoothing schedule MG FGMRES
256 \Theta 256,
# levels Smoothing schedule MG FGMRES
6 J J
9 J J
256 \Theta 256,
# levels Smoothing schedule MG FGMRES
9 J

Table
Manually optimized GMRES smoothing schedule for the two-dimensional model Helmholtz prob-
lem: 'J' denotes Jacobi smoothing and 'D' denotes a direct solve. The FGMRES algorithm uses the
multigrid V-cycle as a preconditioner.
GMRES smoothing still helps on levels which are extremely coarse with regard to
resolution of the waves: in the case performing three GMRES
smoothing steps on level 4 (which corresponds to 1/2 point per wavelength) still
improves convergence.
We remark that the number of outer iterations in all these tests, for both preconditioned
FGMRES and standalone MG, is the same as for the corresponding two-grid
versions of these methods, so we cannot expect faster convergence with respect to the
wave number k. We also note that the number of iterations for standalone multigrid
is very close to that that for FGMRES with multigrid preconditioning. We believe
this is because the relatively large number of GMRES smoothing steps on intermediate
levels eliminates lower frequency errors, and this mitigates the effects of axis
crossings. We will return to this point in Section 3.4
3.3. A Stopping Criterion Based on L 2 -Sections. Hand tuning as in the
previous section is clearly not suitable for a practical algorithm. In this section, we
develop a heuristic for finite element discretizations that automatically determines
a stopping criterion for the GMRES smoother. This technique is based on an idea
introduced in [29].
We briefly introduce some standard terminology for multilevel methods applied
to second order elliptic boundary value problems on a bounded
domain\Omega ae R 2 (see
[30]). We assume a nested hierarchy of finite element spaces
in which the largest space V J corresponds to the grid on which the solution is sought.
We require the L 2 -orthogonal projections defined by
where (\Delta; \Delta) denotes the L 2 -inner product on \Omega\Gamma Let \Phi
denote the
basis of the finite element space V ' of dimension n ' used in defining the stiffness and
mass matrices. By the nestedness property V ' ae V '+1 , there exists an n '+1 \Theta n '
matrix I '+1
whose columns contain the coefficients of the basis \Phi ' in terms of the
basis \Phi '+1 , so that, writing the bases as row vectors,
The stopping criterion we shall use for the GMRES smoothing iterations is based
on the representation of the residual r ' of an approximate solution ~ u ' of the level-'
equation as the sum of differences of L 2 -projections,
which we refer to as residual sections. The following result for coercive problems,
which was proven in [29], shows that the error u
u ' is small if each appropriately
weighted residual section is small:
Theorem 3.1. Assume the underlying elliptic boundary value problem is H 1 -
elliptic and H 1+ff -regular with ff ? 0. Then there exists a constant c independent of
the level ' such that the H
-norm of the error on level ' is bounded by
The boundary value problem (1.1)-(1.3) under consideration is not H 1 -elliptic
and therefore does not satisfy the assumptions of this theorem. We have found,
however, that the bound (3.1) suggests a useful stopping criterion: terminate the
GMRES smoothing iteration on level ' as soon as the residual section (Q
has become sufficiently small. To obtain a formula for the computation of these
sections, assume the residual r ' is represented by the coefficient vector r ' in terms of
the dual basis of \Phi ' . The representation of Q with respect to the dual basis of
\Phi '\Gamma1 is then given by the coefficient vector I
Returning to the representation with respect to the basis \Phi '\Gamma1 requires multiplication
, so that we obtain
I
If the sequence of triangulations underlying the finite element spaces V ' is quasi-
uniform, then the mass matrix of level ' is uniformly equivalent to the identity scaled
by h d , where d denotes the dimension of the domain. For the case
consideration, this means that the Euclidean inner product on the coordinate space
denoted by (\Delta; \Delta) E , when scaled by h 2
' , is uniformly equivalent (with respect to
the mesh size) to the L 2 -inner product on V ' . Therefore, the associated norms satisfy
ch 2
where v ' is the coordinate vector of v ' with respect to \Phi ' . Using this norm equivalence
it is easily shown that
I
I
for some constants c and C uniformly for all levels '. As a result, the residual sections
may be computed sufficiently accurately without the need for inverting mass matrices.
In [29], it was suggested that the GMRES smoothing iteration for a full multigrid
cycle be terminated as soon as the residual section on the given level is on the order
of the discretization error on that level. For the problem under consideration here, we
shall use the relative reduction of L 2 -sections as a stopping criterion, so that roughly
an equal error reduction for all modes is achieved in one V-cycle. On the first level
on which GMRES smoothing is used, we have the additional difficulty that many
eigenvalues may be badly approximated on the next-coarser level. For this reason,
it is better to also smooth the oscillatory modes belonging to the next lower level
and base the stopping criterion on the residual section
use this 'safer' choice on all levels. Numerical experiments with optimal smoothing
schedules have shown the relative reduction of this residual section to scale like kh ' ,
so that we arrive at the stopping criterion
I '
I
I
A complete description of the multigrid V-cycle algorithm starting on the finest
level ' is as follows:
Algorithm 3.1. ~
V-cycle with GMRES smoothing
on coarse levels
~
else
steps of damped Jacobi smoothing to obtain u (1)
else
perform 2 steps of GMRES smoothing to obtain u (1)
endif
steps of damped Jacobi smoothing to obtain ~
else
perform GMRES smoothing until stopping criterion (3.2) is satisfied
or to obtain ~
endif
endif
In the standalone multigrid V-cycle, Algorithm 3.1 is used recursively beginning
with the finest level and iterated until the desired reduction of the relative residual is
achieved on the finest level. In the FGMRES variant, Algorithm 3.1 represents the
action of the inverse of a preconditioning operator being applied to the vector f ' .
3.4. Experiments with Automated Stopping Criterion. We now show how
the multigrid solver and preconditioner perform with the automated stopping criterion
for GMRES smoothing. Each method is applied to the two-dimensional Helmholtz
problem on the unit square with second-order absorbing boundary condition and
random right hand side data. In these tests, we used in (3.2), and we
also imposed an upper bound mmax on the number of GMRES smoothing steps,
terminating the smoothing if the stopping criterion is not satisfied after mmax steps;
we tested two values, levels, where damped
Jacobi smoothing is used, the number of pre-smoothings and post-smoothings was
2.
Standalone MG MG-preconditioned FGMRES
Nnk 2- 4- 8- 16- 32- 64- 2- 4- 8- 16- 32- 64-

Table
Iteration counts for standalone multigrid and multigrid-preconditioned FGMRES for various
fine grid sizes and wave numbers. In all cases, GMRES smoothing is performed on levels for which
kh ? 1=2 and the smoothing is terminated by the L 2 -section stopping criterion or when mmax
smoothing steps are reached.
This count was extrapolated from the maximum of 47 steps that memory constraints permitted.
We present three sets of results. Table 3.5 shows iteration counts for a variety
of wave numbers and mesh sizes. Table 3.6 examines performance in more detail
by showing the automatically generated smoothing schedules for two wave numbers,
Finally, to give an idea of efficiency, Table 3.7 shows an estimate
for the operation counts required for the problems treated in Table 3.6.
Grid # levels Smoothing schedule Iterations
128 \Theta 128 7 J J
128 \Theta 128 7 J J 17 17 11 2 D 13
Grid # levels Smoothing schedule Iterations
128 \Theta 128 7 J J 19 17 11 2 D 9
128 \Theta 128 7 J J 17 17 11 2 D 13
Grid # levels Smoothing schedule Iterations
512 \Theta 512 9 J J 33 37
512 \Theta 512 9 J J 31 37
Grid # levels Smoothing schedule Iterations
512 \Theta 512 9 J J 20 20
512 \Theta 512 9 J J 20 20

Table
Smoothing schedules with automated stopping criterion, for selected parameters.
We make the following observations on these results:
ffl For low wave numbers, the number of iterations of standalone multigrid is
close to that for FGMRES. The difference increases as the wave number in-
creases, especially for the case 20. For large enough k, multigrid fails
to converge whereas MG-preconditioned FGMRES is robust. This behavior
is explained by the results of Section 2.2.2. For large wave numbers, the increased
number of amplified modes eventually causes standalone multigrid to
fail; a larger number of smoothing steps mitigates this difficulty, presumably
by eliminating some smooth errors. The (outer) FGMRES iteration handles
this situation in a robust manner.
ffl The automated stopping criterion leads to smoothing schedules close to those
obtained by hand tuning (see Table 3.4), and correspondingly similar outer
iteration counts.
ffl The operation counts shown in Table 3.7 suggest that MG-preconditioned
FGMRES is more efficient than standalone multigrid even when the latter
Grid MG FGMRES MG FGMRES
64 \Theta 64 13.2 13.3
128 \Theta 128 24.0 22.1
256 \Theta 256 61.2 43.2 1091.2 971.1
512 \Theta 512 196.6 148.1 1418.1 1377.8

Table
Operation counts (in millions) for selected parameters, with
method is effective.
ffl For fixed wave number, outer iteration counts are mesh independent, so that
standard "multigrid-like" behavior is observed. Moreover, because Jacobi
smoothing is less expensive than GMRES smoothing, during the initial stages
of mesh refinement the costs per unknown are increasing at less than a linear
rate.
ffl The growth in outer iteration counts with increasing wave number is slower
than linear in k. The operation counts increase more rapidly, however, because
of the increased number of smoothing steps required for larger wave
numbers.
4. Application to an Exterior Problem. As a final example we apply the
algorithm to an exterior scattering problem for the Helmholtz equation as given in
(1.1)-(1.3). The
domain\Omega consists of the exterior of an ellipse bounded externally by
a circular artificial boundary \Gamma 1 on which we impose the exact nonlocal Dirichlet-to-
Neumann (DtN) boundary condition (see [19]). The source function is forcing
is due to the boundary condition on the boundary \Gamma of the scatterer, given by
with data g(x; representing a plane wave incident at angle ff
to the positive x-axis. The solution u represents the scattered field associated with the
obstacle and incident field g; the resulting total field u+g then satisfies a homogeneous
Dirichlet or Neumann boundary condition on \Gamma, respectively. An angle of incidence
was chosen to avoid a symmetric solution. The problems were discretized
using linear finite elements beginning with a very coarse mesh which is successively
refined uniformly to obtain a hierarchy of nested finite element spaces. The finest
mesh, obtained after five refinement steps, contains 32768 degrees of freedom. Several
combinations of k and h were tested, where in each case kh ! 0:5 on the finest
mesh.

Figure

4.1 shows a contour plot of the solution u of the Dirichlet problem for
8-. The computations make use of the PDE Toolbox of the Matlab 5.3 computing
environment.
The problems were solved using both the standalone and FGMRES-accelerated
versions of multigrid, with GMRES smoothing using the residual section stopping
criterion with stopping criterion requiring residual reduction by a
\Gamma6 as in Section 3, and zero initial guess. In all examples, we used the
maximal number of levels with the exception of the Dirichlet problem for
where we also varied the number of levels from six down to two. The results are
shown in Table 4.1. The table gives the wave number k and the length of the ellipse
2-=k. The third column gives the maximum value of kh on
Dirichlet problem
28
4 100 26
Neumann problem
28

Table
Iteration counts for the exterior scattering problem with Dirichlet or Neumann plane wave data
on the boundary of an ellipse for various wave numbers, grid sizes and numbers of levels.
the finest mesh and the fourth column indicates the number of levels used in each
computation. The last two columns list the iteration counts.
We observe that the preconditioned iteration performs well in all cases, with a
growth in number of iterations slower than linear in k. The standalone multigrid
variant performs less well in comparison, requiring more than 100 steps to converge
in several cases and even diverging in one case. This is particularly the case for the
Neumann problem, where the superiority of the preconditioned variant is even more
pronounced. For the Neumann problems we also notice a slight growth in iteration
counts for fixed k and decreasing h.
5. Conclusions. The results of this paper show that the addition of Krylov
subspace iteration to multigrid, both as a smoother and as an outer accelerating
procedure, enables the construction of a robust multigrid algorithm for the Helmholtz
equation. GMRES is an effective smoother for grids of intermediate coarseness, in that
it appears not to amplify any error modes and in addition tends to have a regularizing
effect on the contribution to the coarse grid correction coming from smoothing on a
given level. The combination of our multigrid algorithm as a preconditioner with
FGMRES is effective in handling the deficiencies of standard multigrid methods for
the Helmholtz equation, and the outer FGMRES acceleration is necessary particularly
for high wave numbers. In addition, results in the paper indicate that grids too coarse
to result in a meaningful discretization of the Helmholtz equation may still provide
some useful information for coarse-grid corrections. Using an automated stopping
criterion based on L 2 -sections of the residual leads to smoothing cycles that are close
to hand-tuned optimal smoothing schedules.
An important aspect of our algorithm is that it consists of familiar building blocks
and is thus easily implemented. For very large wave numbers for which the discretiza-
Fig. 4.1. Contour plot of the solution of the Dirichlet problem with wave number
tion must not only keep kh but also k 3 h 2 small, the grid hierarchy will contain more
grids fine enough to use Jacobi smoothing, thus making the algorithm more efficient.
The result is a multigrid method that appears to converge with a rate independent of
the mesh size h and with only moderate dependence on the wave number k. Finally,
the numerical results show that we are able to effectively solve Helmholtz problems
with wave numbers of practical relevance.



--R


A comparison of two multilevel iterative methods for nonsymmetric and indefinite elliptic finite element equations
Sharp estimates for multigrid rates of convergence with general smoothing and acceleration
An iterative method for the Helmholtz equation


The cascadic multigrid method for elliptic problems
On the combination of the multigrid method and conjugate gradients
The analysis of multigrid algorithms for nonsymmetric and indefinite elliptic problems

Multigrid methods for nearly singular and slightly indefinite

Absorbing boundary conditions for the numerical simulation of waves

Multigrid preconditioners applied to the iterative solution of singularly perturbed elliptic boundary value problems and scattering problems

Finite element method for the Helmholtz equation in an exterior domain

Exact non-reflecting boundary conditions
Analysis and comparison of relaxation schemes in robust multigrid and conjugate gradient methods


A Multigrid Preconditioner for Stabilised Discretisations of Advection-Diffusion Problems
Iterative Methods for Sparse linear Systems
GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems
An observation concerning Ritz-Galerkin methods with indefinite bilinear forms
Some estimates of the rate of convergence for the cascadic conjugate-gradient method
Multigrid techniques for highly indefinite equations
On the performance of Krylov subspace iterations as smoothers in multigrid meth- ods
A new class of iterative methods for nonselfadjoint or indefinite problems
On the multilevel spltting of finite element spaces for indefinite elliptic boundary value problems

--TR

--CTR
Jan Mandel , Mirela O. Popa, Iterative solvers for coupled fluid-solid scattering, Applied Numerical Mathematics, v.54 n.2, p.194-207, July 2005
Peng Li , L. T. Pileggi, Efficient harmonic balance simulation using multi-level frequency decomposition, Proceedings of the 2004 IEEE/ACM International conference on Computer-aided design, p.677-682, November 07-11, 2004
multigrid;helmholtz equation;krylov subspace methods