First-Order System Least-Squares for the Helmholtz Equation.
This paper develops a multilevel least-squares approach for the numerical solution of the complex scalar exterior Helmholtz equation. This second-order equation is first recast into an equivalent first-order system by introducing several "field" variables. A combination of scaled L2 and H-1 norms is then applied to the residual of this system to create a least-squares functional. It is shown that, in an appropriate Hilbert space, the homogeneous part of this functional is equivalent to a squared graph norm, that is, a product norm on the space of individual variables. This equivalence to a norm that decouples the variables means that standard finite element discretization techniques and standard multigrid solvers can be applied to obtain optimal performance. However, this equivalence is not uniform in the wavenumber k, which can signal degrading performance of the numerical solution process as k increases. To counter this difficulty, we obtain a result that characterizes the error components causing performance degradation. We do this by defining a finite-dimensional subspace of these components on whose orthogonal complement k-uniform equivalence is proved for this functional and an analogous functional that is based only on L2 norms. This subspace equivalence motivates a nonstandard multigrid method that attempts to achieve optimal convergence uniformly in k. We report on numerical experiments that empirically confirm k-uniform optimal performance of this multigrid solver. We also report on tests of the error in our discretization that seem to confirm optimal accuracy that is free of the so-called pollution effect.
Introduction
. The scalar Helmholtz equation with exterior radiation boundary
conditions describes a variety of wave propagation phenomena. One such phenomenon
is the electromagnetic scattering of time-harmonic waves, which, since the
advent of stealth technology, has helped propel an interest in e#cient numerical solution
techniques. However, such boundary value problems are challenging because
they are both indefinite and non-self-adjoint. Hence, standard numerical procedures
for solving them generally su#er from poor discretization accuracy and slow convergence
of the algebraic solver. Accuracy can be improved by taking meshes that are
refined enough for the given wavenumber, but such refinements are impractical for
many problems. Also, convergence of the iterative solver can be improved by preconditioning
the algebraic system, but finding an appropriate preconditioner is problematic
and, in general, robust iterative solvers for indefinite and non-self-adjoint problems
are di#cult to design.
# Received by the editors May 29, 1998; accepted for publication (in revised form) February 16,
1999; published electronically April 28, 2000. This work was sponsored by the National Science Foundation
under grant DMS-9706866 and the Department of Energy under grant DE-FG03-93ER25165.
http://www.siam.org/journals/sisc/21-5/33977.html
Center for Applied Scientific Computing, Lawrence Livermore National Laboratory, P.O. Box
808, L-661, Livermore, CA 94551 (lee123@llnl.gov).
Applied Math Department, Campus Box 526, University of Colorado at Boulder, Boulder,
CO 80309-0526 (tmanteuf@boulder.colorado.edu, stevem@boulder.colorado.edu, jruge@boulder.
colorado.edu).
1928 LEE, MANTEUFFEL, MCCORMICK, AND RUGE
Standard multigrid methods are no exception. Helmholtz problems tax multi-grid
methods by admitting certain highly oscillatory error components that yield
relatively small residuals. Because these components are oscillatory, standard coarse
grids cannot represent them well, so coarsening cannot eliminate them e#ectively.
Because they yield small residuals, standard relaxation methods cannot e#ectively
reduce them. Compounding these di#culties is the property that the dimension of
the subspace of troublesome components increases with increasing wavenumber k.
One approach to ameliorate these di#culties (cf. [9], [10], [21], and [23]) is to introduce
ray-like basis functions on the coarser grids, via exponential interpolation and
weighting, and to use multiple coarsening, where several coarse grids are used at a
given discretization level to resolve some of these error components. Introduction
of ray basis functions allows the oscillatory error components to be represented on
coarser grids, and successive coarsening with an increasing but controlled number of
grids allows a full range of these components to be resolved on the coarser levels.
Since the magnitude of the wavenumber k dictates the oscillatory nature of the error
components, then judicious introduction of exponential interpolation and weighting
and of new coarse grid problems as a function of k leads to multigrid schemes whose
convergence is k-uniform.
Direct application of such a multigrid algorithm to the scalar Helmholtz equation
would not generally achieve optimal discretization accuracy that is uniform in k. Because
k-uniform optimality of the discretization and multigrid solver is the central aim
of this paper, the multigrid algorithm will be applied instead to a carefully designed
first-order system least-squares (FOSLS) (cf. [11], [12], [21], and [22]) formulation
of the Helmholtz problem. This FOSLS methodology involves recasting the scalar
equation into a first-order system by introducing "field" variables, deriving boundary
conditions for these new variables, and applying a least-squares approach to the
resulting first-order system boundary value problem.
An important consideration in the numerical solution of Helmholtz problems is
the so-called pollution e#ect. The results in [3], [5], and [19] show that, for large
wavenumber k, even when kh # 1, the dispersive nature of the Galerkin finite element
discretization introduces a large phase lead error term and, in fact, a sharp H 1 error
estimate is contaminated by a pollution term unless k 2 h # 1. A stabilized Galerkin
least-squares finite element method was developed in [3], [5], and [19] to reduce the
pollution e#ect and, hence, to ameliorate the e#ects of this restrictive condition.
In the work presented below, because the FOSLS methodology leads to a minimization
principle, a Rayleigh-Ritz principle can be applied. Thus, standard bases
can be used to achieve O(kh) discretization error bounds in the least-squares norm
(that is, the norm determined by the functional itself) for the FOSLS formulation.
Our numerical results will suggest that this pollution-free performance is achieved by
the FOSLS formulation-not only the least-squares norm but in a scaled H 1 norm as
well.
Thus, the aim of this paper is two-fold: k-uniform multigrid convergence factors
and pollution-free discretization accuracy. Since a fast k-uniformly convergent solver
has been the most di#cult to obtain historically, this will be the primary focus of our
work. However, we will also show in our final numerical experiment that the FOSLS
discretization appears to be pollution free.
Substantial literature is available on computational electromagnetics. Algorithms
that apply to the time-dependent Maxwell equations or div-curl systems in general
include those described in [13], [20], [24], [26], and [27].
FOSLS FOR HELMHOLTZ 1929
Studies based on first-order system formulations of the Helmholtz equation include
[16], [18], [21], and [23]. The approach described in [16] does not include a curl
term in the functional, so its discretizations are somewhat restrictive and its standard
multigrid methods could not perform well with large wavenumbers, as expected.
The functional in [18] incorporates a curl expression but not in a way that achieves
uniformity in the discretization (or the multigrid solver, had they analyzed it). In
fact, except for [21] and [23] which is the basis of the work presented in the present
paper, none of the methods cited above were shown to achieve optimal discretization
accuracy and multigrid convergence.
For general literature on FOSLS, see [8], [11], [12], [21], and [22], which treat
least-squares functionals based on L 2 and L 2
(that is, a combination of L 2 and
products. The least-squares functionals described in the next section are
also based on these norms. In fact, the analytical techniques used in these papers are
also used to derive some of the theoretical results presented below. These techniques
enable the establishment of several uniform norm equivalence results that would not
be possible by way of the more restrictive Agmon-Douglis-Nirenberg theory (cf. [2]).
This paper is organized as follows. Section 2 introduces the functional setting, the
first-order system, and the L 2 and L 2
least-squares formulations. In section 3,
a coercivity result for the L 2
defined over an appropriate Hilbert
space is established using a compactness argument. Unfortunately, this coercivity
estimate does not hold uniformly with respect to the wavenumber, so in section 4 we
introduce an appropriate subspace on which we do prove k-uniform coercivity. These
estimates provide an intuitive basis for the nonstandard multigrid method developed
in section 5. In the final section, we report on numerical tests that demonstrate k-uniform
pollution-free optimality of the discretization technique and multigrid solver.
2. FOSLS formulation. Let D # 2 be a bounded domain, which for technical
reasons (e.g., scaled H 1 equivalence of the L 2 FOSLS functional; cf. Theorems 3.2
and 4.2) will be assumed to have a C 1,1 or convex polygonal boundary # i of positive
measure. Consider the exterior Helmholtz boundary value problem
D,
Here, k is the wavenumber, f is in L 2
loc
#r
denotes the derivative in
the radial direction used in the asymptotic Sommerfeld radiation condition (cf. [14]).
To solve (2.1) numerically, the unbounded domain is truncated and the Sommerfeld
radiation condition is approximated. Thus, let B be a ball containing D (D # B)
with C 1,1 boundary # e , so that the region bounded by # i # e is annulus-like.
Assuming now that f # L
2(# , consider the reduced boundary value problem
in# ,
Here, a first-order approximation to the Sommerfeld radiation condition is used.
Assume henceforth that (2.2) has a unique solution in H
Introducing the "field" variable
k #p, (2.2) can be recast into the first-order
system
in# ,
or
where L is the first-order operator corresponding to (2.3). Note that (2.3) is essentially
a scaled form of Maxwell's equations that has been reduced by assuming a transverse
magnetic wave solution. Variable p gives the electric field and variable u gives the
magnetic field.
To analyze (2.3) and the corresponding least-squares formulations, several functional
spaces are required. With s being a nonnegative integer, let H
s(# denote the
usual Sobolev space of order s with norm # s . With denote the
usual Sobolev trace space with norm # t (cf. [1]). Let # denote the L 2 (#) inner
product. The following spaces and scaled H 1 norm will be needed:
0,d
The boundary condition imposed in the definition of
W(# is to be taken in the
trace sense. Finally, let H
0,d denote the completion of L
2(# in the
respective inverse norms
(, w)
#w# 1,k
FOSLS FOR HELMHOLTZ 1931
and
0,d
(, w)
#w# 1,k
It is well known (cf. [17]) that [H(div; # H
1(# is a Hilbert
space in the norm
where
and
Thus, it is also a Hilbert space in the norm
Using the continuity of the trace operators from H
2 (#) and H(div;
(#), it is easy to see that W is also a Hilbert space in the # k norm
(cf. [22]).
Define the L 2 and L 2
least-squares functionals by
f
and
f
#-1,d
respectively. Both functionals incorporate curl terms that, with essentially no increase
in complexity for F and only moderate increase for G, allow for simplified finite
element discretization and multigrid solution processes.
For brevity, the equivalence results of the next two sections focus on G; we will
state but not prove the analogous H 1,k equivalence results for F (cf. [21] for detailed
proofs of similar results).
3. Nonuniform coercivity:
W(#3 A compactness argument (cf. [7] and [8])
is used here to establish equivalence between the homogeneous part of the L 2
least-squares functional, G(u, p; 0), and the square of the product norm [#] 2
# 1,k .
In what follows, C denotes a generic constant that may change meaning at each
occurrence but is independent of k. The notation C k will be used (in this section
when dependence on k may occur.
Theorem 3.1. There exist constants C and C k such
for all (u, p) # W
1932 LEE, MANTEUFFEL, MCCORMICK, AND RUGE
Proof. Let (u, p) #
W(# . To prove the upper bound, note that the triangle
inequality yields
-1,d
# .
Integration by parts, the triangle and Cauchy-Schwarz inequalities, the trace theorem,
and the #-inequality give
0,d
#w# 1,k
0,d
u, #w # 1
#w# 1,k
0,d
#w# 1,k
0,d
#w# 1,k
0,d
Also, integration by parts and the boundary condition on H 1
give
#w# 1,k
#w# 1,k
#w# 1,k
#u# .
Finally,
#p# -1,d #p#
#p# 1,k .
The upper bound then follows from (3.1)-(3.4).
To prove the lower bound, let
W(# be the completion of
W(# in the #u#
norm. To prove first that
FOSLS FOR HELMHOLTZ 1933
on
W(#5 note that the triangle inequality used twice and (3.4) yield
# .
Hence, (3.5) would follow if it could be shown that
To this end, note that integration by parts producesk 2 #p#
Moreover, since 1
real, then the right-hand side of (3.8) must be real, which,
together with the triangle inequality, complex modulus arithmetic, and the Cauchy-Schwarz
and #-inequalities, yieldsk 2 #p#
so bound (3.5) holds on all of
W(#1 Now, to prove that
1934 LEE, MANTEUFFEL, MCCORMICK, AND RUGE
on
W(# by contradiction, suppose that (3.9) does not hold. Then there exists a
sequence
W(# such that
and
then the Rellich selection theorem (cf. [7])
implies that p j # p for some p in L
2(#3 which, with (3.5) and (3.10), shows that
is a Cauchy sequence in the #u#p# 1,k
has a limit {(u, p)} #
. Also, consider the scaled Helmholtz equation
in# ,
with the corresponding bilinear form
0,d .
Integration by parts gives
Convergence of p j to p in the # 1,k norm, the Cauchy-Schwarz and #-inequalities,
and (3.10) then imply that
lim
FOSLS FOR HELMHOLTZ 1935
or, equivalently, because of the solution uniqueness of (3.11). Finally, by (3.5),
it follows that
so that
which is a contradiction. Thus, (3.9) must hold. The theorem now follows by restricting
this bound to the subspace
The analogous result for F is essentially obtained using a Helmholtz decomposition
of u. Its proof is very similar to the proofs of Thorems 3.5 and 4.5 of [21], and,
thus, will not be included here.
Theorem 3.2. There exist constants C and C k such
for all (u, p) #
W(# .
These theorems show that the homogeneous parts of functionals F and G are
bounded (with constant C) and coercive (with constant 1
Ck ) in
W(# . Unfortunately,
the coercivity constant depends on the wavenumber k, which means that standard
numerical solution processes can degrade in performance as k increases. In the next
section, a subspace of
W(# is introduced on which we prove the coercivity constant
to be independent of k.
4. Uniform coercivity: Z
#2 A fundamental problem with the Helmholtz
equation is that approximate Fourier components with wavenumbers near k produce
relatively small residuals; that is, oscillatory components are in the near nullspace of
(2.1). This degradation of the usual sense of ellipticity is also present in system (2.3),
which means that the homogeneous parts of respective functionals G and F are not
equivalent to the square of the product norms [#] 2
# 1,k and [# 1,k uniformly
in k. However, as we will see, these troublesome functions can be treated specially in
the multigrid coarsening and finite element approximation processes. To guide the
design of these processes, this section will establish equivalence results on a subspace
of
W(# that excludes these near-nullspace components. To this end, note that
1936 LEE, MANTEUFFEL, MCCORMICK, AND RUGE
is a dense subspace of H 1
0,d (see [6]) that corresponds to the range of the inverse
Laplace operator for
in# ,
(that is,
N(#07 The operator in (4.1) is
self-adjoint, so there exists an L
with lim j#
in# ,
To isolate the near-nullspace components of H 1
0,d(#2 let # (0, 1) be given and define
the spaces
k or #2
Z
Here we use span to mean finite linear combinations and overbar to denote closure in
the H 1,k norm.
The subspace Z
# of
W(#0 which is closed in the #| k norm, avoids the
oscillatory near-nullspace error components of system (2.3) that yield small residuals.
Hence, in this subspace, uniform equivalence statements can be made between the
least-squares norms corresponding to G and F (that is, the square roots of their
homogeneous parts) and the respective product norms [#] 2
Theorem 4.1. Let # (0, 1) be given. Then there exists a constant C, independent
of # and k, such
for all (u, p) # Z
# .
Proof. The upper bound follows from that of Theorem 3.1 and the fact that
Z
W(# .
The proof of the lower bound rests on the observation that Z
# is contained in
the closure of
Z 0,#
in the product norm [#] 2
# 1,k . A proof analogous to that of the upper bound
in Theorem 3.1 shows that G(u, p; 0) is continuous on Z
(# in this norm.
It is therefore su#cient to establish the lower bound for (u, p) # Z
0,#3 To this end,
consider the Helmholtz decomposition (cf. [12])
FOSLS FOR HELMHOLTZ 1937
2(# , and # is divergence-free such that n
# e . Substituting this decomposition into G(u, p; 0) and using the L
orthogonality
property of this decomposition and the divergence-free property of # , we have
Also, integration by parts and the Cauchy-Schwarz inequality give
0,d
#w# 1,k
0,d
#w# 1,k
Now, the triangle inequality and (4.4) yield
which gives
But [-#] is positive-definite and self-adjoint on S # and p # S # , so it is easy to
verify the bound
1,k .
Hence, (4.3), (4.5), and (4.6) imply that
for any (u, p) # Z
0,#1 Finally, combining (4.7) with the triangle inequality and the
restriction # < 1, we obtain
(4.
1938 LEE, MANTEUFFEL, MCCORMICK, AND RUGE
The lower bound now follows from (4.7) and (4.8).
We again state the analogous result for F without proof.
Theorem 4.2. Let # (0, 1) be given and k # 1. Then there exists a constant C,
independent of # and k, such
for all (u, p) # Z #
.
Remark 4.1. While the coercivity constant in Theorem 4.2 is uniform in k,
the dimension of Z # grows with k (for fixed #). However, these troublesome
components can still be specially treated in the multigrid coarsening process without
increasing the order of complexity, as we will show in the next section.
5. Nonstandard multigrid. Some of the basic concepts of this section were
taken from or inspired by the work of Achi Brandt (cf. [9] and Brandt and Livshitz [10]).
The L 2 functional F , when it applies, is usually more practical than the L 2
functional G, which involves the somewhat cumbersome negative norms. For this rea-
son, the remaining two sections will focus on minimizing F (u, p; f). A Rayleigh-Ritz
finite element method is used for the discretization process. Let T h be a triangulation
of
domain# into finite elements of maximal length
let W 1 be a finite-dimensional subspace of
W(# having the approximation property
for (v, q) # H
W(# . The discrete fine grid minimization problem is the
following.
. Find
Equivalently, defining the bilinear form
for (u, p), (v, q) #
W(# , the discrete fine grid variational problem is the following.
. Find
f
for all (v 1 ,
A standard multilevel scheme for solving either of these discrete problems is fairly
straightforward. Let
be a conforming sequence of coarsenings of triangulation T h ,
FOSLS FOR HELMHOLTZ 1939
be a set of nested coarse grid subspaces of W 1 , the finest subspace, and
be a suitable (generally local) basis for W j , referred to here as level j. Given an initial
approximation level j, the level j relaxation sweep consists of the following
cycle:
. For each
chosen to minimize
(By #(b, c) is meant (#), where #
f) is a quadratic function in #, this local minimization procedure
is simple and very inexpensive (see [23] and [25]). In fact, this minimization
process is just a block Gauss-Seidel iteration with blocks determined by the choice of
(b recalling that L is the first-order operator corresponding to (2.3),
the block system is expressed as
where
# .
Now, given a fine grid approximation level 1, the level 2 coarse grid
problem is to find a correction
Having obtained the fine grid approximation is corrected according to
Applying this process recursively yields a multilevel scheme in the usual way.
Unfortunately, if the coarser level basis functions are chosen in the natural way,
then this coarse grid correction process may not be very e#ective for the FOSLS
formulation of the Helmholtz equation. Di#culties can arise with error components
for elements in Z # . Specifically, trouble arises from error components that have the
approximate form
sin ## e k(x cos #+y sin #) ,
where # [0, 2#). The main problem is that, on a grid where
of this type are highly oscillatory and, hence, poorly approximated by standard
1940 LEE, MANTEUFFEL, MCCORMICK, AND RUGE
coarse grids, yet they produce small residuals and are therefore poorly reduced by
relaxation. A compounding problem is that there are, in principle, infinitely many
of these components since # can be virtually any angle in [0, 2#). More precisely,
as k increases, so does the dimension of Z # , which consists of elements that are approximately
of the form in (5.6). Fortunately, these problems can be circumvented
by introducing exponential interpolation and multiple coarsening into the multilevel
scheme ([10] and [23]).
Exponential interpolation. Exponential interpolation is used to approximate
the troublesome oscillatory error components of approximate form (5.6). But the ray
e k(x cos #+y sin #) is the problematic factor of this form, so a more useful characterization
of these error components is given by
s(x, y)e k(x cos #+y sin #) ,
where s(x, y) is a smooth function satisfying appropriate boundary conditions.
j be the components of the product space W j , which is assumed
for concreteness to consist only of continuous piecewise bilinear functions.
(These spaces generally di#er over # only at the boundary.) Also, let level j correspond
to a grid that is fine enough relative to k that components (5.7) can be
adequately approximated with bilinear functions, although only marginally so in the
sense that this approximation just begins to deteriorate on level j + 1. Now, exact
coarse level correction from level j is assumed to be e#ective at eliminating the
smooth components left by relaxation on level j - 1. At level however, the
ability of bilinear functions to approximate components (5.7) deteriorates enough to
contaminate two-level performance between levels j. Some heuristics show
that this occurs when 2 To approximate these exponential components
from levels j basis functions are simply rescaled by the ray
elements sin #) . To be more specific, fix # [0, 2#). For each
l,# denote entry number # of element number # of the standard
basis B l and, on level l, define the level j ray element by
l,# ,
l is the element of W #
that agrees with e k(x cos #+y sin #) d at the
level j node values; that is, E #
(d) is the level j nodal interpolant of e k(x cos #+y sin #) d.
This has the e#ect of altering the usual bilinear interpolation formula that relates
coarse to fine node values so that the rays are better approximated. In fact, this
ensures that the level j interpolant of the given ray is in the range of interpolation on
all coarser levels determined by the corresponding #. Figure 5.1 sketches the real part
of a one-dimensional level 3 basis function d 3 and a typical ray basis element E #
determined by an oscillatory level 1 ray (assuming for simplicity that
Note that exponential interpolation is not much more costly than bilinear inter-
polation: a given function d on level l > j is bilinearly interpolated up to level j, then
simply rescaled by e k(x cos #+y sin #) at the level j nodes.
Multiple coarsening. Consider now the case that level j is where approximation
of components (5.7) just begins to deteriorate. It may be enough to introduce one ray
(e.g., with in the sense that an exact correction from level
might then adequately correct "smooth" level j error. However, it is perhaps more
e#ective here to introduce a few coarse grids, one for each of, say, four rays (e.g., with
FOSLS FOR HELMHOLTZ 1941
ray h
Fig. 5.1. One-dimensional real parts of a typical standard d 3 basis function and a E #
basis element.
critically, as each level is coarsened, scale k2 j-1 h doubles and
essentially twice as many rays become oscillatory. To be more specific, suppose that
are angles corresponding to two neighboring level j rays in the sense that
these two rays were used to coarsen level perhaps others, but none were
needed with angles between # 1 and # 2 . Now, on any of the level
will be one for each ray, that is, for each coarsening), the ray with angle #1 +#2will
be poorly reduced in relaxation and poorly approximated by any of the level j
grids. This means that twice as many rays must somehow be used in coarsening to
level which in turn will again require ray doubling to level
This can be accomplished e#ectively by introducing two separate ray coarsenings
on a level that may itself have come from ray coarsening. Suppose that level
was created from a level j ray with angle #, and that the task is now to coarsen level
using rays of angles #. One way to accomplish this is to introduce the level
ray #
j+2,# . But now #
oscillates on level j scale, so applying
an operator to this function will generally require level j evaluations. Alternatively,
these ray coarsenings can be e#ciently introduced using intermediate level
Suppose the level j rays #
j+1,# have already been introduced on level j + 1. Since
e k[x cos(#)+y
then letting E #
j denote the level j interpolant of the ray perturbation
e k{x[cos(#)-cos #]+y[sin(#)-sin #]}
means that the perturbed rays
can be computed directly from level appealing to level j. In particular,
since
is a level
j+2,# can be constructed by exponentially
interpolating the level Figure 5.2 for a sketch
1942 LEE, MANTEUFFEL, MCCORMICK, AND RUGE
of this process for creating a level 3 ray basis function from levels 1 and 2.) Hence,
given an operator M and assuming that
has been constructed for the ray basis of level j +1, then M#
j+2,# can be constructed
by exponentially weighting the M#
's:
M#
c  d #
c  M#
which is a weighting of elements of (5.9) with exponential weights c  . Recursively
applying this exponential weighting process, analogous coarser level calculations can
be computed without appealing to level j.
having solved a coarse grid problem, its solution can be exponentially interpolated
up to level j using the angle perturbations in a recursive way. For example, if
j+2 is the solution of the level problem corresponding to perturbation
#, then s #
j+2 can be interpolated up to level j by first exponentially interpolating
up to level
then exponentially interpolating this result to level j using #. Since exponential interpolation
consists of a standard bilinear interpolation and an exponential scaling, then
this successive interpolation process is easy to implement. Moreover, since the exponential
interpolation operator is the adjoint of the exponential weighting operator,
then exponential weighting is also easy to implement.
Angles (#). The angles cannot be arbitrary, but must be chosen such that
an optimal number of exponential components are adequately approximated by the
perturbed ray elements. To describe this selection procedure, assume that rays are
first introduced on level j + 1. Then, to obtain adequate approximation, a simple but
tedious analysis [21] shows that the perturbation angle # must be chosen so that the
following approximation is within the level j discretization error:
j+2,# .
For example, for
2 , the level 3 angle
perturbations must be #, as depicted in Figure 5.4. The angle perturbations
are simply halved at each successively finer level.
Nonstandard multigrid scheme and computational cost. Exponential in-
terpolation/weighting and multiple coarsening are essentially all that is required to
FOSLS FOR HELMHOLTZ 1943
d
ray 2h
"perturbed" ray hFig. 5.2. Level 3 ray element determined from a perturbation corresponding to angle # of the
level 2 ray element.
Level (j-1)
Level j
Level (j+1)
Level (j+2)
Fig. 5.3. Multiple coarsening.
develop an e#ective nonstandard multigrid scheme for the Helmholtz equation. With
these components, ray basis elements are naturally introduced, and this nonstandard
multigrid scheme has the same form as the standard multigrid scheme defined by (5.2)
and (5.5). To see this, again assume that rays are first introduced on level
that all perturbation angles are determined before the multigrid cycling begins. Here,
we choose j to be the largest integer satisfying
since this guarantees essentially 8 grid points per wavelength in the piecewise bilinear
approximation on level j. Now, for levels l # j, the nested spaces
still consist of products of piecewise bilinear functions, but, for levels l #
must introduce the product ray element spaces
l := span #r
1944 LEE, MANTEUFFEL, MCCORMICK, AND RUGE
Level 3 angles
Level 2 angles
Fig. 5.4. Angle choices (uniformly distributed) for level 2 and level 3 rays.
where t l is the total number of rays introduced on level l, # r is the angle of the
rth ray, and #r
is the corresponding perturbed vector ray basis
introduced on level l. We first introduce four rays on level double that
number for each successively coarser level, so we have t . Note that
l
only if # s is a perturbation of angle # r (that is, only if As before, for
levels l # j, relaxation is defined by (5.2) with (b l , c l ) # the standard basis for W l and,
for levels l < j, the level (l problem is to find (u l+1 , p l+1 ) # W l+1
such that
(v l+1 ,q l+1 )#W l+1
where l ) is the level l approximation. However, for the ray levels l > j, relaxation
is again defined by (5.2); but now with (b l , c l ) # the product ray basis element #r
l,# ,
and, for levels l # j, the two level (l problems corresponding to each
r , are to find (u l+1 ,
(v l+1 ,q l+1 )#W #s
where
l is now the approximation on level l with angle # r . Note that
each level l with angle # r is corrected by two level (l (actually, four
when This gives the structure of the multigrid process. The scheduling
that determines how the various levels are visited is best described by studying the
schematic in Figure 5.5. We will continue to refer to this process as a V -cycle, although
it looks very di#erent from a standard V-cycle below level j.
For this nonstandard multigrid scheme (which can be interpreted as a multilevel
multiplicative Schwarz method; cf. [28]), the level j coarse grid problem involves
subspace corrections from the standard finite element space W j and the ray subspaces
l
l > j. Also, relaxation on the ray spaces is simply a block Gauss-Seidel iteration
FOSLS FOR HELMHOLTZ 1945
Level (j-1)
Level
Level (j+1)
Level (j+2)
Fig. 5.5. V-cycle for the nonstandard multigrid scheme.
(5.3) with (b l , c l ) #r
l,#
. Since L#r
l,#
can be computed e#ciently using exponential
weighting, so can the elemental block matrices or stencils
l, , L#r
l,# .
In fact, these stencils can be computed and stored prior to any multigrid cycling so
that the same relaxation routine can be used for all levels and angles.
Now, assuming that the stencils have been computed and stored before any multi-grid
cycling occurs, the computational cost of such a multiple coarsening algorithm
is generally not excessive. In fact, one such multilevel V-cycle requires only several
fine grid work units (that is, the cost of a fine grid relaxation). For example, consider
the worst case, when four rays are introduced on level 2. (The need for more rays
would signal a level h so coarse relative to k that the discretization accuracy would
be almost meaningless; in fact, the need for four rays on level 2 is already a sign that
level 1 accuracy is quite poor.) Then, for level j # 2, there would be a total of 2 j
level j problems. Also, for each j, if it is assumed that the stepsize is doubled at each
coarsening, then the number of fine grid work units required for relaxation is 2 -2j+2 .
Omitting the cost of ray interpolation (which is substantial but proportional to the
number of points) and other more minor operations, the total number of fine level
work units for a nonstandard V (1, 0)-cycle is then
which is only about double the cost of a standard V (1, 0)-cycle.
6. Numerical examples. For simplicity and to facilitate comparisons, the examples
treated here are only for the unit square. The lack of a Dirichlet boundary is
in fact a more severe test of the methodology in some sense. The results will include
multigrid performance for both the standard and the nonstandard schemes described
1946 LEE, MANTEUFFEL, MCCORMICK, AND RUGE

Table
Standard coarsening: V (1, 0)-cycles.

Table
Standard coarsening: two-level cycles.
in the previous section. We also assess accuracy of the discretization based on ray-type
elements.
To facilitate assessment of the iterative solver, we first restrict ourselves to the homogeneous
boundary value problem (2.2) with The exact solution is of course
but this is useful for measuring asymptotic algebraic convergence factors of
stationary linear iterative methods such as multigrid methods since it avoids the limiting
stagnation caused by machine representation. The unit square is partitioned
into rectangles, and the functional F is minimized over the space of piecewise continuous
bilinear functions (possibly involving rays) that satisfy the boundary condition
on #. In most cases, the finest grid is chosen to satisfy
8 , and
rays are first introduced on level
8 (that is, typically when
The relaxation scheme is nodal Gauss-Seidel relaxation (that is, block Gauss-Seidel
where a block consists of all of the unknowns (u h
# ) at each node #). To assess the
worst-case convergence factors, twenty multigrid V(1,0)-cycles (that is, one relaxation
before and none after coarsening) were performed starting from a random initial guess
). The convergence measure # is then defined as
where the superscript in parentheses denotes the iteration number. Tables 6.1 and 6.2,
respectively, summarize the results of twenty V (1, 0)- and two-level cycles for standard
bilinear elements. Table 6.1 confirms that convergence degrades as a function of k
due to the troublesome oscillatory error components present on the coarser levels.

Table

6.2 shows that Poisson-like (that is, are obtained when kh is not
too large and, comparing it with Table 6.1, indicates that rays should be introduced
as soon as kh > 1
Tables

6.3 and 6.4 show that multigrid convergence factors improve substantially
by introducing ray functions. Here, "exact" rays (that is,
are intro-
FOSLS FOR HELMHOLTZ 1947

Table
Exact ray (#
0)-cycles.

Table
Exact ray (#
cycles.
.3

Table
Approximate ray (#
0)-cycles.

Table
Approximate ray (#
delayed until kh > 1
0)-cycles.
duced when kh > 1using the algorithm described in the previous section. These
factors are fairly uniform in k and compare favorably to the two-level factors for
standard coarsening. Moreover, they degrade only slightly when the exact rays are
approximated using the recursive scheme described in the previous section (see Table
6.5).
To test the e#ect of introducing the rays late in the multigrid solver, we reran
these examples but delayed the use of rays until kh > 1. As expected, the rates
degrade noticeably, as Tables 6.6 and 6.7 show.
As a measure of discretization accuracy, we ran tests to see if our FOSLS scheme
has any pollution e#ects. We specified the exact solution
5 +y sin #) , which
solves (2.2) with except that the boundary condition is inhomogeneous:
5 +y sin #)
5 +y sin #) on #.
1948 LEE, MANTEUFFEL, MCCORMICK, AND RUGE

Table
Approximate ray (#
delayed until kh > 1
0)-cycles.

Table
Discretization error in the least-squares functional norm.

Table
Relative discretization error in the scaled H 1,k norm.
kh\h 11111
We then used our FOSLS scheme to approximate p and the scaled gradient u =k #p. To isolate possible pollution, we chose various h, then varied k so that kh
is constant 1/64). The errors were measured by comparing the
discrete approximation (obtained after applying several multigrid cycles) to the exact
solution in a relative sense in the least-squares and product H 1,k norms. Tables 6.8
and 6.9 essentially show constant discretization errors for constant kh, with about a
factor of 2 decrease as kh is halved, which is consistent with the assertion that our
approach is pollution-free and O(kh), respectively.



--R

New York
Estimates near the boundary for solutions of elliptic partial di
A generalized finite element method for solving the Helmholtz equation in two dimensions with minimal pollution
The partition of unity method
Is the pollution e
Expansions in Eigenfunctions of Selfadjoint Operators
Finite Elements Theory
A least-squares approach based on a discrete minus one inner product for first order systems
Stages in developing multigrid solutions



Finite element method for the solution of Maxwell's equations in multiple media
Inverse Acoustic and Electromagnetic Scattering Theory
Partial Di
On numerical methods for acoustic problems
Finite Element Methods for Navier-Stokes Equations

Finite element solution of the Helmholtz equation with high wave number.



Multilevel first-order system least-squares (FOSLS) for Helmholtz equations

Multilevel Projection Methods for Partial Di
A mixed method for approximating Maxwell's equations
On electric and magnetic problems for vector fields in anisotropic nonhomogeneous media
Iterative methods by space decomposition and subspace correction
--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
Jan Mandel, An iterative substructuring method for coupled fluid-solid acoustic problems, Journal of Computational Physics, v.177 n.1, p.95-116, March 20, 2002
first-order system least-squares;helmholtz equation;nonstandard multigrid