The Perfectly Matched Layer in Curvilinear Coordinates.
In 1994 Brenger showed how to construct a perfectly matched absorbing layer for the Maxwell system in rectilinear coordinates.  This layer absorbs waves of any wavelength and any frequency without reflection and thus can be used to artificially terminate the domain of scattering calculations.  In this paper we show how to derive and implement the Brenger layer in curvilinear coordinates (in two space dimensions).  We prove that an infinite layer of this type can be used to solve time harmonic scattering problems.  We also show that the truncated Brenger problem has a solution except at a discrete set of exceptional frequencies (which might be empty).  Finally numerical results show that the curvilinear layer can produce accurate solutions in the time and frequency domain.
Introduction
.
Finite element, finite difference or finite volume discretizations of scattering
problems are faced with the problem of truncating the infinite domain problem
without introducing excessive error into the computed solution. Usually
an artificial boundary is introduced on which a special absorbing boundary
condition is used that reduces the reflection of waves incident on the artificial
boundary. The key is to find boundary conditions that do not ruin the stability
of the computation, are not too expensive and are acceptably accurate.
It is difficult to satisfy these three competing criteria. Low order boundary
conditions are usually stable, but in order to obtain sufficient accuracy the
absorbing boundary must be far from the scatterer making the computation
costly. Higher order conditions are more complex to implement and sometimes
have stability problems (for a survey of absorbing boundary conditions
see [16, 14]).
Another approach to terminating the computational domain is to surround
the scatterer by an absorbing layer. In the past such layers have had to be
unacceptably large to control spurious reflections from the surface and interior
of the layer (see for example [20]). However, in 1994, B'erenger [2] showed how
to perturb the Maxwell system to provide a perfectly matched absorbing layer
(abbreviated to PML for the rest of this paper) for electromagnetic scattering
RR n3049
4 Francis Collino et Peter Monk
problems. His approach constructs an absorbing half-space such that waves
incident from the standard Maxwell medium are not reflected from the half
space independent of their frequency and direction of propagation. In the
PML, the wave is absorbed and decays exponentially with distance into the
layer. Thus the PML layer itself can be truncated to form a narrow absorbing
slab with a low reflection for any angle of incidence. The scatterer can then be
surrounded by perfectly matched slabs to truncate the scattering domain. This
approach is easy to implement and is very effective. It has already been the
subject of numerous of papers [18, 1, 5, 12, 11, 13, 25, 26, 23, 21, 15, 7, 4, 24, 3].
However, to our knowledge, there has been very little theoretical analysis of
the method.
Of particular interest from the point of view of this paper is the paper of
Chew and Wheedon [6] in which they show that the B'erenger PML can be
constructed via a complex change of variables (see also [8]). We shall develop
this viewpoint in the paper and show how to extend this idea to include a PML
with a variable absorption coefficient (as is used in practice) and we shall show
that method can be used in curvilinear coordinates (see also [8]). We shall
show how the coordinate change is related to evaluating a continuation of the
solution into the complex plane such that the solution becomes evanescent.
This justifies the change of variables approach.
In this paper we shall show how to view the PML as a complex change of
variables in the time harmonic scattering problem. Using this observation we
shall construct a PML layer in cylindrical coordinates (in fact in two space
dimensions) for the time dependent and time harmonic problems. For some
scatterers the use of circular absorbing boundaries can improve the efficiency
of the scheme by decreasing the area that must be meshed. We shall demonstrate
how the construction of the PML is related to special function expansions
of the solution and also provide a partial theoretical analysis of the time harmonic
problem. Our methods extend to three dimensions and to ellipsoidal
coordinates. These aspects will be reported in future papers.
INRIA
The PML in Curvilinear Coordinates 5
Perfectly matched layers and change of variable

As discussed in the introduction, in [2] B'erenger proposes a new absorbing layer
model called the perfectly matched layer (PML), which has the astonishing
property of generating no reflection at the interface between the free medium
and the artificial absorbing medium.
In this section we shall show how B'erenger's PML can be viewed as a
complex change of variables applied to the Maxwell system. In particular
for the 2D Maxwell's equations, the design of the PML in the region x ? 0
is as follows. The propagation of a transverse electric wave is governed by
@H z
@t
@y
@x
@t
@x
@t
@y
(1)
The construction of an absorbing layer in the right half space (i.e. in the region
performed in two steps. First, we split the magnetic component into
two parts according to
and rewrite equations (1) as
@H zy
@t
@y
@H zx
@t
@x
@t
@x
@t
@y
At this stage, equations (2) and (3) are equivalent to the initial problem.
In the second step, a damping factor oe(x), introduced
in each equation in those places where the normal derivative operator @ x
appears and the other parts of the equations remain unchanged:
@H zy
@t
@y
@H zx
@t
@x
@t
@H z
@x
@t
@H z
@y
RR n3049
6 Francis Collino et Peter Monk
System (4) is B'erenger's PML model.
If we look for the time harmonic solutions with frequency !, we get
@y
@x
H z
@x
H z
@y
H zy
We can rewrite this system of equations with the only unknown being "
H z . We
have
H zy
@y
@x
or
H z
@
@x
H z
@x
In the region where oe is zero, we recover the well known property that the
magnetic field of a non-stationary 2D transverse electric wave satisfies the
Helmholtz equation. Inside the layer, we remark that the PML model consists
in making the simple substitution
@
@x
\Gamma!
@
@
@x
Z xoe()d (7)
in the initial equation, whereas the tangential derivatives are left unchanged.
Equation (6) can be easily solved (even if oe is not constant), we get
z
R x(1+ ioe()
)d dk y ;
y
z
The solutions appears as a superposition of right-going evanescent plane waves.
The waves propagate without any reflection even if the variations of oe are
strong. In particular, at the interface reflection occurs as the wave
INRIA
The PML in Curvilinear Coordinates 7
penetrates the PML. After propagating a distance ffi in the y direction the k y
component of the signal is damped by a factor
y
for the propagating waves (i.e. jk y
y
for the evanescent waves (i.e. jk y j ? !).
In case of a finite PML (it is necessary to truncate the PML to obtain a
problem that can be solved numerically), a Dirichlet boundary condition is
added at the right extremity, say ffi. The solutions can now be sought as
z
x
x
)d
e ikx
R ffi(1+i oe()
R ffi(1+i oe()
)d
z
located above the support of the initial conditions. These solutions
are obtained by combination of plane waves which vanish at are
continuous at The expression of the coefficient of reflection R(k y ) is
obtained by writing down the continuity of the tangential component of the
electric field
which gives
R ffi(1+i ioe()
d
R ffi(1+i oe()
d
RR n3049
8 Francis Collino et Peter Monk
or
R ffi(1+i oe()
We recover B'erenger's result for the expression of the coefficient of reflection
induced by a PML of width ffi. A direct inspection of this formula would
lead us to pick oe to be the largest possible which would allow us to obtain
a reflection as weak as desired. However, this nice property is no longer true
when a discretization scheme for solving numerically the equations is used. As
the matter of fact, see [8], the discretization of the equations makes the layers
slightly imperfectly matched. The numerical dispersion induces a reflection at
the interface which is an increasing function of oe\Deltax. Thus, oe has to be chosen
small enough to control this numerical reflection coefficient and large enough
to obtain a low reflection coefficient (due to the Dirichlet condition on the wall
at
In square domains, B'erenger's PML for 2D transverse electric waves is
constructed with the help of two damping functions oe x (x) and oe y (y) satisfying
oe s
oe s
Equation (4) is replaced by
@H zy
@t
@y
@H zx
@t
@x
@t
@x
@t
@y
Following the same developments as above, we obtain the stationary solution
as the solution of
@
@y
H z
@y
@
@x
H z
@x
0:
A direct inspection shows that the Green's function for this equation is given
by
iH (1)'
~
INRIA
The PML in Curvilinear Coordinates 9
where
Z xoe x (s)ds; ~
Z yoe y (s)ds
This is one key to understanding the nice properties of the PML model. This
model is simply obtained using a complex change in variables applied to the
original equation. The axes x; y 2 R are moved into two paths in the upper
complex plane. The existence of an analytical extension of the Green's function
in this upper complex plane allows us to properly define "
H z . The exponential
decay of this extension ensures the damping effect of this transformation on
the solution.
We end up this section by giving a variational formulation for the harmonic
B'erenger's PML model with Dirichlet boundary conditions. The problem is
posed in D
n\Omega where\Omega is a perfectly conducting
scatterer entirely contained in (] \Gamma The problem is
to find "
H z such that
d x d y
@
@x
dx
Hz
@x
@
@y
dy
Hz
@y
H z
@
\Gammag s 2
where g is a given function and
oe x (x); d y
oe y (y):
If
the variational formulation of the above problem is simply
Z\Omega
d x d y
Z\Omega /
d y
d x
H z
@x
@x
Z\Omega /
d x
d y
H z
@y
@y
Z
@\Omega gF ds
RR n3049
Francis Collino et Peter Monk
This formulation can be discretized by a finite element method in the usual
way way. Concerning the existence and uniqueness of the solution, a result
similar to those we will obtain for B'erenger's PML in curvilinear coordinates
(see Theorems 1 and 2), can be obtained.
In conclusion of this preliminary section, we have shown that the PML
model can be interpreted by a change in variables and that the corresponding
model can be put in a variational form suitable for a finite element method.
With this in mind, we can now begin our study of a PML model in curvilinear
coordinates.
3 The time-harmonic scattering problem
In this paper we will develop a curvilinear PML for a simple electromagnetic
scattering problem using the change of variable technique introduced in the
previous section.
Let\Omega be a smooth bounded domain in the plane, and assume
that the boundary
of\Omega is a perfectly conducting scatter. The problem we shall
consider is to compute the transverse electromagnetic field scattered from the
boundary of \Omega\Gamma We shall start with the frequency domain problem. Let "
denote the magnetic field scattered by
H satisfies the Helmholtz
equation in the exterior of \Omega\Gamma
where the wave number
is the frequency of the incident field, ffl
is the dielectric constant and  the permeability of the surrounding medium.
On the perfectly conducting boundary, Neumann boundary data is given:
@
= g on
where  is the unit outward normal
to\Omega and g is a given smooth function
derived from the incoming or incident field. Finally the scattered field must
satisfy the Sommerfeld radiation condition uniformly in all directions:
@ae
INRIA
The PML in Curvilinear Coordinates 11
is the radial coordinate. For sufficiently smooth g
problem has a unique solution H (R 2 n \Omega\Gamma (see for example
[9]).
Let us assume
that\Omega is contained in the interior of the circle of radius
a 0 . Then for the solution H can be written using separation of
variables:
a
n is the Hankel function of first kind and order zero, and
ae exp(i'). This series converges uniformly for r  a 0 [9].
To obtain the B'erenger PML layer in cylindrical coordinates, we want to
extend the solution given by (20) to the complex plane such that the wave
becomes evanescent as the imaginary part tends to infinity. Suppose ~
~
ae r and ~
ae i are real. The desired asymptotic behavior of the
Hankel functions is
Cm
ik~ae
exp (i!( ~
exp (\Gammak ~
aej !1:
This implies that the correct half-plane is
With this observation in mind, we define the complex radius ~
ae. We choose
a ? a 0 and a real (artificial) absorption coefficient parametrized by ae
satisfying a and oe(ae) ? 0 for ae ? a, then
~
a
ds if ae  a
ae if a
We also assume that
lim
a
Note that
d~ae
dae
if ae  a
RR n3049
12 Francis Collino et Peter Monk
It will prove convenient to define oe by
a
ds if ae  a
Then ~
ae and ae are related by d as follows
~
Now we can define the B'erenger solution "
n\Omega as follows:
We need to derive a boundary value problem satisfied by "
n\Omega
satisfies the standard Helmholtz equation. Across the circle
normal derivative are continuous. Thus it
remains to derive an equation in jxj ? a 0 . In this region "
H B is a sum of Hankel
functions and hence satisfies the Helmholtz equation (17) but using ~
ae in place
of ae
~
ae
@
@~ae
~
ae
@
@~ae
~
ae
We can start by using the chain rule:
@~ae
@ae
@ae
@~ae
=d
@ae
Hence using this result and (21) we can rewrite (23) in terms of ae. We obtain:
dae
/d
@
@ae
dae
d
@
@ae
dae
which we rewrite slightly as
@
@ae
dae
d
@
@ae
d
dae
INRIA
The PML in Curvilinear Coordinates 13
It is convenient (from the point of view of implementation in a standard finite
element code) to transform this equation back to rectangular coordinates to
obtain
where A is the following matrix function:
d
d
sin 2 (') cos(') sin(')
d
d
d
d
sin
d
If ae  a 0 we define A = I (the identity matrix), and The problem of
computing the B'erenger solution "
H B is thus the problem of finding "
that
r
@
= g on
uniformly bounded as jaej !1 (28)
Here the boundedness assumption has replaced the standard Sommerfeld condition
at infinity (because the solution is assumed to be decaying at infinity).
We can now state and prove our first theorem:
Theorem 1 There exists a unique solution of the B'erenger problem (26)-(28)
given by (22).
Remark. This theorem shows that the solution of the curvilinear B'erenger
time-harmonic scattering problem is the desired physical solution of the standard
scattering problem in the neighborhood of the scatterer.
Proof: Existence is not a problem since we have constructed the boundary
value problem (26)-(28) so that (22) (derived from the known solution of the
standard scattering problem) is a solution.
To prove uniqueness we can modify the standard proof as follows. We set
then using the fact that (26) is just
the Helmholtz equation in D a 0
we can use Green's theorem applied to D a 0
to
show that if
Z
@r
@r
ds
RR n3049
14 Francis Collino et Peter Monk
But for r ? a
hence is a sum of Hankel functions:
(a
Since we assume that oe is chosen so that =(k~ae) ! 1 as ae ! 1 the large
argument asymptotics of the Hankel functions (note shows that
H (2)
unbounded as ae ! 1. Hence to satisfy the
radiation condition (28) we must require b using
the Hankel function expansion in (29) shows that
where we have used the fact that ~
ae on a 0  ae  a. Using the non-vanishing
of the Wronskian of the two Hankel functions we conclude that a
n. Hence "
and a standard continuation argument
shows that "
n\Omega (see for example [9]). This completes the proof.The boundary value problem (26)-(28) is still not suitable for computation
since it is posed on an infinite domain. However truncation of this problem is
not as delicate as for the standard scattering problem since the solution decays
exponentially outside the circle ae = a. For example we can impose the Dirichlet
boundary condition on a suitably large radius circle containing the scatter (in
three dimensions this would be equivalent to specifying that the tangential
component of the magnetic field vanishes on the artificial boundary),
Let
On the outer boundary we impose the Dirichlet boundary condition
Then the first boundary value problem we shall approximate numerically is to
H B such that (26) is satisfied in D b with boundary data given by (27) on
@\Omega and by (30) on the artificial boundary
As we shall see, better numerical behavior can be obtained if we use an
absorbing boundary condition on the outer boundary. So to further improve
INRIA
The PML in Curvilinear Coordinates 15
the efficiency of our numerical experiments we use a simple absorbing boundary
condition on the artificial boundary. We impose the Sommerfeld radiation
condition suitable for the evanescent Hankel function:
@ae
Then the field "
the Robin problem (26), (27) and (31).
The obvious questions are whether the problems outlined above have unique
solutions and whether their solution is close to the solution of the infinite
domain B'erenger problem (in particular we want to know if the solution is close
to the true scattering solution close to \Omega\Gamma/ In particular, does the absorbing
boundary condition improve the method. We can provide a partial answer to
the first of these questions as the following theorem shows but we have been
unable to prove the approximation property.
Before stating the theorem we define the constant M ? 0 by
We will prove our theorem by allowing complex valued ! in the following set
ae C which contains the positive real axis with the exception of
physically interesting values of ! lie):
r =M
Recall that
Theorem 2 Suppose
. Then:
1. The boundary value problem (26), (27) and (30) on D b has a unique
solution "
H B in H 1 (D b ) for every real k except possibly for a discrete set
of values of k.
2. The same conclusion holds if "
the Robin problem
(26), (27) and (31).
RR n3049
Francis Collino et Peter Monk
Remark. We have been unable to rule out the possibility that a unique
solution fails to exist for some real k.
Proof: The proof of the two conclusions of the theorem are almost identical,
so we will prove only the case of the Robin problem which is more interesting
in practice.
The proof rests on the analytic Fredholm theorem (see [10]). First we pose
the B'erenger problem (26), (27) and (31) as a weak problem on the space
(for the Dirichlet problem we would use the subspace
Z
Z
@\Omega uv ds and !
Z
uv ds
then we seek "
Now we define two operators
for all u; v 2 H 1 (D b ). Both operators are well defined by the Riez representation
theorem. Note that because d, d and k depend on !, the operators A
and B depend on !.
Since the coefficients defining A are analytic at each ! 2 U , the operator
A is analytic in !. Next we show that A \Gamma1 exists and is continuous by using
the Lax-Milgram lemma. The continuity of A is obvious so it remains to show
that A is strongly coercive. But
and so strong coercivity is proved if there is a positive constant ffi (possibly
depending on !) such that
INRIA
The PML in Curvilinear Coordinates 17
For any u 2 H 1 (D b )
Z
Z
a 0 !jxj!b
But from the definition of the B'erenger layer:
Z
a 0 !jxj!b
a 0
d
d
@ae
d
A d' dr:
It thus suffices to show that
min
where ffi is independent of ae but can depend on !.
We shall prove the above inequality for !(d=d), the other inequality being
similar. Note that if
d
d
Hence !(d=d) is a monotonic function of ! r and so
d
d
min
min
min
r =M then
d
d
r
r
Hence !(d=d)  Having proved the strong coercivity of A, we
can now conclude that A is invertible for any ! 2 U .
The operator B is compact, because if
the following boundary value problem:
@w
@
@\Omega
@w
@r
= \Gammaikdu on
RR n3049
Francis Collino et Peter Monk
Since d and d are functions in L1 (D b ), standard elliptic regularity theory
implies that jjwjj H
1(\Omega\Gamma but since H
2(\Omega\Gamma is compactly embedded
in H 1
(\Omega\Gamma we conclude that B is compact. Furthermore B is an operator valued
analytic function of k for all k.
Now let G 2 H 1
(\Omega\Gamma be defined by
then the B'erenger problem is equivalent to finding "
1(\Omega\Gamma such that
A \Gamma1 is bounded, A \Gamma1 B is compact (and an operator valued analytic
function of k) and so the analytic Fredholm alternative theorem is applicable.
It remains to decide between the alternatives of this theorem. To do this
shall prove that a solution of the B'erenger problem exists when Once
this is proved, the theorem is proved since the analytic Fredholm theorem
implies that exists except possibly for a discrete set of exceptional
wavenumbers.
the functions d and d are real and positive since
Hence the matrix A is real and uniformly positive definite (when
dd ? 0. The weak problem is to find "
But in this case the bilinear form associated with this problem is coercive since
if
Hence the Lax-Milgram lemma shows that a unique solution exists for (34).2
INRIA
The PML in Curvilinear Coordinates 19
4 The Time Dependent Problem.
B'erenger first proposed his PML for the time dependent Maxwell system. Here
we describe how to generalize B'erenger's PML to allow for an annular absorbing
layer in the time domain. The scattering problem is again to predict
transverse mode scattering from an infinite perfectly conducting cylinder occupying
the
region\Omega described in the previous section (although, as we have
seen, the details of the scatterer are unimportant from the point of view of
developing the PML equations).
We assume that at the fields vanish for x an
incident wave interacts with the scatterer to produce a scattered electromagnetic
field where the electric field
and the magnetic field
H M satisfy the two dimensional Maxwell system
@y
in R 2
n\Omega
@x
in R 2
n\Omega
@t
in R 2
n\Omega
with zero initial data and boundary data appropriate for a perfectly conducting
obstacle:
where fl is a given function (which is determined by the incident electromagnetic
field).
As we have seen in the previous section, the curvilinear B'erenger PML
can be derived by writing the frequency domain equations in curvilinear co-ordinates
and changing variables. We start by writing the Maxwell system in
cylindrical coordinates. For a ae the electric field
magnetic field H M satisfy
ae
@t
=ae
@t
@ae
RR n3049
Francis Collino et Peter Monk
@t
@
@ae
ae
To derive the B'erenger medium, we move to the frequency domain. If we
assume that E M
ae exp(\Gammai!t) and similarly for the remaining variables we
obtain the time harmonic Maxwell system:
ae =ae
@ae
@ae
ae
@'A
For the time dependent problem, we proceed analogously to the previous
section. The B'erenger system is just the standard Maxwell system, written
using ~
ae in place of ae (the fields now carry the superscript B to indicate that
we are considering a B'erenger medium):
~
ae
@~ae
ae@ @
@~ae
ae
Note that if we use (35) and (36) to eliminate "
ae and "
' from (37) we obtain
(25). As in the case of the time harmonic problem, we want to derive the
equivalent of this system using real coordinates. Using (24) implies that (36)
can be written:
@ae
We can easily rewrite (35) using the fact that ~
dae to obtain:
ae =ae
INRIA
The PML in Curvilinear Coordinates 21
It remains to rewrite (37). We define the auxiliary B'erenger variable "
H (ae) by
d
@
@ae
Expanding the derivative, we can rewrite (37) as
ae
@
@~ae
ae
We can then change variables back to ae and add and subtract (1=dae)(@=@ae) "
to both sides to obtain
d
d
@ae
dae@ @
@ae
ae
Hence
H (ae) \Gammaae@ @
@ae
ae
@'A
and this becomes
@ae
ae
Collecting the equations (38), (39), (40) and (41) we obtain the following time
harmonic B'erenger system:
ae =ae
@ae
ae@ @
@ae
ae
@'A
@
@ae
RR n3049
22 Francis Collino et Peter Monk
Finally, we can move back to the time domain to obtain:
ae
@t
ae
@t
@ae
@t
@
@ae
ae
@H (ae)
@t
@ae
5 Numerical Results
5.1 A Special Case of the Time Harmonic Problem
In the special case of a concentric circular scatterer we can solve the standard
scattering problem and B'erenger problem using special functions. We will use
this case to examine the improvement gained by using the Robin problem over
the Dirichlet problem. Let us suppose that the
scatter\Omega is a circle of radius
a 0 .
First let us consider the standard scattering problem of computing "
H such
that (17), (18) and (19) are satisfied in the case when the
scatterer\Omega is a disk
of radius a 0 and the boundary data is
@ae
on
This problem has the solution
a
where
a
INRIA
The PML in Curvilinear Coordinates 23
If we truncate the domain using a circle of radius b, use a B'erenger layer in
use the Dirichlet condition on ae = b we must solve
r
@ae
@ae
Let use denote by ~
aej ae=b , then the above problem has the solution
a B
exp(in')
where
a B
When the problem is not uniquely solvable. Of course the solution
H B is independent of the details of the B'erenger layer and depends only on
~
ae b . This is not the case when numerical methods are used to discretize the
B'erenger layer since there will be spurious reflections from within the layer
and the layer must be designed to minimize these (as well as reach the desired
value of ~
ae b ).
We can evaluate the error in the B'erenger solution on the circle ae = a 0 by
computing the relative error:
In practice we truncate the series using \Gamma10  n  10. Arbitrarily we choose
We plot the error as a function of ~
ae in

Figure

1. When
~
ae is real (i.e. when there are eigenvalues for the annular domain and
RR n3049
Francis Collino et Peter Monk
these show up clearly as peaks in the error surface. As the imaginary part of ~
ae
increases the error decreases and can be made arbitrarily small if the imaginary
part is large enough. In Figure 2 we show the error as a function of the real
part of ~
ae for fixed imaginary part. By deciding on a desired error level, one
can find an appropriate value of ~
ae b .
In the subsequent sections we actually use the Robin boundary condition
on the outer boundary. In this case (for the circular scatterer of radius a 0 ) we
solve the following problem:
r
@ae
@ exp(ikx)
@ae
@ae
This problem has the solution
a B
exp(in')
where
a B
When the problem is not uniquely solvable. We could not prove existence
or uniqueness for this B'erenger problem in general but our calculations
suggest that for the annulus existence and uniqueness only fail at
the problem becomes a pure Neumann problem). This accounts for the error
spike towards ~
shown in Figure 3. Compared to Figure 1 we can see
that there are no other eigenvalues for the range of ~
ae considered here. The
error curves shown in Figure 4 (compared to Figure 2) show that the error
obtained using the Robin boundary condition is much smaller than the error
obtained using Dirichlet boundary conditions. This justifies our use of the
Robin condition in the numerical results presented later in this section.
INRIA
The PML in Curvilinear Coordinates 2513500.40.800.20.4
Real part
Imaginary part

Figure

1: Here we plot the relative L 2 error defined by (46) between the exact
solution of the time harmonic scattering problem for a circle and the B'erenger
solution using a Dirichlet outer boundary. We show the error as a function
of ~
ae b . When ~
ae b is real there are eigenvalues for the annular domain and these
show up as peaks in the error surface. As the imaginary part of ~
ae b increases,
the error decreases. Away from the real axis we see no evidence of lack of
existence or uniqueness for the B'erenger problem.
RR n3049
26 Francis Collino et Peter Monk
Real part
Imaginary
Real part

Figure

2: Here we plot the relative L 2 error defined by (46) between the exact
solution of the time harmonic scattering problem for a circle and the B'erenger
solution using a Dirichlet outer boundary. We show the error as a function of
In the top panel =(~ae b and in the bottom panel =(~ae b
the error is reduced to approximately 3%.
INRIA
The PML in Curvilinear Coordinates 2713500.40.800.20.4
Real part
Imaginary part

Figure

3: Here we plot the relative L 2 error defined by (46) between the exact
solution of the time harmonic scattering problem for a circle and the B'erenger
solution using a Robin outer boundary. We show the error as a function of ~
ae b .
For ~
small the error is large but for a sufficiently large real or imaginary part
the error can be made arbitrarily small. Apart from when
evidence of lack of existence or uniqueness for the B'erenger problem.
RR n3049
28 Francis Collino et Peter Monk
Real part
Imaginary
Imaginary

Figure

4: Here we plot the relative L 2 error defined by (46) between the exact
solution of the time harmonic scattering problem for a circle and the B'erenger
solution using a Robin outer boundary. We show the error as a function of
In the top panel =(~ae b and in the bottom panel =(~ae b
Compared to Figure 2 we see that the use of a Robin condition on the outer
boundary improves the error in the B'erenger solution. For example we can
reduce the error to approximately 3% by using =(~ae b
when we use the Dirichlet boundary condition.
INRIA
The PML in Curvilinear Coordinates 29
5.2 Time harmonic problems
The purpose of this section is to show that the Cartesian and cylindrical B'eren-
ger PML layers can be used to compute near field solutions of the Helmholtz
equation (and hence by a suitable transformation the far field). We shall only
show the results of one numerical experiment since we will provide a more
detailed numerical test of the curvilinear method for the time domain problem
(where the B'erenger method is most popular and where the behavior of the
Cartesian layer is well known).
For this test, we compute the field scattered from a perfectly conducting
metal obstacle in transverse polarization. This corresponds to Neumann data
imposed on the metal wall. The incoming wave is "
exp(ikx). The
scatterer is contained in the box \Gamma2
the speed of light to be unity so wave number is 6:2832. The grid
used is shown in Figure 5 where the outer radius of the circle is
the maximum diameter of elements in the mesh is 0:279. Note that
the grid is not aligned with Cartesian or radial coordinate system, which is
convenient from the point of view of using a standard mesh generator. Cubic
isoparametric elements are used to discretize the problem. For the Cartesian
case we use a simple Newumann boundary condition on the outer boundary
(it would be more usual in this case to use a rectangular outer boundary but
we want to make comparisons with a single grid). For the curvilinear PML
we use the modified Sommerfeld radiation condition being used on the outer
boundary (ie., we are solving the Robin problem (26), (27) and (31)). We
choose oe(s; In all cases we
choose oe
In the Cartesian case we use
oe(x; 2:5) and d y
This implies that the layer is half a distance unit from the scatterer. For the
curvilinear case we use

Figure

6 shows contours of the real part of "
computed by a capacitance
matrix technique that matches the finite element solution to an integral equation
solution outside the grid (thus handling the infinite domain accurately)
RR n3049
Francis Collino et Peter Monk
inlet.nopo
1373 POINTS
COIN
COIN
5.500 6.548

Figure

5: The mesh used in the numerical experiments reported on the time
harmonic wave equation.
INRIA
The PML in Curvilinear Coordinates 31
Matched FE/Nystrom
Real Part

Figure

The real part of the scattered field computed using a matched finite
element - integral equation code.
RR n3049
Francis Collino et Peter Monk
Berenger, Cartesian PML
Real Part
-0.577853

Figure

7: The real part of the scattered field computed using the Cartesian
B'erenger absorbing layer.
INRIA
The PML in Curvilinear Coordinates 33
Berenger, curvilinear PML
Real Part
-0.577828

Figure

8: The real part of the scattered field computed using the curvilinear
B'erenger absorbing layer.
RR n3049
34 Francis Collino et Peter Monk
[19].

Figure

7 shows the real part of the B'erenger solution "
using the Cartesian
layer, and 8 shows the real part of the B'erenger solution "
using the
curvilinear layer. Clearly the near field (for example in the scatterer) computed
using both B'erenger layers is similar to the capacitance matrix solution
in

Figure

6. But, as expected, the B'erenger solution dies away rapidly in the
absorbing layer. Furthermore, in both cases, the contours of the solution show
no abrupt curvature changes when the B'erenger layer is entered which indicate
that the layer is "perfectly" matched.

Figure

9 shows the real part of the scattered field along the x axis computed
using the matched finite element and integral equation approach. Figure
shows the corresponding result for the Cartesian and curvilinear
B'erenger media. The solutions are similar for jxj ! 2:5 (this is for coordinate
values from 2.5 to 7.5 on the one dimensional plots) but the B'erenger result
dies away quickly once the B'erenger medium is reached. These results give
numerical support to the claim that the rectilinear and curvilinear B'erenger
method can be used to compute time harmonic solutions.
5.3 Time Dependent System
In this section we shall investigate in some detail the curvilinear B'erenger
method for the time dependent problem. Figure 12 shows a schematic of the
grid used for the numerical tests of the curvilinear PML. Inside the circle of
a we use isoparametric quadrilateral N'ed'elec elements without
lumping [22]. Outside for a use isoparametric quadrilateral
N'ed'elec elements with an extra unknown associated with the centroid of the
elements. The curved segments of the mesh (for example in the B'erenger layer)
are approximated by piecewise linear approximations. This means that there
is some geometric error in approximating the B'erenger layers.
The arrangement of degrees of freedom for the elements is shown in Figure
13. The magnetic degrees of freedom (for H or H (ae) ) are associated with the
centroid of the elements and we use mapped piecewise constant basis functions
for the fields. The electric degrees of freedom for E ae and E ' are associated
with edges of the mesh.
INRIA
The PML in Curvilinear Coordinates 35
Matched FE/Nystrom
Magnitude
-0.005759950.5684321.142620.0 2.5 5.0 7.5 10.0

Figure

9: The real part of the scattered field computed using the matched
finite element - integral equation approach. Here we plot the solution along
the line on the x axis corresponds to the left outer boundary of
the domain). The break in the solution is due to the perfectly conducting
scatterer.
RR n3049
36 Francis Collino et Peter Monk
Berenger, Cartesian PML
Magnitude
-0.577853

Figure

10: The real part of the scattered field computed using the Cartesian
B'erenger absorbing medium. Here we plot the solution along the line
The break in the solution is due to the perfectly conducting scatterer.
INRIA
The PML in Curvilinear Coordinates 37
Berenger, Curvilinear PML
Magnitude
-0.577828
-0.004380230.5690681.142520.0 2.5 5.0 7.5 10.0

Figure

11: The real part of the scattered field computed using the curvilinear
B'erenger absorbing medium. Here we plot the solution along the line
The break in the solution is due to the perfectly conducting scatterer.
RR n3049
38 Francis Collino et Peter Monk
Circle r=a
Circle r=a'
Annular
Standard Maxwell
Zone
Zone
Beringer
Circle r=b

Figure

12: A schematic of the grid used. Inside the circle ae = a up to the scatterer
we use standard linear N'ed'elec elements on isoparametric quadrilaterals.
Outside four sided isoparametric B'erenger elements are used.
INRIA
The PML in Curvilinear Coordinates 39
(a) (b)

Figure

13: A schematic showing the arrangement of degrees of freedom for the
elements. The electric field unknowns are associated with edges in the mesh.
The magnetic field unknowns are associated with centroids of the elements.
(a) Standard N'ed'elec element, (b) B'erenger element.
On the outer edge of the domain (when we use the first order
absorbing boundary condition:
ffl. Note that this boundary condition is note the same as
the Sommerfeld absorbing boundary condition used in the frequency domain
calculations, but it is easier to implement for finite element methods.
Now let us detail the discrete problem. Let D b denote the bounded domain
including the B'erenger layer. Let  h be a mesh of linear isoparametric
quadrilaterals of maximum edge length h approximating the domain D b as
described above. Then the finite element space U h satisfies U h ae
and is constructed using the standard space of isoparametric linear edge ele-
ments. Let U \Gammag where \Gamma is the surface of the
scatterer (in the results shown later the boundary is the circle of radius 0.1).
h be the corresponding set of mapped piecewise constant functions containing
r \Theta U h . Let \Deltat ? 0 represent the time step (which must be sufficiently
small compared to h for stability). Then we
RR n3049
Francis Collino et Peter Monk
h and H 1=2
h approximate the fields
at respectively (the field H 1=2
h is computed from E 0
h and H 0
by a Runge-Kutta step applied to the Maxwell system) and H (ae);1=2
(we denote by (:; :) the L
inner product for vector or
scalar functions and use standard x and y components of the electric field):
\Deltat
ffloe

n \Theta
r \Theta H n+1=2
where (if
Next H (ae);n+1=2
computed via@ @ H (ae);n+3=2
@
@ae
';h denotes the ' component of E n+1
h . This computation is only
carried out in the B'erenger layer (the derivative of E n+1
';h is well defined by
virtue of having an annular grid in the B'erenger layer).
Finally H n+3=2
INRIA
The PML in Curvilinear Coordinates 41
Note that we have not used a centered difference for the B'erenger term here.
Numerical experiments show that this choice gives a better solution than the
standard symmetric average of B'erenger variables, but is likely to be method
dependent.
5.4 Scattering from a circle.
We consider an infinite domain problem of scattering of a plane wave off a
circular perfect conductor. We assume that the incident plane wave is given
by
where g is given by
This wave is incident on a perfectly conducting circular cylinder of radius 0.1m
centered at the origin. We choose
Using special function theory, an exact solution is available for this problem
[17].
For the numerical problem we take an annular domain with inner radius
0.1m and an outer radius of 0.3m. On
which describes how the incident field scatters off of a perfect conductor. The
grid used has mesh points distributed uniformly in the (r; ') plane, i.e. the
mesh points are (r
and
RR n3049
42 Francis Collino et Peter Monk

Figure

14: The base mesh used in the numerical experiments reported later
on the time dependent Maxwell system. Additional layers of cells are added
to form the B'erenger layer.
INRIA
The PML in Curvilinear Coordinates 43
A typical grid is shown in Figure 14 with N This is the
base mesh used for the calculations shown next.
The B'erenger boundary layer is implemented by adding extra layers of cells
outside the base mesh. We use a total of N layer layers. The radial size the these
cells is the 0:2=N r (ie. the same as for the cells in the base mesh). The time
step used is \Deltat = 0:01 and the mesh size (measured as the longest edge in
a quadrilateral) ranges from approximately near the scatterer to
on the outer layer of the base mesh. The solution is computed for
7.
In

Figure

15 we show a plot of H(x; against t at the point (x;
(\Gamma0:206; \Gamma0:0104). This is a grid point for the magnetic field. We can compute
a solution that is not degraded by mesh termination by computing on a grid
in which 150 extra layers have been added to the base mesh and in which the
B'erenger variables are set to zero. In this case the relative L 2 error is 0:18%.
Using the simple first order ABC (with 9 extra layers) gives a relative error
of 15:5%. We have tried three layers shown in Table 1. In each case oe is
computed by
where oe 0 is the value in the first B'erenger layer, i is the index of the layer, fl
is the growth factor, and h l is the radial mesh size in the layer.
From

Table

2 we can see that our time domain results produce solutions
that have less than a tenth of the error of the first order absorbing boundary
condition results. Unfortunately, with the time step of \Deltat = 0:01 the error in
the B'erenger solution is about four times the error in the underlying finite element
method. With a smaller time step this can be reduced to approximately
twice the underlying error suggesting that the phase error in the finite element
method is polluting the results.
6 Conclusion.
In this paper we have presented a method for implementing the B'erenger PML
in curvilinear coordinates. Our theory and numerical experiments suggest that
such a layer can be used in the frequency domain with a standard finite element
method. In the time domain the layer can also be used and gives more accurate
RR n3049
44 Francis Collino et Peter Monk
Time t
Magnetic
Field
Exact Solution
First Order Absorbing BC
Berenger Layer #1
Berenger Layer #2
Berenger Layer #3

Figure

15: A graph of H(x; against t at (x;
show the result of using a three different B'erenger layers with 9 shells in each
layer. The results are compared to an exact solution and the solution computed
using a first order absorbing boundary condition.
INRIA
The PML in Curvilinear Coordinates 45
Index of shell layer 1
9 0.2097152 0.4194304 0.8589935

Table

1: This table gives the value of h l oe for each ring in the three different
B'erenger layers used in the paper.
Thickness layer 1 layer 2 layer 3
6 1.03 0.92 (0.54) 1.03
9 0.74 (0.42) 0.92 (0.54) 0.80 (0.48)

Table

2: This table shows the relative L 2 error (as a percentage) in the solution
at layers with various thicknesses. The
structure of each layer is given in Table 1, and the thickness is the number
of rings (N layer )added to the base mesh. The finite difference error is 0.18%
(computed on a very large grid such that reflections from the outer boundary
do not pollute the solution). Numbers in parenthesis are for \Deltat = 0:001.
RR n3049
Francis Collino et Peter Monk
results than the first order absorbing condition. These results indicate that
it is worthwhile to test curvilinear B'erenger layers in three dimensional time
harmonic computations which we are currently doing.
Much remains to be done. The theoretical analysis of the B'erenger system
is incomplete. The presence (or absence hopefully) of exceptional frequencies
needs to be investigated further, and an error estimate for the layer needs to
be derived. Finally other discretization schemes (for example the curvilinear
method in [24]) could be investigated to find the one best suited to the B'erenger
layer.
INRIA
The PML in Curvilinear Coordinates 47



--R

A comparison of the Berenger perfectly matched layer and the Lindman higher-order ABC's for the FDTD method
A Perfectly Matched Layer for the Absorption of Electromagnetics Waves
Perfectly matched layer for the FDTD solution of wave-structure interaction problems
Proof of perfectly matched layer conditions in three dimensions
Modified Berenger PML absorbing boundary condition for FD-TD meshes
A 3D perfectly matched medium from modified Maxwell's equations with stretched coordinates

Perfectly matched layer for the paraxial equation.
Integral Equation Methods in Scattering Theory

Behavior of Berenger's ABC for evanescent waves
Implementation of Berenger absorbing boundary conditions in TLM by interfacing FDTD perfectly matched layers
Generalized perfectly matched layer-An extension of Berenger's perfectly matched layer boundary condition
Non reflecting boundary conditions
Modification of Berenger's perfectly matched layer for the absorption of electromagnetics waves in layered media
Initial boundary value problems for hyperbolic systems
The theory of electromagnetism
Validation and extension to three dimensions of the Berenger PML absorbing boundary condition for FD-TD meshes
An analysis of the coupling of finite element and Nystrom methods in acoustic scattering
Approximation of radiation boundary condi- tions

A comparison of three mixed methods for the time dependent Maxwell equations
Perfectly matched absorbing boundary conditions based on anisotropic lossy mapping of space

Numerical implementation and performance of perfectly matched layer boundary condition for waveguide structures

--TR

--CTR
Lin Zschiedrich , Roland Klose , Achim Schdle , Frank Schmidt, A new finite element realization of the perfectly matched layer method for Helmholtz scattering problems on polygonal domains in two dimensions, Journal of Computational and Applied Mathematics, v.188 n.1, p.12-32, 1 April 2006
Marcus J. Grote , Christoph Kirsch, Nonreflecting boundary condition for time-dependent multiple scattering, Journal of Computational Physics, v.221 n.1, p.41-62, January, 2007
Bjrn Sjgreen , N. Anders Petersson, Perfectly matched layers for Maxwell's equations in second order formulation, Journal of Computational Physics, v.209 n.1, p.19-46, 10 October 2005
Bradley Alpert , Leslie Greengard , Thomas Hagstrom, Nonreflecting boundary conditions for the time-dependent wave equation, Journal of Computational Physics, v.180 n.1, p.270-296, July 20, 2002
A. Bendali , Y. Boubendir , M. Fares, A FETI-like domain decomposition method for coupling finite elements and boundary elements in large-size problems of acoustic scattering, Computers and Structures, v.85 n.9, p.526-535, May, 2007
Frdric Nataf, A new approach to perfectly matched layers for the linearized Euler system, Journal of Computational Physics, v.214 n.2, p.757-772, 20 May 2006
Fang Q. Hu, A perfectly matched layer absorbing boundary condition for linearized Euler equations with a non-uniform mean flow, Journal of Computational Physics, v.208 n.2, p.469-492, 20 September 2005
A. Bermdez , L. Hervella-Nieto , A. Prieto , R. Rodrguez, An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems, Journal of Computational Physics, v.223 n.2, p.469-488, May, 2007
Adrianus T. de Hoop , Robert F. Remis , Peter M. van den Berg, The 3D wave equation and its Cartesian coordinate stretched perfectly matched embedding - A time-domain Green's function performance analysis, Journal of Computational Physics, v.221 n.1, p.88-105, January, 2007
Dan Givoli , Beny Neta, High-order non-reflecting boundary scheme for time-dependent waves, Journal of Computational Physics, v.186 n.1, p.24-46, 20 March
perfectly matched layer;finite element;absorbing layer